Norm MPEG-1 Layer III (Mp3) Filter banks distortion

Question

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

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

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?).

Answer 1

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:

2- 32 channel filter bank:

3- A pure tone input:

4- An reconstructed back output: