Matlab Exercise #5 - Spectral factorization

In this problem we will consider a practical example of spectral factorization. The following Matlab function is meant to return the Daubechies minimum-phase filter D_k for any k, following the procedure illustrated on pages 130-131 of the textbook. The code provided below computes the coefficients of the polynomial R(z), following the notation of the book. (Note: you can dowload the code at the usual http address).

(a) Complete the daub function by implementing a spectral factorization algorithm for R(z) and by computing the corresponding Daubechies filter's coefficients. Compare your results for k=2 and k=4 to the filter values given on page 131 and 129 respectively; they should be identical. (Hint: consider the orthonormality requirement when it comes to normalizing the coefficients.)

(b) The proposed technique, which is based on finding the roots of a polynomial, can become numerically unstable for large values of k. Matlab warns you of that: which is the first value of k for which Matlab prints a warning message ?

You can verify yourself the consequences of this numerical problem: plot, superimposed, the values of $\vert\langle d_{k}[n],d_{k}[n+2m]\rangle\vert$for m=1 to 5 and k=8 to 18, where d_k[n] is the D_k Daubechies filter. Provide comments on the plot.

function f = daub(k)

% Length of R(z), which is also the order of the filter
r = 2*k-1;

% Let's compute the coefficients of (1 + z^-1)^k * (1 + z)^k
% Remember that polynomial multiplication is equivalent to convolution

p = [1 1]; 	%  This can represent both (1 + z^-1) and (1 + z)

for n=1:k-1
	p = conv(p, [1 1]);
end
pz = conv(p, p);

% Now we want to build a convolution matrix A so that AR , with 
%  R = [r(k) r(k-1) ... r(0) ... r(k)]', gives us the coefficients of
%  the even powers of P(z).

% Time-reverse (not necessary due to symmetry) and zero-pad pz
pz = [zeros(1,r) pz zeros(1,r)];
l = length(pz)

A = zeros(r,r);
for n=1:r
	A(n,:) = pz(l-2*n-r+1:l-2*n);
end

% Now we solve AR = B, where B is zero for all odd powers except for z^0

B = zeros(r,1); B((r+1)/2) = 1;
R = A\B;