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 optimizations 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 (complementary) filter bank : % ----------------------------------------------------------------- numbands = 32; % number of banks n=0:L-1; hhs = zeros(numbands,L); % bank 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 = interpolation 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: ANALYSIS STAGE : % ----------------------------------------------------------- yyd = zeros( numbands, floor(N/numbands)); % decimated outputs.. M = numbands; for k=1:1:numbands 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 SYNTHESIS filterbank 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)))));