2
$\begingroup$

So, I'm working with the filter banks used in Mp3, in Matlab. I am obtaining quite a huge distortion in the output and I don't really know whether I'm making something wrong.

I've reviewed the literature and I'd say I'm making it properly, so I don't really understand what I'm making wrong.

So I'm basing the work on a filter with 512 coefficients that can be found on page 9-13, here: https://www.mp3-tech.org/programmer/docs/jacaba_main.pdf (I think it comes from a standard so no big deal there, I guess).

I want to make a filter bank with 32 branches. I want to, based on that Low Pass filter prototype, calculate the 32 branches that will cover the entire bandwidth, so I make:

(note "h" is the variable with the filter coefficients)

branches=32;
hh=zeros(branches-1,length(h));
n=0:511;
for k=0:1:branches-1 
   hh(k+1,:)=h.*cos(((2*k+1)*pi*(n-16))/(2*32));
end

That one for the generation of the analysis filters coefficients and the following one for the synthesis filters:

branches=32;
hh=zeros(branches-1,length(h));
n=0:511;
for k=0:1:branches-1 
   hh(k+1,:)=h.*cos(((2*k+1)*pi*(n+16))/(2*32));
end

As per the way of applying the filters to the input signal and to reconstruct the output signal:

yy=zeros(branches,floor(length(wav_in)/32)+1)
for k=1:1:branches
  yy(k,:)=decimate(filter(hh(k,:),1,wav_in),branches);
end

This will provide a matrix with all the outputs of the shifted analysis filters (yy).

Now, in order to reconstruct the signal, with the synthesis filters: (Note it's a different function so we will not have problems with the variables names)

yy = zeros(1,length(xd(1,:))*32)
for k=1:1:branches
  yy=yy+(filter(hh(k,:),1,upsample(xd(k,:),branches)));
end

The filter responses look good (it seems like that to me, at least): Analysis Filter bank response

Synthesis Filter bank reponse

If I use as an input a 1 kHz tone, I have:

enter image description here

The output is "After". As you can see the signal is very distorted. In general, I am following the approach in "How sound is processed in MP3 player?" by T. Dutoit and N. Moreau.

I think it is not working and honestly, I don't have a single clue on what I'm doing wrong. (Maybe how I'm applying the filter banks?).

$\endgroup$
  • $\begingroup$ Can you please also post the low-pass protoype filter impulse response coefficients ? As it's hard to copy from the pdf document and i don't want to use a homemade filter there... $\endgroup$ – Fat32 Oct 29 '17 at 13:53
  • $\begingroup$ Hi! Sure, you can grab the mat file from here: drive.google.com/open?id=0B3aZza3bav9gWWwwSEpYNGgxTE0 $\endgroup$ – f.gallardo Oct 29 '17 at 15:35
  • $\begingroup$ Can you edit your source code to omit redundant/confusing parentheses? cos(((((2*k+1)*pi*(n-16))/(2*32)))) is plain bad code. Same goes for things like doing a convolution of complicated expressions just to get its length. The convolution of two sequences has a defined length, and there's no reason to actually first upsample something to a defined length, and then convolve it with something of a defined length, just to calculate length() of that. I'll have to ask you to make your code acceptably well-written before asking about it. $\endgroup$ – Marcus Müller Oct 29 '17 at 16:41
  • $\begingroup$ Thanks for your comment @MarcusMüller. As per the coise content , I'll update it. As per the length of the convolution. You are right, I was kind of lazy there. BTW: Although the length of a conv. of two sequences has a defined length (the sum of both lengths -1), filter function makes length(out)=length(in). $\endgroup$ – f.gallardo Oct 29 '17 at 17:10
  • $\begingroup$ So, this can make the length line easier: yy=zeros(1,length(xd(1,:))*32) Anyway, you are right: the code shall always be as readable as possible. $\endgroup$ – f.gallardo Oct 29 '17 at 17:11
2
$\begingroup$

I have tested your code.

There are problems with the decimation and upsample stages. When you run the analysis and synthesis filters at the undecimated rate, then you get the perfect reconstruction (plus a filter shift)

First of all in case of applications which are sensitive to group delay, you should better use the conv / filter function and then compansate for the group delay.

Also, the decimate function applies a lowpass anti-aliasing filter, which kills the signal content at the $k$-th channel considering the fact that its output is a bandpass signal rather than lowpass. Hence you should avoid the anti-aliased decimation but instead employ a bandpass decimation strategy. A direct subsampling by 32 would provide a convincing result.

Same goes with interpolation part. After you expand your signal, you should select the proper band of image spectrum by a suitable bandpass filter, rather than a lowpass filter from a conventional interpolator. A direct expansion may convince as its output will already be bandpass filtered later by the polyphase filterbank.

That being said, your main problem is the document that you based your code. That's far too out of signal processing context. I suggest you read the MPEG audio standard book of BOSI. But it's up to you of course.

The following is a working corrected MATLAB / OCTAVEcode:

% S0 - Load the prototype lowpass filter impulse response:
% --------------------------------------------------------
load h2.mat;         % h[n] is the prototype lowpass filter... of length 512
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;


% S1 - Create the 32 band filter bank by cosine modulation
% --------------------------------------------------------
numbands = 32;                  % number of bands
n=0:L-1;

hha=zeros(numbands,L);         % band of filters
for k=0:1:numbands-1 
   hha(k+1,:) = h.*cos( ( (2*k+1)*pi*(n-16) ) / (2*numbands) );
end

figure
plot(linspace(-1,1,4*L),20*log10(abs(fftshift(fft(hha(1,:),4*L)))));
title('ANALYSIS FILTER BANK');
hold on
for k=1:numbands
    plot(linspace(-1,1,4*L),20*log10(abs(fftshift(fft(hha(k,:),4*L)))));
end

% S2 - Design the synthesis filterbank:
% -------------------------------------
numbands = 32;                     % number of banks
n=0:L-1;

hhs = zeros(numbands,L);           % band of filters
for k=0:1:numbands-1 
   hhs(k+1,:) = h.*cos( ( (2*k+1)*pi*(n+16) ) / (2*numbands) );
end


% S3 - Generate a test input signal:
% ----------------------------------
N = 2*1024;      % NOTE: According to the MPEG specification a determined LENGTH must be applied.
wav_in = sin(0.05*pi*[0:N-1]);        % pure sine tone

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


% S4 - Apply a wav file to it,  get every output of the every sub-band:
% --------------------------------------------------------------------
yyd = zeros( numbands, floor(N/numbands));   % decimated outputs..
yyu = zeros( numbands , N);

for k=1:1:numbands
    temp = conv( wav_in,hha(k,:));
    yyu(k,:) = temp(L/2+1:L/2+N);
    yyd(k,:) = yyu(k,1:numbands:end);
end


% S6 - Apply synthesis on the decimated signal
% --------------------------------------------
ys = zeros(1, N);
yi = zeros(1,N);

for k=1:1:numbands
    yi(1:numbands:end) = yyd(k,:);
    temp = conv(yi,hhs(k,:));
    ys = ys + temp(L/2+1:L/2+N);
end
ys = numbands*ys;

% Display the result:    
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)))));

With the outputs being:

1- The prototype filter:

enter image description here

2- 32 channel filter bank:

enter image description here

3- A pure tone input:

enter image description here

4- An reconstructed back output:

enter image description here

$\endgroup$
  • 1
    $\begingroup$ You are right! Thanks! I was killing the signal with the LPF in the decimate part. Now it works. Thanks a lot for the rest of the comments. $\endgroup$ – f.gallardo Oct 30 '17 at 5:46

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.