# Polyphase Filter decomposition. It is not working

So, I want to do a polyphase implementation of a filter bank (I had problems with it since the beginning and I already asked for help to the community in the past, indeed: Norm MPEG-1 Layer III (Mp3) Filter banks distortion). Thanks to the support I was able to make it work (without polyphase implementation).

Nonetheless, I tried to do all the analysis filters and synthesis filters by means of polyphase filters banks as follows:

From:

To:

The idea is, basically, based on the first picture, with the low pass filter prototype modulated by a cosine to make frequency shifts of the filter response to cover all the bandwidth (with 32 filters), decomposing each particular filter response into 32 polyphase branches.

(So 32 branches, obtained by means of shifting the prototype by multiplying a cosine. Each of these 32 branches is implemented by means of decomposing each analysis and synthesis filters into 32 polyphase branches).

So, once the first part works (i.e. The filter bank not implemented by means of polyphase filter banks), I do the polyphase decomposition as follows:

1. Polyphase decomposition:

• h are the modulated coefficients of each filter of the bank.
• hh is the polyphase filter bank coefficients
• the $i^{\rm th}$ row of the matrix represents the coefficients of the $i^{\rm th}$ filter of the polyphase bank
• M are the polyphase filter branches, 32 in our example
hh=zeros(M,length(h)/M);
for l=0:1:M-1
n=l+1:M:length(h);
hh(l+1,:)=h(n);
end
end

2. The shift ($Z^{-1}$) is performed by a function that performs as follows ("delay" is the delay to apply. e.g. if 0 samples are to be delayed, then the output shall be equal to the input):

• in_data: Input sequence
• out_data: Output sequence.
delay_samp=delay+1;
out_data=zeros(1,length(in_data));
out_data(delay_samp:length(in_data))=in_data(1:length(in_data)-(delay_samp-1));

3. Ok, now, I apply the filters as follows:

Analysis:

yy=zeros(branches,length(hh(1,:))+fix(length(in_signal)/branches));

for k=1:branches

%1. Shift
shifted_input=shift_data(in_signal,k-1);

%2. Downsample
temp=zeros(1,fix(length(in_signal)/branches));
downsampled_data=shifted_input(1:branches:end);

temp(1:length(downsampled_data))=downsampled_data;
shifted_downsampled_input=temp;

%3. Filter
temp=conv(hh(k,:),shifted_downsampled_input);
yy(k,:)=temp;

end

And Synthesis:

yy=zeros(branches,branches*(length(hh(1,:))+length(in_signal)-1));

for k=1:branches

%1. Filter
temp_conv=conv(hh(k,:),in_signal);

%2. Upsample
temp=zeros(1,length(temp_conv)*branches);
temp(1:branches:end)=temp_conv;%upsample(temp,branches);

%3. Shift.
yy(k,:)=shift_data(temp,branches-(k-1));

end

The result is very distorted:

The output with a 1KHz tone as an input:

I guess that I am making something wrong with the upsampling and downsampling processes as there are spectral replicas all over the second example, but I have no clue on what exactly is wrong.

I checked on the internet and I also checked the reference referred in the other post. I mean this book. But I think there is a problem with the way I am coding this not with the theoretical back.

• ...Let me thank you for your huge help on this one and on the other question (and also for the tips you gave me on using "conv" instead of "filter" and the reference book you referred me to). So: Thanks a lot!! Commented Dec 8, 2017 at 17:34

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)))));

• Code above uses MATLAB R12 (from year 2000). On a more recent version, you should modify 3D arrays hhap and hhsp accordingly. Commented Jun 12, 2023 at 17:31