4
$\begingroup$

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: enter image description here

To: enter image description here

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: enter image description here

The output with a 1KHz tone as an input: enter image description here

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.

This is the low pass filter prototype.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Fat32 Dec 7 '17 at 22:12
  • $\begingroup$ It's working now, isn't it? :-))) $\endgroup$ – Fat32 Dec 8 '17 at 15:51
  • 1
    $\begingroup$ 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. ... $\endgroup$ – f.gallardo Dec 8 '17 at 17:34
  • 1
    $\begingroup$ ...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!! $\endgroup$ – f.gallardo Dec 8 '17 at 17:34
  • $\begingroup$ 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(). $\endgroup$ – Fat32 Dec 8 '17 at 19:08
2
$\begingroup$

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)))));
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Fat32 Dec 8 '17 at 21:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.