我已经测试了你的代码。
抽取和上采样阶段存在问题。当您以未抽取的速率运行分析和合成滤波器时，您将获得完美的重建（加上滤波器移位）
首先对于对群延迟敏感的应用，最好使用conv / filter功能，然后补偿群延迟。
此外，抽取函数应用低通抗混叠滤波器，考虑到它的输出是带通信号而不是低通信因此，您应该避免使用抗锯齿抽取，而是采用带通抽取策略。32 的直接二次抽样将提供令人信服的结果。k
插值部分也是如此。扩展信号后，您应该通过合适的带通滤波器而不是来自传统插值器的低通滤波器来选择适当的图像频谱带。直接扩展可能会令人信服，因为它的输出稍后已经被多相滤波器组带通滤波。
话虽这么说，您的主要问题是您基于代码的文档。这远远超出了信号处理的范围。我建议你阅读 BOSI 的 MPEG 音频标准书。但这当然取决于你。
以下是经过修正的 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)))));
输出为：
1-原型过滤器：
2- 32 通道滤波器组：
3-纯音输入：
4-重建的反向输出：
