【完全重建QMF滤波器组的设计】 滤波器组
完全重建QMF 滤波器组的设计
摘 要
本文从信号处理的角度分析了两通道完全重建正交镜像滤波器组(PR QMFB) 理论，由于重建过程中各种失真带来的误差影响，完全重建QMFB 具有一定的困难。为此，本文提出了一种基于matlab 的完全重建QMFB 的方法，并通过仿真获得了完全重建QMFB 的最佳系数，从而使重建的QMFB 性能进一步提高，接近完全重建的目标，并给出了优化前后的模拟仿真结果。从仿真结果中可以看出，这种方法得到的重建信号误差相对较小，具有较高的精度，且实现起来具有直观、简单、快捷的优点。
关键词：滤波器设计 完全重建 QMFB matlab
1、引言
多抽样率技术的应用如今越来越广泛，它可以实现一个信道中的多路通信：可以用于数据压缩（如分频带编码后进行存储和传输）；可以用于图像处理和语音处理；可以用于窄带数字滤波的实现；可以用于保密系统等，多抽样率滤波器组已成为信号处理领域强有力的工具。完全重建QMF 滤波器组最大的优点是在对信号进行抽取后，可以根据每个子带的不同特征分别进行处理，而插值和合成环节又能消除信号失真的各种因素，因此被广泛用于语音处理、图像处理、国防通信和小波变换中。目前，完全重建QMF 滤波器组的设计有多种优化设计方法，如特征值法、最小二乘法[1]、遗传法、多项式分解法[2]等等，都将信号失真降到了很小的范围，但是这些方法计算复杂、参数不容易确定、程序编写较难。而利用matlab 的信号处理功能和运算功能，可以快速、有效的设计完全重建QMF 滤波器组，且具有较高的精度。
2、两通道正交镜像滤波器组理论
一个两通道正交镜像滤波器组如图1所示，在分析滤波器组一侧，输入信号（设为宽带信号）被分成K 个子频带信号（窄带信号），通过抽取可降低采样率；在综合滤波器一侧，通过零值内插和带通滤波可以重建原来的信号[3]。
图 1两通道正交镜像滤波器组
ˆ(nT 1) 与分析滤波器组原来输入的信号如果综合滤波器组最后输出的信号x x (nT 1) 有如下关系：
ˆT ) =c [（2.1）x (-n 0n ) 1 T ] x (n 1
ˆ(nT 1) 是对x (nT 1) 的完全重建（Perfect 其中，c 和n 0均为固定的常数。称x
Reconstruction ），简称PR 。也就是说，输出信号是输入信号的延迟样本，在幅度上它们只相差一个固定的倍数c ，在时间上有一个固定的延迟n 0T 1。对应的频域关系为：
n 0ˆ(z ) =c -z X
（2.2） X ( ) z
对于一个给定的信号，经过分析滤波器后，再进行抽取、编码、传输，可以通过零值内插、综合滤波器滤波、求和运算得到恢复和重建。但是重建后的信号并不能与原始信号完全相同，两者之间存在着误差，主要包括：
（1）混叠失真。由抽取和内插产生的混叠和镜像带来的误差，导致分析滤波器组和综合滤波器组的频带不能完全分开；
（2）幅度失真。由于分析和综合滤波器组的频带在通带内不是全通函数，其幅频特性波纹产生的误差；
（3）相位失真。由滤波器相频特性的非线性所产生的误差；
（4）量化失真。由编、解码产生的误差，与量化噪声相似，这类误差无法完全消除，只能设法减小[4]。
因此，在设计QMF 组时，就需要综合考虑如何减小和消除上述的各类误差。消除混叠失真一种简单形式采取：
G 0(z ) =H 0(z ) （2.3） G 1(z ) =-H 1(z ) （2.4）
当两通道无混叠滤波器组的分解滤波器满足：
H 1(z ) =H 0(-z ) （2.5）
时，该滤波器组为正交镜像滤波器组（quadrature mirror filter bank，QMFB ）。且滤波器的幅度特性满足
H 1(e jw ) =H 0(e j (π+w ) ) （2.6）
由（2.3）、（2.4）、（2.5）式可知
G 0(z ) =H 0(z ) G 1(z ) =-H 1(z ) =-H 0(-z ) （2.7）可以看出，两通道分析滤波器组和综合滤波器组都由H 0(z ) 决定，如果H 0(z ) 是一个低通滤波器，H 1(z ) 、 G 0(z ) 、G 1(z ) 则分别是高通、低通、高通滤波器。
3、完全重建QMFB 遇到的问题和解决办法
在完全重建QMFB 过程中，希望设计的滤波器通带尽量平、过渡带尽量窄，且阻带尽可能快速衰减。但是由假设得到的H 0(z ) 和H 1(z ) 不能满足这些要求，因此它们没有实用的意义。可以看出，在式H 1(z ) =H 0(-z ) 指定分析滤波器组的条件下想要做到既无混叠失真又无幅度失真是不可行的，也就是做不到QMFB 的完全重建。实际上，尽管PR 是目的，但是滤波器组的核心作用是子带分解。在QMFB 设计过程中，希望H 0(z ) 和H 1(z ) 能把输入信号分成两个子带信号，且频谱尽量不重叠，这就使得对H 0(z ) 和H 1(z ) 的通带和阻带性能要求较高。
解决办法：
（1）用FIR QMF 滤波器组，去除相位失真的前提下，尽可能的减小幅度失真，近似实现完全重建；
（2）用IIR QMF滤波器组，去除幅度失真，不考虑相位失真，近似实现完全重建；
（3）修正QMF 滤波器H 1(z ) =H 0(-z ) 的关系，去考虑更合理的形式，从而实现完全重建[5]。
综上分析可知，实现完全重建QMFB 并不简单，一般只能做到近似重建，近似程度取决于设计的优化。
4、完全重建QMFB 的设计
目前，完全重建QMFB 的设计有多种优化方法，如特征值法、最小二乘法、遗传法、多项式分解法等等，这些方法都将信号失真降到了很小的范围，但是计算相对复杂、参数不易确定、实现相对较难。
利用matlab 实现完全重建QMFB 的设计，只需要知道各滤波器的阶数N 和滤波器h 0的通带截止频率ω，就可以得到完全重建QMFB 的分析、综合滤波器组的时域形式h 0、h 1、g 0、g 1，误差较小且能达到良好的精度。其中，N 必须为奇数，ω必须小于0.5。
完全重建QMFB 需要找到合适的N 和ω的值，使得重建效果最好。为此，首先选定ω，对应不同的输入信号x (n ) ，改变N 的大小求出均方误差mse ，通过比较得到最优的N 值；其次，固定N ，对应不同的输入信号x (n ) ，改变ω的大小求出均方误差mse ，通过比较得到最优的ω值。至此，找到了使得完全重建QMFB 效果最好的参数N 和ω。
5、完全重建QMFB 的仿真实现
步骤一：确定N 值。
选定ω=0.45，以不同的输入信号x 1(n ), x 2(n ), x 3(n ) ，改变N 的大小求出相应的mse ，通过比较得到最优的N 值。
x 1(n ) =δ(n -1) +δ(n -2) +2δ(n -5) +2δ(n -6) +2δ(n -7) +1.5δ(n -16) +1.5δ(n -17)
+1.5δ(n -18) +δ(n -23) +δ(n -24) +3δ(n -32) +3δ(n -33) +3δ(n -34) x 2(n ) =δ(n ) +δ(n -1) +δ(n -2) +2δ(n -8) +2δ(n -9) +2δ(n -10) +3δ(n -15) +3δ(n -16)
+3δ(n -17) +4δ(n -23) +4δ(n -24) +4δ(n -25) +3δ(n -32) +3δ(n -33) +3δ(n -34) +2δ(n -40) +2δ(n -41) +2δ(n -42) +δ(n -48) +δ(n -49) +δ(n -50)
x 3(n ) =sin(nT ), T =0.2,1≤n ≤500 仿真结果如表1：
表1 均方误差图 (ω=0.45，N 变化)
从表1中可以看出，对输入信号一般有：N>127时，不能用matlab 实现QMFB 的设计；N=41时mse 最小，重建效果最好。 步骤二：确定ω值。
选定N=41，以不同的输入信号x 1(n ), x 2(n ), x 3(n ) ，改变ω的大小求出相应的mse ，通过比较得到最优的ω
值。 仿真结果如表2：
表2 均方误差图 (N=41，ω变化)
从表2 中可以看出，当ω=0.43时，重建效果最好，ω>0.48时，不能实现。这样，便获得了最佳的N 和ω。
根据上述结果，选取各滤波器的阶数N=41，滤波器h 0的通带截止频率ω=0.43，输入信号x (n ) =sin(nT ), T =0.2,1≤n ≤500， 仿真结果如下：
（1）滤波器H 0(z ) 和H 1(z ) 的幅度响应
100-10-20-30
幅度, d B
-40-50-60-70-80-90-100
ω/π
图 2
（2）幅度响应关系误差（理论公式如下）
H 1(e jw ) +H 0(e j (π+w ) ) ≈1,
2
2
10log H 1(e jw ) +H 0(e j (π+w ) ) ≈0
22
0.040.03
0.020.01
误差, d B
0-0.01-0.02-0.03-0.04
00.10.20.30.4
0.5ω/π
0.60.70.80.91
图 3
（3）输入信号
50
100
150
200
250
300
350
400
450
500
图 4
（4）理想输出信号（红色线）与重建的输出信号（蓝色线）
100
200
300
400
500
600
图 5
（5）理想输出信号与重建输出信号的偏差
-5
[***********]
[1**********]0
图 6
为了检验完全重建的程度，以理想输出信号与重建输出信号的偏差及其均方误差mse 的大小来判定，e 和mse 的值越小，完全重建的效果就越好。在各滤波器的阶数
N=41，滤波器h 0的通带截止频率ω=0.43，输入信号
x (n ) =sin(nT ), T =0.2,1≤n ≤500的条件下，均方误差mse=2.8074e-010，且从图中可
以看出其理想输出信号与重建输出信号的偏差相对较小，因此，可以实现QMFB 的近似完全重建。
在QMFB 的设计仿真中，发现1
Magnitude Response (dB)
M a g n i t u d e (d B )
Normalized Frequency (⨯π rad/sample)
图 7
h 0的幅度响应
综上可知，可以用matlab 实现完全重建QMFB 的设计，具有直观、简便、快捷的优点，极大的减轻了工作量，得到的均方误差也较小，具有相对较高的精度，最优参数一般可以选取N=41，至此，就完成了两通道正交镜像滤波器组(QMFB)
ω=0.43。
的设计。
6、结束语
从两通道正交镜像滤波器组(QMFB)理论出发，可以得到理想意义下的QMFB 设计，但是这种设计的复杂性高，具有一定的困难，且没有实际使用的意义。而利用matlab 可以快速有效的设计正交镜像滤波器组，且设计简单快捷。通过变化参数N 和 ，得到了实际意义下的取值范围及对应的幅度特性和最优参数值，可以大大节约时间、减少工作量、有效利用资源。
7、总结
通过多抽样率信号处理这门课程的学习，我学到了很多以前没有学过的知识，多抽样率系统的理论，包括整数倍抽取和内插、分数倍抽样率转换；多抽样率系统中的网络结构和实现，包括多抽样率系统的恒等变换、多项结构和多级实现；滤波器组理论以及多抽样率技术的应用。其中，我对正交镜像滤波器组的完全重建，即(PR QMFB) 进行了深入的研究，从两通道正交镜像滤波器组理论的研究，到matlab 设计寻找实际意义下的取值范围和最优参数，最后仿真得到具有较高精度的对比结果。经过了理论学习——提出问题——解决问题——发现问题——解决办法——matlab 设计——仿真实现，在整个过程中不仅学到了有关这方面的理论知识，更重要的是培养了独立思考和解决问题的能力，对今后的学习有很大的帮助，从中收获颇多。
参考文献
[1] 赖晓平. FIR数字滤波器的递推最小二乘设计算法[J]，Sep.1999,260-264.
[2] 石光明，焦李成. 多项式分解理论与低时延完全重构两通道滤波器组的设计[J]. 电子与信息学报，2002,24（7），910-915.
[3] 陶然，张惠云，王越. 多抽样率数字信号处理理论及应用[M]. 北京：清华大学出版社2007
[4] 胡广书. 现代信号处理教程[M]. 北京: 清华大学出版社，2004.142-144,161-169
[5] 张晨，李双田. 设计QMF 组的一种新算法及基于GA 的优化，信号处理[J]. 2005,4,120-125
源代码汇总：
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]);
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);
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);
%%%%%%%%%%理想输出信号与重建输出信号%%%%%%% figure(5);
axis([0,500,-1,1]);
plot(xshift,"r");hold on;
plot(xidle,"-");
axis([0,600,-1.1,1.1]);
%%%%%%%理想输出信号与重建输出信号的偏差%%%%%% 理想输出信号与重建的输出信号的偏差
figure(6);
plot(xshift-xidle);