# Filter design by convex iteration

(Difference between revisions)
 Revision as of 04:03, 24 August 2010 (edit)← Previous diff Revision as of 04:46, 24 August 2010 (edit) (undo)Next diff → Line 77: Line 77:

N = 32;                                                                                                    N = 32;
-                                                           delta1 = 0.01;                              +                                                              delta1 = 1.01;
delta2 = 0.01;                                                                                             delta2 = 0.01;
+                                                              w = linspace(0,pi,N);

-                                                           cvx_begin                                   +                                                              cvx_begin
-                                                           variable h(N,1);                            +                                                              variable h(N,1);
-                                                                                                       +                                                              variable r(2*N-1,1);
-                                                           G = h*h';                                   +                                                              variable g(N^2,1);
-                                                           minimize(norm(diag(G),inf));                +                                                              variable G(N^2,N^2);
-                                                           r = xcorr(h);                               +                                                              variable R(2*N-1,1);
-                                                           R = fftcc(r);                               +
-                                                           1/delta1^2 <= R <= delta1^2;    %%% range?  +                                                              for I=1:N
+                                                              g((I-1)*N+1:I*N,1) = circshift(h,[0,I-1]);
+                                                              end
+                                                              G = g*g';
+                                                              %   G >= 0;
+                                                              minimize(max(diag(G)));
+                                                              for I=1:N
+                                                              r(I) = 1/I*trace(diag(ones(1,N),I-N)*G);
+                                                              end
+                                                              for I=N+1:2*N-1
+                                                              r(I) = 1/(I-N)*trace(diag(ones(1,N),I-N)*G);
+                                                              end
+                                                              for I=1:N
+                                                              temp = 0;
+                                                              for J=1:N-1
+                                                              temp = temp + 2*r(J)*cos(w(J)*J);
+                                                              end
+                                                              R(I) = r(1) + temp;
+                                                              end
+                                                              for I=1:N/2
+                                                              1/delta1^2 < R(I);
+                                                              R(I) < delta1^2;
+                                                              end
+                                                              for I=N/2+1:N
+                                                              R(I) < delta2^2;
+                                                              end
+
+                                                              R > 0;
+
+
cvx_end                                                                                                    cvx_end

## Revision as of 04:46, 24 August 2010

$LaTeX: H(\omega) = h(0) + h(1)e^{-j\omega} + \cdots + h(N-1)e^{-j(N-1)\omega}$  where  $LaTeX: h(n)\in\mathbb{C}^N$ denotes impulse response.

For a low pass filter, frequency domain specifications are:

$LaTeX: \begin{array}{ll} \frac{1}{\delta_1}\leq|H(\omega)|\leq\delta_1\,, &\omega\in[0,\omega_p]\\ |H(\omega)|\leq\delta_2\,, &\omega\in[\omega_s\,,\pi] \end{array}$

To minimize peak magnitude of $LaTeX: h\,$ , the problem becomes

$LaTeX: \begin{array}{cll} \textrm{minimize} &\|h\|_\infty\\ \textrm{subject~to} &\frac{1}{\delta_1}\leq|H(\omega)|\leq\delta_1\,, &\omega\in[0,\omega_p]\\ &|H(\omega)|\leq\delta_2\,, &\omega\in[\omega_s\,,\pi] \end{array}$

But this problem statement is nonconvex.

$LaTeX: g \triangleq \left[

 \begin{array}{c} h(n) \\ h(n-1) \\ \vdots \\ h(n-N) \\ \end{array} \right]\in\mathbb{C}^{N^2}

$

is defined by concatenation of time-shifted versions of $LaTeX: h\,$ .

Then

$LaTeX: G\triangleq gg^\mathrm{H}=\,\left[\begin{array}{*{20}c} h(n)h(n)^{\rm H} & h(n)h(n-1)^{\rm H} & h(n)h(n-2)^{\rm H} & \cdots & h(n)h(n-N)^{\rm H}\\ h(n-1)h(n)^{\rm H} & h(n-1)h(n-1)^{\rm H} & h(n-1)h(n-2)^{\rm H} & \cdots & h(n-1)h(n-N)^{\rm H}\\ h(n-2)h(n)^{\rm H} & h(n-2)h(n-1)^{\rm H} & h(n-2)h(n-2)^{\rm H} & \ddots & h(n-2)h(n-N)^{\rm H}\\ \vdots & \vdots & \ddots & \ddots & \vdots\\ h(n-N)h(n)^{\rm H} & h(n-N)h(n-1)^{\rm H} & h(n-N)h(n-2)^{\rm H} & \cdots & h(n-N)h(n-N)^{\rm H} \end{array}\right]\in\mathbb{C}^{N^2\times N^2}$

is a positive semidefinite rank 1 matrix.

Summing along each of $LaTeX: 2N-1\,$ subdiagonals gives entries of the autocorrelation function $LaTeX: r\,$ of $LaTeX: h\,$ .

In particular, the main diagonal of $LaTeX: G\,$ holds squared entries of $LaTeX: h\,$ .

Minimizing $LaTeX: \|h\|_\infty$ is therefore equivalent to minimizing $LaTeX: \|\textrm{diag}(G)\|_\infty$ .

Define $LaTeX: I_n\triangleq\,$...

By spectral factorization, $LaTeX: R(\omega)=|H(\omega)|^2\,$ , an equivalent problem is:

$LaTeX: \begin{array}{cll} \textrm{minimize}_{G\,,\,r}&\|\textrm{diag}(G)\|_\infty\\ \textrm{subject~to} & R(\omega) = r(0)+\!\sum\limits_{n=1}^{N-1}2r(n)\cos(\omega n)\\ &\frac{1}{\delta_1^2}\leq R(\omega)\leq\delta_1^2 &\omega\in[0,\omega_p]\\ & R(\omega)\leq\delta_2^2 & \omega\in[\omega_s\,,\pi]\\ & R(\omega)\geq0 & \omega\in[0,\pi]\\ & r(n)=\frac{1}{n}\textrm{trace}(I_{n-N}G) &n=1\ldots N\\ & r(n)=\frac{1}{n-N}\textrm{trace}(I_{n-N}G) &n=N+1\ldots 2N-1\\ & G\succeq0\\ & \textrm{rank}(G) = 1 \end{array}$

Excepting the rank constraint, this problem is convex.

N = 32;
delta1 = 1.01;
delta2 = 0.01;
w = linspace(0,pi,N);

cvx_begin
variable h(N,1);
variable r(2*N-1,1);
variable g(N^2,1);
variable G(N^2,N^2);
variable R(2*N-1,1);

for I=1:N
g((I-1)*N+1:I*N,1) = circshift(h,[0,I-1]);
end
G = g*g';
%   G >= 0;
minimize(max(diag(G)));
for I=1:N
r(I) = 1/I*trace(diag(ones(1,N),I-N)*G);
end
for I=N+1:2*N-1
r(I) = 1/(I-N)*trace(diag(ones(1,N),I-N)*G);
end
for I=1:N
temp = 0;
for J=1:N-1
temp = temp + 2*r(J)*cos(w(J)*J);
end
R(I) = r(1) + temp;
end
for I=1:N/2
1/delta1^2 < R(I);
R(I) < delta1^2;
end
for I=N/2+1:N
R(I) < delta2^2;
end

R > 0;

cvx_end