# Filter design by convex iteration

$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 h(n)h(n)^{\rm H} \in\mathbb{C}^{N\times N}$

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_0\triangleq I\,$ and define $LaTeX: I_n\,$ as a zero matrix having vector $LaTeX: \mathbf{1}$ along the $LaTeX: n^{\rm{th}}\,$ superdiagonal when $LaTeX: n\,$ is positive or $LaTeX: \mathbf{1}$ along the $LaTeX: n^{\rm{th}}\,$ subdiagonal when $LaTeX: n\,$ is negative.

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)=\textrm{trace}(I_n^{\rm T}G) &n=0\ldots N-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


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 $LaTeX: r$ such that the problem is feasible.

But it still doesn't work...