完全重建QMF滤波器组的设计

作业一：利用FFT分析信号

分析两个信号：一个信号在前端，另一个信号在后端（相当于第一个信号时移）
x1：正态分布的中心点在10，x2：正态分布的中心点在90
matlab实现：

clear;
clc;
close all;

N=100;
n=[0:1:N-1];
x1=normpdf(n,10,1);
x2=normpdf(n,90,1);

%绘制原始时域信号波形
figure(1)
subplot(2,1,1)
stem(x1)
title('信号1的时域波形')
subplot(2,1,2)
stem(x2)
title('信号2的时域波形')

Y1=fft(x1);
Y2=fft(x2);
figure(2)
subplot(3,1,1)
stem(abs(Y1));
title('Y1的幅度谱');
subplot(3,1,2)
stem(abs(Y2));
title('Y2的幅度谱');
subplot(3,1,3)
stem(abs(Y1)-abs(Y2));
title('Y1-Y2');

figure(3)
subplot(2,1,1)
stem(angle(Y1));
title('Y1的相位谱');
subplot(2,1,2)
stem(angle(Y2));
title('Y2的相位谱');

figure1：两原始时域信号

figure2：fft后的幅度谱

figure3：fft后的相位谱

对比两图fft后的幅度谱和相位谱，发现二者幅度谱相同，figure2中图三，二者的差在10^-16上下变化，说明相差很小，只是相位发生变化

作业二：完全重建QMF滤波器组的设计

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

二、误差分析
对于一个给定的信号，经过分析滤波器后，再进行抽取、编码、传输，可以通过零值内插、综合滤波器滤波、求和运算得到恢复和重建。但是重建后的信号并不能与原始信号完全相同，两者之间存在着误差，主要包括：
（1）混叠失真。由抽取和内插产生的混叠和镜像带来的误差，导致分析滤波器组和综合滤波器组的频带不能完全分开；
（2）幅度失真。由于分析和综合滤波器组的频带在通带内不是全通函数，其幅频特性波纹产生的误差；
（3）相位失真。由滤波器相频特性的非线性所产生的误差；
（4）量化失真。由编、解码产生的误差，与量化噪声相似，这类误差无法完全消除，只能设法减小。

三、设计完全重建QMFB
在完全重建QMFB过程中，希望设计的滤波器通带尽量平、过渡带尽量窄，且阻带尽可能快速衰减。实际上，无法做到完全重建QMFB，尽管PR是目的，但是滤波器组的核心作用是子带分解。在QMFB设计过程中，希望H0(z)和H1(z)能把输入信号分成两个子带信号，且频谱尽量不重叠，这就使得对H0(z)和H1(z)的通带和阻带性能要求较高。
解决办法：
（1）用FIR QMF滤波器组，去除相位失真的前提下，尽可能的减小幅度失真，近似实现完全重建；
（2）用IIR QMF滤波器组，去除幅度失真，不考虑相位失真，近似实现完全重建；
（3）修正QMF滤波器H1(z)=H0(-z)的关系，去考虑更合理的形式，从而实现完全重建。
综上分析可知，实现完全重建QMFB并不简单，一般只能做到近似重建，近似程度取决于设计的优化。

四、重建QMFB的实现
利用matlab实现完全重建QMFB的设计，只需要知道各滤波器的阶数N和滤波器h0的通带截止频率w，就可以得到完全重建QMFB的分析、综合滤波器组的时域形式h0、h1、g0、g1，误差较小且能达到良好的精度。其中，N必须为奇数，w必须小于0.5。
完全重建QMFB需要找到合适的N和w的值，使得重建效果最好。为此，首先选定w，对应不同的输入信号x(n)，改变N的大小求出均方误差mse，通过比较得到最优的N值；其次，固定N，对应不同的输入信号x(n)，改变的大小求出均方误差mse，通过比较得到最优的w值。至此，找到了使得完全重建QMFB效果最好的参数N和w。

五、matlab代码实现

close all;
N=41;
w=0.43;
[h0,h1,g0,g1]=firpr2chfb(N,w);
[H1z,w]=freqz(h0,1,512);
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]); 
title('滤波器h0与h1的幅度响应');
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'); 
title('幅度响应关系误差');
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);
hhp=mfilt.firdecim(2,h1);
glp=mfilt.firinterp(2,g0);
ghp=mfilt.firinterp(2,g1);
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);
title('输入信号input');
%%%%%%%%%%理想输出信号与重建输出信号%%%%%%%
figure(5);
axis([0,500,-1,1]); 
plot(xshift,'r');hold on;
plot(xidle,'b');
title('理想输出信号与重建输出信号');
%%%%%%理想输出信号与重建输出信号的偏差%%%%%%
figure(6);
plot(xshift-xidle);
title('理想输出信号与重建输出信号的偏差');

六、运行结果
figure1滤波器H0(z)和H1(z)的幅度响应

其中，实线为H0，低通滤波器。虚线为H1，高通滤波器

figure2幅度响应关系误差（理论公式如下）

figure3 h0的幅度

figure4输入信号

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

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

通过figure5和figure6可看出，理想输出信号与实际重建输出信号的误差非常小。