文章目录
傅里叶仿真分析/正交镜像滤波器组分析及实现
一、傅里叶仿真分析
1、分析
由实验结果可知:
- 时域移位会造成相位谱变化,幅度谱不变
- 时域越窄,频域越宽;时域越宽,频域越窄
2、实现代码
clear;
clc;
close all;
x1=zeros(1,20);
x2=zeros(1,20);
x3=ones(1,20);
x1(1,2:6)=[1 2 3 2 1];
x2(1,end-5:end-1)=[1 2 3 2 1];
x3(1,10)=0;
x1k=fft(x1);
x2k=fft(x2);
x3k=fft(x3);
n=0:length(x1)-1;
figure()
% 时域波形
subplot(331);
stem(n,x1);
title("x1序列时域波形");
subplot(332);
stem(n,x2);
title("x2序列时域波形");
subplot(333);
stem(n,x3);
axis([0 20 0 1.1]);
title("x3序列时域波形");
% 幅度谱
subplot(334);
stem(n,abs(x1k));
title("x1序列频域幅度谱");
subplot(335);
stem(n,abs(x2k));
title("x2序列频域幅度谱");
subplot(336);
stem(n,abs(x3k));
axis([0 20 0 1.1]);
title("x3序列频域幅度谱");
% 相位谱
subplot(337);
stem(n,angle(x1k));
title("x1序列频域相位谱");
subplot(338);
stem(n,angle(x2k));
title("x2序列频域相位谱");
subplot(339);
stem(n,angle(x3k));
title("x3序列频域相位谱");
二、正交镜像滤波器组分析及实现
1、重建信号时存在的问题
(1)混叠失真。由抽取和内插产生的混叠和镜像带来的误差,导致分析滤波器组和综合滤波器组的频带不能完全分开;
(2)幅度失真。由于分析和综合滤波器组的频带在通带内不是全通函数,其幅频特性波纹产生的误差;
(3)相位失真。由滤波器相频特性的非线性所产生的误差;
(4)量化失真。由编、解码产生的误差,与量化噪声相似,这类误差无法完全消除,只能设法减小。
2、两通道正交镜像滤波器组理论
一个两通道正交镜像滤波器组如图1所示,在分析滤波器组一侧,输入信号(设为宽带信号)被分成K个子频带信号(窄带信号),通过抽取可降低采样率;在综合滤波器一侧,通过零值内插和带通滤波可以重建原来的信号。
消除混叠失真的方法:
G
0
(
z
)
=
H
0
(
z
)
G_{0}\left( z\right) =H_{0}\left( z\right)
G0(z)=H0(z)
G 1 ( z ) = − H 1 ( z ) G_{1}\left( z\right) =-H_{1}\left( z\right) G1(z)=−H1(z)
当两通道无混叠滤波器组的分解滤波器满足:
H
1
(
z
)
=
H
0
(
−
z
)
H_{1}\left( z\right) =H_{0}\left( -z\right)
H1(z)=H0(−z)
时,该滤波器组为正交镜像滤波器组(quadrature mirror filter bank,QMFB)。且滤波器的幅度特性满足
∣
H
1
(
e
j
w
)
∣
=
∣
H
0
(
e
j
(
π
+
w
)
)
∣
\left\vert H_{1}\left( e^{jw}\right) \right\vert =\left\vert H_{0}\left( e^{j\left( \pi +w\right) }\right) \right\vert
∣∣H1(ejw)∣∣=∣∣∣H0(ej(π+w))∣∣∣
由上式可知
G
0
(
z
)
=
H
0
(
z
)
,
G
1
(
z
)
=
−
H
1
(
z
)
=
−
H
0
(
−
z
)
G_{0}\left( z\right) =H_{0}\left( z\right) ,G_{1}\left( z\right) =-H_{1}\left( z\right) =-H_{0}\left( -z\right)
G0(z)=H0(z),G1(z)=−H1(z)=−H0(−z)
可以看出,两通道分析滤波器组和综合滤波器组都由H0(z)决定,如果 是一个低通滤波器,H1(z)、G0(z) 、G1(z) 则分别是高通、低通、高通滤波器。
目标:在QMFB设计过程中,希望H0(z)和H1(z)能把输入信号分成两个子带信号,且频谱尽量不重叠。
这就使得对H0(z)和H1(z)的通带和阻带性能要求较高。
3、已有方法
- 特征值法
- 最小二乘法
- 遗传法、
- 多项式分解法等
问题:计算相对复杂、参数不易确定、实现相对较难。
4、完全重建QMFB的仿真实现
步骤一:确定N值。
选定w=0.45(固定 ),以不同的输入信号,改变N的大小求出相应的mse,通过比较得到最优的N值。
步骤二:确定w值。
选定N=41,以不同的输入信号,改变w的大小求出相应的mse,通过比较得到最优的w值。
根据上述结果,选取各滤波器的阶数N=41,滤波器h0的通带截止频率w=0.43
仿真结果如下:
(1)滤波器H0(z)和H1(z)的幅度响应
由图可知,H0(z)和H1(z)两个滤波器以w=0.5为中心对称,把输入信号分成两个子带信号。
(2)幅度响应关系误差
由图可知,利用该方法实现可使误差减小至可忽略水平。
(3)输入信号和输出信号比较
输入信号 | 理想重建信号 |
---|---|
- 重建信号和输入信号相比有一定时间的延迟
- 重建信号的幅度与输入信号的幅度相差一个固定倍数
(4)理想输出信号与重建输出信号的偏差
二者偏差都在10-3数量级,故可以将重建输出信号近似代替理想输出信号,实现近似的完全重建。
5、实现代码
clear;
clc;
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]);
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);