# Filter design by convex iteration

### From Wikimization

Line 34: | Line 34: | ||

- | Define <math> | + | Define <math>I_0\triangleq I\,</math> and define <math>I_n\,</math> as a zero matrix having vector <math>\mathbf{1}</math> along the <math>n^{\rm{th}}\,</math> superdiagonal when <math>n\,</math> is positive or <math>\mathbf{1}</math> along the <math>n^{\rm{th}}\,</math> subdiagonal when <math>n\,</math> is negative. |

By spectral factorization, <math>R(\omega)=|H(\omega)|^2\,</math> , | By spectral factorization, <math>R(\omega)=|H(\omega)|^2\,</math> , | ||

Line 46: | Line 46: | ||

& R(\omega)\leq\delta_2^2 & \omega\in[\omega_s\,,\pi]\\ | & R(\omega)\leq\delta_2^2 & \omega\in[\omega_s\,,\pi]\\ | ||

& R(\omega)\geq0 & \omega\in[0,\pi]\\ | & R(\omega)\geq0 & \omega\in[0,\pi]\\ | ||

- | & r(n)= | + | & r(n)=\textrm{trace}(I_n^{\rm T}G) &n=0\ldots N-1\\ |

- | + | ||

& G\succeq0\\ | & G\succeq0\\ | ||

& \textrm{rank}(G) = 1 | & \textrm{rank}(G) = 1 |

## Revision as of 22:47, 24 August 2010

where denotes impulse response.

For a low pass filter, frequency domain specifications are:

To minimize peak magnitude of , the problem becomes

But this problem statement is nonconvex.

So instead,

is a positive semidefinite rank 1 matrix.

Summing along each of subdiagonals gives entries of the autocorrelation function of .

In particular, the main diagonal of holds squared entries of .

Minimizing is therefore equivalent to minimizing .

Define and define as a zero matrix having vector along the superdiagonal when is positive or along the subdiagonal when is negative.

By spectral factorization, , an equivalent problem is:

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

But the above doesn't work.

Debug: try the following:

N = 16; delta1 = 1.01; delta2 = 0.01; w = linspace(0,pi,N); DFT = real(fft(eye(2*N))); cvx_begin variable r(2*N-1,1); r = [0; r]; R = fftshift(DFT*r); for I=1:N 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

The above just to see if there exists an such that the problem is feasible.

But it still doesn't work...