Filter design by convex iteration
From Wikimization
(Difference between revisions)
Line 77: | Line 77: | ||
<pre> | <pre> | ||
N = 32; | N = 32; | ||
- | delta1 = | + | 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 = | + | variable g(N^2,1); |
- | minimize( | + | variable G(N^2,N^2); |
- | r = | + | variable R(2*N-1,1); |
- | + | ||
- | 1/delta1^2 < | + | 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 | ||
</pre> | </pre> |
Revision as of 05:46, 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 = 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