# Filter design by convex iteration

### From Wikimization

(Difference between revisions)

Line 76: | Line 76: | ||

<pre> | <pre> | ||

- | N = | + | 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); | ||

- | + | ||

+ | 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 | ||

- | + | 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 | ||

- | + | 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 | ||

- | + | ||

- | + | R = fftshift(DFT*r); | |

- | + | ||

- | + | for I=1:N % R is of length of 2N | |

- | + | ||

- | + | ||

- | + | ||

- | for I=1:N | + | |

1/delta1^2 < R(I); | 1/delta1^2 < R(I); | ||

R(I) < delta1^2; | R(I) < delta1^2; | ||

end | end | ||

- | for I=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

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, a new vector

is defined by concatenation of time-shifted versions of .

Then

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 ...

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