完全重建QMF滤波器

任务一

输入两个波形，一个脉冲在前，一个脉冲在后，观察其傅里叶变换后区别。

a = zeros(1, 100);
a(1:10) = 1; %冲击在前的信号
subplot(231);
stem(a);
title('原信号');
subplot(232);
f = fft(a);
stem(abs(f));%冲击在前信号的幅度
title('幅度谱');
subplot(233);
stem(angle(f)) %冲击在前信号的相位
title('相位谱');

a = zeros(1, 100);
a(90:100) = 1; %冲击在后的信号
subplot(234);
stem(a);
title('原信号'); %冲击在后信号的幅度
subplot(235);
f = fft(a);
stem(abs(f));
title('幅度谱');
subplot(236);
stem(angle(f)); %冲击在后信号的相位
title('相位谱');

结果：

可以看到除了相位，二者的幅度并没有区别。其实不难理解，毕竟在时域上的移动对应到频域就是相位的变化。

任务二

一个两通道正交镜像滤波器组如下图所示，在分析滤波器组一侧，输入信号（设为宽带信号）被分成K个子频带信号（窄带信号），通过抽取可降低采样率；在综合滤波器一侧，通过零值内插和带通滤波可以重建原来的信号。

将一个信号拆成两个信号，可以减小其动态范围，此后再进行量化编码时可以减小量化误差。同时可以间隔取值，从而降低每一通道的抽样点数，使得数据量保持不变。

为了消除混叠失真，令
$$G_0(z)=H_0(z)$$
$$G_1(z)=-H_1(z)$$
那么就有
$$H_0(z)=G_0(z)=-G_1(z)$$
其中，若$H_0(z)$为低通，则$H_1(z),G_0(z),G_1(z)$分别为高通、低通、高通滤波器。

代码如下：

N=41;
w=0.43;
[h0,h1,g0,g1]=firpr2chfb(N,w); %wo-channel FIR filter bank for perfect reconstruction
[H1z,~]=freqz(h0,1,512); %返回n点的频率相应向量
H1_abs=abs(H1z);H1_db=20*log10(H1_abs);
[H2z,w]=freqz(h1,1,512);
H2_abs=abs(H2z);H2_db=20*log10(H2_abs);
%%%%%%%%%%滤波器h0和h1的幅度响应%%%%%%%%%%
figure(1); 
plot(w/pi,H1_db,'-',w/pi,H2_db,'--'); 
axis([0,1,-100,10]); 
grid 
xlabel('\omega/\pi');ylabel('幅度,dB'); 
sum1=H1_abs.*H1_abs+H2_abs.*H2_abs; 
d=10*log10(sum1);
%%%%%%%%%%%%幅度响应关系误差%%%%%%%%%%%%%
figure(2) 
plot(w/pi,d);grid; 
xlabel('\omega/\pi');ylabel('误差,dB'); 
axis([0,1,-0.04,0.04]); 
%%%%%%%%%%%%%x1(n)%%%%%%%%%%%%%%%%%%%%%
x=zeros(1,500);
x(2)=1;x(3)=1;
x(6)=2;x(7)=2;x(8)=2;
x(17)=1.5;x(18)=1.5;x(19)=1.5;
x(24)=1;x(25)=1;
x(33)=3;x(34)=3;x(35)=3;
%%%%%%%%%%%%%%x2(n)%%%%%%%%%%%%%%%%%%%%
x=zeros(1,500);
x(1)=1;x(2)=1;x(3)=1;
x(9)=2;x(10)=2;x(11)=2;
x(16)=3;x(17)=3;x(18)=3;
x(24)=4;x(25)=4;x(26)=4;
x(33)=3;x(34)=3;x(35)=3;
x(41)=2;x(42)=2;x(43)=2;
x(49)=1;x(50)=1;x(51)=1;
%%%%%%%%%%%%%%x3(n)%%%%%%%%%%%%%%%%%%%%
n=1:500;
T=0.2;
x=sin(n*T);
hlp=mfilt.firdecim(2,h0);  % returns a direct-form FIR polyphase decimator object
hhp=mfilt.firdecim(2,h1);  % returns a direct-form FIR polyphase decimator object
glp=mfilt.firinterp(2,g0); % returns a FIR polyphase interpolator object Hm with aninterpolation factor of 2 and gain equal to 2.
ghp=mfilt.firinterp(2,g1); % returns a FIR polyphase interpolator object Hm with an interpolation factor of 2 and gain equal to 2.
x0=filter(hlp,x);
x0=filter(glp,x0);
x1=filter(hhp,x);
x1=filter(ghp,x1);
xidle=x0+x1;
xshift=[zeros(1,N) x(1:end-N)];
e=xidle-xshift;
mes=sum(abs(e).^2)/length(e)
fvtool(h0)
%%%%%%%%%%%%输入信号%%%%%%%%%%%%%%%%%%
figure(4);
plot(x);
%%%%%%%%%%理想输出信号与重建输出信号%%%%%%%
figure(5);
plot(xshift(1:40),'r');hold on;
plot(xidle(1:40),'-');
%%%%%%%理想输出信号与重建输出信号的偏差%%%%%%
figure(6);
plot(xshift-xidle);

结果：
滤波器$h_0$和$h_1$的幅度响应

使用一对镜面对称的滤波器实现子代分解。

幅度响应关系误差

可以看到误差非常小。

$h_0$的幅度响应

输入信号

理想输出信号（红色线）与重建的输出信号（蓝色线）

这里由于差别非常小，基本上看不出差别。故使用放大的局部：

plot(xshift(1:40),'r');hold on;
plot(xidle(1:40),'-');

结果

可以发现误差大概在$10^{-5}$量级。在本例中，相对于原信号，输出信号的幅度并未发生改变，只是发生了时延。

理想输出信号与重建输出信号的偏差

误差大概在$10^{-5}$量级，几乎可以算作是无失真。