Filter design by convex iteration

From Wikimization

(Difference between revisions)
Jump to: navigation, search
Line 76: Line 76:
<pre>
<pre>
-
N = 32;
+
N = 16;
delta1 = 1.01;
delta1 = 1.01;
delta2 = 0.01;
delta2 = 0.01;
w = linspace(0,pi,N);
w = linspace(0,pi,N);
 +
DFT = real(fft(eye(2*N))); % generate DFT matrix but only use the real components
 +
% https://ccrma.stanford.edu/~jos/st/DFT_Matrix.html
cvx_begin
cvx_begin
-
variable h(N,1);
 
variable r(2*N-1,1);
variable r(2*N-1,1);
-
variable g(N^2,1);
 
variable G(N^2,N^2);
variable G(N^2,N^2);
-
variable R(2*N-1,1);
+
 
 +
r = [0; r]; % otherwise r is not even
 +
-
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)));
minimize(max(diag(G)));
 +
 +
G2 = G(1:N,1:N); % I don't really need to use the entire G.
 +
% So just use the 1st block: G2
for I=1:N
for I=1:N
-
r(I) = 1/I*trace(diag(ones(1,N),I-N)*G);
+
temp = diag(ones(1,I),-(N-I));
 +
[i,j,k] = find(temp);
 +
index = [i j ones(length(i),1)];
 +
index = [index; [N, N, 0]];
 +
temp2 = spconvert(index);
 +
r(I) == trace(temp2*G2);
end
end
 +
for I=N+1:2*N-1
for I=N+1:2*N-1
-
r(I) = 1/(I-N)*trace(diag(ones(1,N),I-N)*G);
+
temp = diag(ones(1,2*N-I),I-N);
 +
[i,j,k] = find(temp);
 +
index = [i j ones(length(i),1)];
 +
index = [index; [N, N, 0]];
 +
temp2 = spconvert(index);
 +
r(I) == trace(temp2*G2);
end
end
-
for I=1:N
+
-
temp = 0;
+
R = fftshift(DFT*r);
-
for J=1:N-1
+
-
temp = temp + 2*r(J)*cos(w(J)*J);
+
for I=1:N % R is of length of 2N
-
end
+
-
R(I) = r(1) + temp;
+
-
end
+
-
for I=1:N/2
+
1/delta1^2 < R(I);
1/delta1^2 < R(I);
R(I) < delta1^2;
R(I) < delta1^2;
end
end
-
for I=N/2+1:N
+
for I=N+1:2*N
R(I) < delta2^2;
R(I) < delta2^2;
end
end

Revision as of 12:23, 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.

So instead, a new vector

LaTeX: g \triangleq \left[ </p> <pre> \begin{array}{c} h(n) \\ h(n-1) \\ \vdots \\ h(n-N) \\ \end{array} \right]\in\mathbb{C}^{N^2} </pre> <p>

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 = 16;
delta1 = 1.01;
delta2 = 0.01;
w = linspace(0,pi,N);
DFT = real(fft(eye(2*N)));  % generate DFT matrix but only use the real components
                            % https://ccrma.stanford.edu/~jos/st/DFT_Matrix.html

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

  r = [0; r];   % otherwise r is not even
  
  
  minimize(max(diag(G)));

  G2 = G(1:N,1:N);       % I don't really need to use the entire G.  
                         % So just use the 1st block: G2
  for I=1:N
    temp = diag(ones(1,I),-(N-I));
    [i,j,k] = find(temp);
    index = [i j ones(length(i),1)];
    index = [index; [N, N, 0]];
    temp2 = spconvert(index);
    r(I) == trace(temp2*G2);
  end
   
  for I=N+1:2*N-1
    temp = diag(ones(1,2*N-I),I-N);
    [i,j,k] = find(temp);
    index = [i j ones(length(i),1)];
    index = [index; [N, N, 0]];
    temp2 = spconvert(index);
    r(I) == trace(temp2*G2);
  end
  
  R = fftshift(DFT*r);
  
  for I=1:N                 % R is of length of 2N
      1/delta1^2 < R(I);
      R(I) < delta1^2;
  end
  for I=N+1:2*N
      R(I) < delta2^2;
  end

  R > 0;
  
  
cvx_end
Retrieved from "http://www.convexoptimization.com/wikimization/index.php/Filter_design_by_convex_iteration"
Personal tools