The following is a working code that uses 32-component polyphase decomposition of the associated 32-channel anslysis and synthesis filterbanks. As I have already commented, the speed gain is not dramatic in this cae due to short signal and filter lengths. However further architectural improvements as well as coding optimizaitonoptimizations can provide better results.
% S0 - Load the prototype lowpass filter impulse response h0[n]:
% --------------------------------------------------------------
load h2.mat; % h[n] is the prototype lowpass filter of length 512
L = length(h);
% S1 - Create the 32 x 512 analysis filter-bank hha[k,n] by cosine modulation from protoype :
% -----------------------------------------------------------------------------------------
numbands = 32; % number of banks (channels)
n=0:L-1;
hha=zeros(numbands,L); % bank of filters hha[k,n] = 32 x 512 array.
for k=0:1:numbands-1
hha(k+1,:) = h.*cos( ( (2*k+1)*pi*(n-16) ) / (2*numbands) );
end
% S2 - Create the 32-polyphase components hhap[k,m,n] , for each one of 32 analysis filters hha[k,n]:
% ---------------------------------------------------------------------------------------------------
numpoly = numbands; % polyphase component number = decimation ratio = number of channels
hhap = zeros(numbands,numpoly, L/numpoly); % hhap = 32 x 32 x 512/32 , 3D ANALYSIS filter bank array
M = numpoly; % polyphase system decimation ratio
for k=1:numbands
for m = 1:numpoly
hhap(k,m,:) = hha(k,m:M:end); % create the m-th polyphase component of k-th channel filter
end
end
% S3 - Design the 32 x 512 synthesis (cmplementarycomplementary) filter bank :
% -----------------------------------------------------------------
numbands = 32; % number of banks
n=0:L-1;
hhs = zeros(numbands,L); % bankdbank of filters
for k=0:1:numbands-1
hhs(k+1,:) = h.*cos( ( (2*k+1)*pi*(n+16) ) / (2*numbands) );
end
% S4 - Obtain the 32-polyphase components hhsp[k,m,n] , for each one of 32 synthesis filters hhs[k,n]:
% ----------------------------------------------------------------------------------------------------
numpoly = numbands; % polyphase component number = decimationinterpolation ratio = number of channels
hhsp = zeros(numbands,numpoly, L/numpoly); % hhap = 32 x 32 x 512/32 , 3D ANALYSIS filter bank array
M = numpoly; % polyphase system decimation ratio
for k=1:numbands
for m = 1:numpoly
hhsp(k,m,:) = hhs(k,m:M:end); % create the m-th polyphase component of k-th channel filter
end
end
% S5 - Generate the test input signal
% -----------------------------------
N = 2*1024;
wav_in = cos(0.01791*pi*[0:N-1]); % pure sine tone
% S6 - Apply test signal to the filterbank, using the: filterbandANALYSIS approachSTAGE :
% --------------------------------------------------------------------------
yyd = zeros( numbands, floor(N/numbands)); % decimated outputs..
M = numbands;
for k=1:1:numbands
%temp = conv( wav_in,hha(k,:)); % THESE STEPS are now implemented as a POLYPHASE filter
temp = conv([wav_in(1:M:end),0] , hhap(k,1,:));
for m=2:M
temp = temp + conv([0,wav_in(M-m+2:M:end)],hhap(k,m,:));
end
yyd(k,:) = temp(L/(2*M)+1 : L/(2*M)+N/numbands);
end
% S7 - Apply synthesisSYNTHESIS filterbanksfilterbank on the decimated signal :
% ----------------------------------------------------------
ys = zeros(1, N);
for k=1:numbands
temp = zeros(1, N+L-1);
for m = 1:numpoly
temp(m:numbands:end-31) = conv( yyd(k,:) , hhsp(k,m,:) );
end
ys = ys + temp(L/2+1:L/2+N);
end
ys = numbands*ys;
% SX - DISPLAY RESULTS:
% ---------------------
L = length(h);
figure,subplot(2,1,1)
stem([0:L-1],h);title('The Prototype Lowpass Filter');
subplot(2,1,2)
plot(linspace(-1,1,4*L),20*log10(abs(fftshift(fft(h,4*L)))));
grid on;
figure
plot(linspace(-1,1,4*L),20*log10(abs(fftshift(fft(hha(1,:),4*L)))));
hold on
for k=2:numbands
plot(linspace(-1,1,4*L),20*log10(abs(fftshift(fft(hha(k,:),4*L)))));
end
title('32 CHANNEL FILTERBANK');
figure,subplot(2,1,1)
plot(wav_in);title('input signal')
subplot(2,1,2)
plot(linspace(-1,1,4*N),20*log10(abs(fftshift(fft(wav_in,4*N)))));
figure,subplot(2,1,1)
plot(ys);title('Synthesized Back');
subplot(2,1,2)
plot(linspace(-1,1,4*N),20*log10(abs(fftshift(fft(ys,4*N)))));