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:

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.

• So you are doing this for computatinal efficiency purposes. But a naive software implemenation may not provide whay you looked for. Indeed an Octave implementation would be even slower, with such lengths of 1024 signal samples, 512 filter samples and 32 polyphase components, due to overheads? Dec 7 '17 at 22:12
• It's working now, isn't it? :-))) Dec 8 '17 at 15:51
• Hell yes! :) Thanks! On the other hand (shitty comment on my side ;-P ), it seems to like the "conv" function changed or something. An error was raised (both in Matlab and Octave) about the size of the matrices. It seems like there was a problem with the matrix where the coefficients were stored. I sorted out that by strictly assigning the values of the coefficients of interest to an empty array, that way Matlab/Octave "conv" did not complain about that. Still, I don't understand what I was doing wrong in my code but I'll use your example to understand my problem. ... Dec 8 '17 at 17:34
• ...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!! Dec 8 '17 at 17:34
• I haven't checked your code so I cannot tell which part of it caused the error but you can take my implementation as a basis and develop your own version. About the conv-filter functions: indeed they are the same (they produce the same outputs) except the fact that conv() computes the full length ($L_y = L_x + L_h-1$) result of the mathematical conv operator $y[n] = x[n] * h[n]$, whereas the filter(), computes solution of the LCCDE representation of an LTI system, and produces only the first $L_x$ samples of the output...Once you are ok with the code, you can replace conv() with filter(). Dec 8 '17 at 19:08

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)))));
• NOTICE: I've developed the code on a machine with Matlab 6.0 (yes it's) and it works ok (that's why I posted it) but now I just chekced in on a recent version (2015) and I see that the 3D array hhap(32,32,16) and hhsp(32,32,16) causes a problem in the conv call, so the code should be modified accordingly. The necessary modification will not be big task though, it can effect the performance. Nevertheless the true performance of a polyphase filter can better be oberseved on a pure C/C++ kind of implementation. Dec 8 '17 at 21:20