% Polyphase filter implementation

% Determine the filters.
[h0,h1,f0,f1] = orthfilt(dbwavf('db2'));
M = length(h0);

% Load a signal.
load noisdopp
x = noisdopp;
N = length(x);

% Change the signal into polyphase form.
xeven = dyaddown(x,1);  % Even part
xodd = dyaddown(x,0);   % Odd part
xodd = [0 xodd(1:N/2-1)];
X = [xeven; xodd];

% Construct the polyphase matrix.
H = zeros(2,2,M/2);
H(1,1,:) = dyaddown(h0,1);  % h0,even[n]
H(1,2,:) = dyaddown(h0,0);  % h0,odd[n]
H(2,1,:) = dyaddown(h1,1);  % h1,even[n]
H(2,2,:) = dyaddown(h1,0);  % h1,odd[n]

H(:,:,1)
H(:,:,2)

% Run the polyphase filter.
Y = polyfilt(H,X);

% Plot the results.
L = N/2;
n = 0:L-1;
clf
subplot(2,1,1)
plot(n,Y(1,1:L))
axis tight
xlabel('Sample number')
ylabel('Lowpass')
title('Output from polyphase filter')
subplot(2,1,2)
plot(n,Y(2,1:L))
axis tight
xlabel('Sample number')
ylabel('Highpass')
pause

% Compute the results using the direct approach.
y0 = dyaddown(conv(x,h0),1);
y1 = dyaddown(conv(x,h1),1);

% Now compare the results.
clf
subplot(2,1,1)
plot(n,Y(1,1:L)-y0(1:L))
axis tight
xlabel('Sample number')
ylabel('Lowpass difference')
title('Difference in outputs produced by polyphase and direct forms')
subplot(2,1,2)
plot(n,Y(2,1:L)-y1(1:L))
axis tight
xlabel('Sample number')
ylabel('Highpass difference')
pause

% Plot the determinant of the polyphase matrix as a function of frequency.
R = 32;
W = 2/R*(-R/2:R/2-1);
H0even = fftshift(fft(H(1,1,:),R));
H0odd = fftshift(fft(H(1,2,:),R));
H1even = fftshift(fft(H(2,1,:),R));
H1odd = fftshift(fft(H(2,2,:),R));
delta = zeros(1,R);
delta(:) = H0even .* H1odd - H0odd .* H1even;
clf
plot(W,abs(delta),'x-')
axis([-1 1 0 1.5])
xlabel('Angular frequency (normalized by pi)')
ylabel('Magnitude of determinant')
title('Determinant of the polyphase matrix')
pause

% Verify that the filter is orthogonal i.e. Hp'(w*) Hp(w) = I
A11 = zeros(1,R);
A11(:) = abs(H0even).^2 + abs(H1even).^2;
A12 = zeros(1,R);
A12(:) = conj(H0even).*H0odd + conj(H1even).*H1odd;
A21 = zeros(1,R);
A21(:) = conj(H0odd).*H0even + conj(H1odd).*H1even;
A22 = zeros(1,R);
A22(:) = abs(H0odd).^2 + abs(H1odd).^2;
clf
subplot(4,1,1)
plot(W,A11,'x-');
axis([-1 1 0 1.5])
ylabel('A11')
title('Variation of A = Hp''(w*) Hp(w) with frequency')
subplot(4,1,2)
plot(W,A12,'x-')
axis([-1 1 0 1.5])
ylabel('A12')
subplot(4,1,3)
plot(W,A21,'x-')
axis([-1 1 0 1.5])
ylabel('A21')
subplot(4,1,4)
plot(W,A22,'x-')
axis([-1 1 0 1.5])
xlabel('Angular frequency (normalized by pi)')
ylabel('A22')