MATLAB在数字信号处理中的应用

1、第4章 MATLAB在数字信号处理中的应用 MATLAB 中的信号处理工具箱提供了很多数字信号处理中需要用到的函数及解决方法。在数字信号处理中所学习到的问题，诸如滤波器的设计、自适应滤波，维纳滤波、卡尔曼滤波等理论都可以通过MATLAB仿真得到实现和验证。本章提供了12个MATLAB在数字信号处理课程学习中的应用实例，可作为学习数字信号处理课程的参考。第4章 MATLAB在数字信号处理中的应用知知 识识 架架 构构 1实验目的 (1)掌握脉冲响应不变法设计IIR数字滤波器的具体设计方法。 (2)熟悉脉冲响应不变法设计低通滤波器的仿真。 2实验原理 脉冲响应不变法是从滤波器的脉冲响应出发，使数字

2、滤波器的单位脉冲响应序列h(n)模仿模拟滤波器的冲击响应ha(t)，使h(n)正好等于ha(t)的采样值，即：4.1 IIR 4.1 IIR 带通滤波器设计带通滤波器设计第4章 MATLAB在数字信号处理中的应用 T为采样周期。若以Ha(s)及H(z)分别表示ha(t)的拉氏变换及h(n)的Z变换，即：4.1 IIR 4.1 IIR 带通滤波器设计带通滤波器设计第4章 MATLAB在数字信号处理中的应用 根据采样序列Z变换与模拟信号拉氏变换的关系得: 上式表明，采用脉冲响应不变法将模拟滤波器变换为数字滤波器时，它所完成的S平面到Z平面的变换，正是以前讨论的拉氏变换到Z变换的标准变换关系。 脉冲

3、响应不变法特别适用于用部分分式表达的传递函数，模拟滤波器的传递函数若只有单阶极点，且分母的阶数高于分子阶数，则可表达为部分分式形式：4.1 IIR 4.1 IIR 带通滤波器设计带通滤波器设计第4章 MATLAB在数字信号处理中的应用 其拉氏反变换为： 对ha(t)采样就得到数字滤波器的单位脉冲响应序列： 再对h(n)取Z变换，得到数字滤波器的传递函数： 第二个求和为等比级数之和：4.1 IIR 4.1 IIR 带通滤波器设计带通滤波器设计第4章 MATLAB在数字信号处理中的应用 要收敛的话，必有： 比较部分分式形式的Ha(s)和上式H(z)可以看到，把S平面上的极点变换到Z平面上对应的极点

4、，而Ha(s)与H(z)中部分分式所对应的系数不变。如果模拟滤波器是稳定的，则所有极点都在S左半平面那么变换后H(z)的极点也都在单位圆内，因此数字滤波器保持稳定。 所以有： 3仿真思路 在MATLAB中，可以用下列函数辅助设计IIR数字滤波器。 (1)利用cheb1ord可以确定低通原型和切比雪夫滤波器的阶数和截止频率。 (2)num，den=cheby1(N，Wn)，num，den=cheby2(N，Wn) (切比雪夫1 型和2 型)可以进行滤波器的设计。 (3)利用impinvar可以完成脉冲响应不变法的模拟滤波器到数字滤波器的转换。 4程序代码4.1 IIR 4.1 IIR 带通滤波器

5、设计带通滤波器设计第4章 MATLAB在数字信号处理中的应用 clear wp=6*pi*103;ws=9*pi*103;ap=1,as=15; Fs=30*103; wp1=wp/Fs;ws1=ws/Fs; N,WC=cheb1ord(wp,ws,ap,as,s); b,a=cheby1(N,ap,WC,s) ； bz,az=impinvar(b,a,Fs)； w0=wp1,ws1； Hx=freqz(bz,az,w0)； H,W=freqz(bz,az)； dbHx=-20*log10(abs(Hx)/max(abs(H)； plot(W,abs(H); xlabel(相对频率);ylab

6、el(幅频); grid4.1 IIR 4.1 IIR 带通滤波器设计带通滤波器设计第4章 MATLAB在数字信号处理中的应用 5运行结果与分析 结论：由图4.2可知，切比雪夫型滤波器的振幅特性在通带内是等波纹的，在阻带内是单调的。4.1 IIR 4.1 IIR 带通滤波器设计带通滤波器设计第4章 MATLAB在数字信号处理中的应用6思考题(1) 带通滤波器的幅度有什么特点？图4.2中滤波器的下降沿怎么样才能更陡峭？(2) 用本节所学的知识设计切比雪夫型滤波器。图 4.2： 1实验目的 (1)熟悉用双线性变换法设计IIR数字滤波器的原理与方法。 (2)熟悉带阻数字滤波器设计方法。 (3)通过观

7、察对实际心电图信号的滤波作用，获得数字滤波的感性知识。 2实验原理 1)用双线性变换法设计IIR数字滤波器数字低通技术指标为：4.2 IIR 4.2 IIR 带阻滤波器设计带阻滤波器设计第4章 MATLAB在数字信号处理中的应用 模拟低通的技术指标为：查巴特沃斯归一化低通滤波参数表可得：第4章 MATLAB在数字信号处理中的应用设计巴特沃斯低通滤波器。阶数N计算如下。数字滤波器系统函数H(z)为：4.2 IIR 4.2 IIR 带阻滤波器设计带阻滤波器设计 A=0.09036 B1=1.2686, C1=0.7051 B2=1.0106, C2=0.3583 B3=0.9044, C3=0.2

8、155 可见H(z)是由3个二阶滤波器H1(z)、H2(z)、H3(z)级联组成的，如图4.3所示。4.2 IIR 4.2 IIR 带阻滤波器设计带阻滤波器设计第4章 MATLAB在数字信号处理中的应用 图 4.3 滤波器H(z)的组成图2) 带阻IIR数字滤波器实现过程 (1)按一定规格将数字滤波器的技术指标转为模拟低通滤波器的技术指标。 (2)根据转换后的技术指标使用滤波器阶数函数，确定最小阶数N 和截止频率Wc。 (3)利用最小阶数N产生模拟低通滤波器原型。 (4)利用截止频率Wc把模拟低通原型转化为模拟带阻滤波器。 (5)利用冲激响应不变法或双线性不变法把模拟滤波器转换成数字滤波器。

9、3仿真思路 (1)用双线性变换法设计一个低通IIR数字滤波器。4.2 IIR 4.2 IIR 带阻滤波器设计带阻滤波器设计第4章 MATLAB在数字信号处理中的应用 设计参数为：在通带内频率低于0.2rad时，最大衰减小于1dB，在阻带内0.3，频率区间上，最小衰减大于15dB。 (2)所设计的滤波器对实际心电图信号采样序列进行仿真滤波处理，观察滤波前后的图形。 (3)设计一个带阻IIR数字滤波器，其具体要求是：通带的截止频率wp1=650Hz，wp2=850Hz；阻带的截止频率ws1=700Hz，ws2=800Hz，带通内的最大衰减rp=0.1dB，阻带内的最小衰减为rs=50dB，采样频率

10、为Fs=2000Hz。4.2 IIR 4.2 IIR 带阻滤波器设计带阻滤波器设计第4章 MATLAB在数字信号处理中的应用4程序代码1) 心电图信号滤波代码 x=-4,-2,0,-4,-6,-4,-2,-4,-6,-6,. -4,-4,-6,-6,-2, 6,12, 8, 0,-16,. -38,-60,-84,-90,-66,-32,-4,-2,-4, 8,. 12,12,10,6,6,6,4,0,0,0,. 0,0,-2,-4,0,0,0,-2,-2,0,. 0,-2,-2,-2,-2,0; %存有高频干扰心电图信号序列 k=1; close all; Figure(1) subplot

11、(2,2,1) axis(0 56 -100 50);4.2 IIR 4.2 IIR 带阻滤波器设计带阻滤波器设计第4章 MATLAB在数字信号处理中的应用n=0:55;stem(n,x,.);hold on;n=0:60;m=zeros(61);plot(n,m);xlabel(n);ylabel(x(n);title(心电图信号采样序列x(n);B=0.09036 2*0.09036 0.09036;A=1.2686 -0.7051;A1=1.0106 -0.3583;A2=0.9044 -0.2155;while(k0.7,break,end %查找行号 end for l=1:4; w

12、pz2=0.4 0.3625 0.325 0.3; %低频的4 个高通截止数字频率 wsz2=0.25; rp2=3;rs2=60;第4章 MATLAB在数字信号处理中的应用 N2,wc2=buttord(wpz2(l),wsz2,rp2,rs2); %高通数字滤波器 Bz2,Az2=butter(N2,wc2,high); y2=filter(Bz2,Az2,x); %滤波后的单频信号 z2=max(abs(y2); %滤波后信号在离散域所取的最大值 subplot(3,3,m+6);plot(n,y2);grid;title( 高通后的波形);xlabel(n); ylabel(y) if

13、 z20.7,break,end; %查找列号 end TNr=TNr+tm(k,(5-l)*10(3-m); %将3位电话号码表示成一个3位数 end disp(接受断检测到的号码为：) disp(TNr) 5运行结果与分析 从图4.17中可知，3个数字的低频部分频谱相同， 高频部分不同，跟假设匹配。通过对滤出频率分析，查找表，可以得出输出的数字。第4章 MATLAB在数字信号处理中的应用 图4.17数字1、2、3的时域和滤波后的波形 6思考题 双音多频滤波器常用来解决什么问题？第4章 MATLAB在数字信号处理中的应用 1实验目的 (1)知道直接型滤波器转换为级联型滤波器的原理。 (2)了

14、解直接型滤波器转换为级联型滤波器后，级联型网络结构的优点。 (3)掌握直接型滤波器转换为级联型滤波器的实现方式及级联型滤波器函数输出的编程原理。 2实验原理 为了对连续的或离散的信号进行滤波处理，就必须构造出合适的实际结构。对于同样的系统函数H(z)或H(s)往往有不同的实现方案。常用的有直接形式和级联形式以及并联形式。下面主要讨论直接形式4.8 4.8 直接型与级联型滤波器的比较直接型与级联型滤波器的比较第4章 MATLAB在数字信号处理中的应用 到级联形式的转换以及二者的比较。 直接形式用延迟元件和乘法器以及加法器以给定的形式直接实现差分方程，设M=N=4，则差分方程为：4.8 4.8 直

15、接型与级联型滤波器的比较直接型与级联型滤波器的比较第4章 MATLAB在数字信号处理中的应用 在这种形式中将系统函数H(z)写成实系数二阶子系统的乘积形式。首先把分子和分母多项式的根解出，然后把每一对共轭复根或任意两个实根组合在一起，得到二阶子系统。在以下原理中假设N 为偶 数，于是可把上式转化为下面的形式：4.8 4.8 直接型与级联型滤波器的比较直接型与级联型滤波器的比较第4章 MATLAB在数字信号处理中的应用上式中的参量满足：二阶子系统为： 在工程实际中，一般把上式所示的结构称为双二阶环节(Biquad),它的输入是第(k 1) 个双二阶环节的输出，同时第k个双二阶环节的输出为第(k

16、+1)个双二阶环节的输入。 3仿真思路 本程序由1个主函数和3个子函数组成，主函数完成由下面的差分方程描述的直接型滤波器到级联型滤波器的转换。 15y(n)+10y(n1)+ 3y(n2)5y(n3)2y(n4) =x(n)4x(n1)+10 x(n2)25x(n3) +16x(n4) 子函数dir2cas实现将直接形式转化为级联形式，子函数stepseq(n0,n1,n2)，实现阶跃函数u(n-n0)，n1nn2。子函数casfilter是编制实现级联形式滤波器的函数。 级联结构中每一个二阶网络决定一对零点和一对极点，可以通过灵活地调整对应零、极点的系数来4.8 4.8 直接型与级联型滤波器

17、的比较直接型与级联型滤波器的比较第4章 MATLAB在数字信号处理中的应用 改变一对零、极点的位置，相对于直接型结构，其优点是调整方便，此外，级联结构中后面的网络输出不会流到前面，运算误差的积累相对于直接型也小。 4程序代码 1)主程序 b=1 -4 10 -25 16; a=15 10 3 -5 -2; b0,B,A=dir2cas(b,a) stepin=stepseq(0,0,5); format long; hcas=casfilter(b0,B,A,stepin)4.8 4.8 直接型与级联型滤波器的比较直接型与级联型滤波器的比较第4章 MATLAB在数字信号处理中的应用hdir=f

18、ilter(b,a,stepin)format;figure(1)subplot(2,1,1)stem(hcas); hold onplot(0 6,0 0);axis(0 6 -1.8 1);title(级联型的阶跃响应); hold off;subplot(2,1,2)stem(hdir);hold on;plot(0 6,0,0);axis(0 6 -1.8 1);title(直接型的阶跃响应); hold off;第4章 MATLAB在数字信号处理中的应用2) stepseq 函数 function x,n=stepseq(n0,n1,n2) %Generate x(n)=u(n-n0

19、);n1=n n0;3) casfilter 函数 function y=casfilter(b0,B,A,x) K,L=size(B); N=length(x); w=zeros(K+1,N); w(1,:)=x; for i=1:K w(i+1,:)=filter(B(i,:),A(i,:),w(i,:); end y=b0*w(K+1,:);第4章 MATLAB在数字信号处理中的应用4) dir2cas 函数 function b0,B,A=dir2cas(b,a) b0=b(1); b=b0/b0; a0=a(1); a=a/a0; b0=b0/a0; M=length(b); N=l

20、ength(a); if NM b=b zeros(1,N-M); else if NM a=a zeros(1,N-M); else NM=0; end K=floor(N/2);第4章 MATLAB在数字信号处理中的应用B=zeros(K,3);A=zeros(K,3);if K*2=M b=b 0; a=a 0;endbroots=cplxpair(roots(b);aroots=cplxpair(roots(a);For i=1:2:2*K; Brow=broots(i:1:i+1,:); Brow=real(poly(Brow); B(fix(i+1)/2),:)=Brow; Aro

21、w=aroots(i:1:i+1,:); Arow=real(poly(Arow); A(fix(i+1)/2),:)=Brow;End第4章 MATLAB在数字信号处理中的应用 5运行结果与分析 b0 = 0.0667 B = 1 0 0 1 0 0 A = 1 0 0 1 0 0 由图4.18可知，这两种形式结构的滤波器本质上是完全相同的，但相比于直接型结构，级联形式的滤波器可以更灵活地调整零、极点的系数来改变一对零、极点的位置，调整方便，而且级联结构中后面的网络输出不会流到前面，运算误差的积累也很小，因而得到了广泛应用。第4章 MATLAB在数字信号处理中的应用 图4.18 级联型滤波器

22、和直接型滤波器的阶跃响应 6思考题 (1)级联型滤波器和直接型滤波器各有什么优缺点？ (2)级联型滤波器是如何实现的？第4章 MATLAB在数字信号处理中的应用 1实验目的 (1)进一步了解自适应滤波原理。 (2)学习LMS自适应算法及其MATLAB仿真。 2实验原理 自适应滤波器由参数可调的数字滤波器和自适应算法两部分组成。输入信号x(n)通过参数可调数字滤波器后产生输出信号y(n)，将其与参数信号d(n)进行比较，形成误差信号e(n)。e(n)通过某种自适应算法对滤波器参数进行调整，最终使e(n)的均方值最小。 最小均方误差LMS准则的目的在于使滤波器输出与期望信号误差的平方的统计平均值最

23、小。LMS自适应横向滤波器如图4.19所示。4.9 4.9 自适应滤波自适应滤波LMSLMS算法算法第4章 MATLAB在数字信号处理中的应用 图4.19 最小均方误差准则的自适应横向滤波器原理图该自适应滤波器的输入矢量为：4.9 4.9 自适应滤波自适应滤波LMSLMS算法算法第4章 MATLAB在数字信号处理中的应用加权矢量为为：滤波器的输出为：y(n)相对于滤波器期望输出d(n)的误差为： 根据最小均方误差准则，最佳的滤波器参量应使得性能函数均方误差为最小。 假定输入信号x(n)和期望响应d(n)是联合平稳过程，那么在时刻n的均方误差是加权矢量的二次函数，其表示式为：4.9 4.9 自适

24、应滤波自适应滤波LMSLMS算法算法第4章 MATLAB在数字信号处理中的应用 均方误差是权向量的二次函数，它是一个上凹的抛物面，具有唯一的最小值，调解权向量使得均方误差最小，相当于沿抛物面下降寻找最小值。可以用梯度法来求该最小值，对权向量W 求导得到均方误差的梯度为： 在性能曲面上最佳权矢量对应点的梯度等于零，即：4.9 4.9 自适应滤波自适应滤波LMSLMS算法算法第4章 MATLAB在数字信号处理中的应用 该方程称为正则方程，由此解出最佳权向量称为维纳解。 利用上式求解，需要精确地知道输入信号和期望信号的先验统计知识，而且还要对矩阵求逆运算。最陡下降法可避免求逆运算，它通过递推的方式寻

25、求加权矢量的最优值，是LMS算法的理论基础。首先设置一个W的初值W(0)，可以想象，沿减小的方向 调整W，可以找到最佳权矢量。因为梯度方向是增加最快的方向，所以负梯度方向就是减少最快的方向。 最小均方算法是一种很有用很简单的估计梯度的方法，其突出特点是计算量小、易于实现，且不要求脱线计算。 LMS最核心的算法是使用平方误差代替均方误差，即4.9 4.9 自适应滤波自适应滤波LMSLMS算法算法第4章 MATLAB在数字信号处理中的应用所以将上式代入最陡下降法迭代计算权矢量的公式得：则LMS算法的基本关系式为： 3仿真思路 设计一个二阶加权自适应横向滤波器，对一个正弦信号进行滤波。实验通过设置不

26、同的收敛因子，由MATLAB程序图形观察滤波效果。讨论的重要性。 4程序代码 clear all ; fs =10000; t=0：1/fs ：1; sn =sin(2 * pi * t); %产生初始信号 n= randn (size (t); %产生高斯噪声 xn = sn + n; %信号加噪声 w = 0 0.5; %设置权初值 u=0.00026; %设置收敛因子 for i =1:length(t)-1第4章 MATLAB在数字信号处理中的应用 for i =1:length(t)-1 y (i+1)=n(i:i+1)*w; %噪声通过滤波器输出为y e (i+1)=xn (i+1

27、)-y(i+1); w = w+2*u*e(i+1)*n(i:i+1); %权的变化公式 end subplot (3,1,1) plot (t, xn ) %输出信号加噪声图形 title (带噪声原始信号) grid； subplot (3,1,2) plot (t , sn) title(原始正弦信号) subplot(3,1,3) plot (t , e) title(滤波结果)第4章 MATLAB在数字信号处理中的应用 5 5运行结果与分析运行结果与分析 收敛因子决定收敛速度及稳定性。图4.20设置的收敛因子为0.00026，图4.21设置的收敛因子为0.0026，图4.22设置的收

28、敛因子为0.026。由图可见，收敛因子的选择对滤波器的性能有很大影响。图4.20，滤波效果不错，但是收敛较慢；图4.21，收敛速度加快，但滤波效果不如图4.20；图4.22 收敛速度很快，但滤波效果太差。可见，选择合适的收敛因子对于滤波性能有很大影响。收敛因子的选择，其实就是在滤波性能及速度之间折中。 6思考题 (1)简述LMS 算法的原理。 (2)根据上述结果，收敛因子该如何选择？第4章 MATLAB在数字信号处理中的应用 图4.20 =0.00026 图4.21 =0.0026第4章 MATLAB在数字信号处理中的应用 图4.22 =0.026 时正弦加噪信号的滤波第4章 MATLAB在数

29、字信号处理中的应用 1实验目的 (1)了解频率采样定理是数字信号处理中的重要理论。 (2)掌握频率域采样会引起时域周期化的概念。 (3)掌握频率域采样定理及其对频率采样点数选择的 指导作用。 2实验原理 1) 频率域采样定理的要点 对信号x(n)的频谱在0,2上等间隔采样N点， 得到:4.10 4.10 频率采样定理频率采样定理第4章 MATLAB在数字信号处理中的应用 则N点的离散傅里叶反变换得到的序列就是原序 列x(n)以N为周期进行周期延拓后的主值序列，公式 为：4.10 4.10 频率采样定理频率采样定理第4章 MATLAB在数字信号处理中的应用 由上式可知，频率域采样点数N必须大于等

30、于时域离散信号的长度M(即NM)，才能使时域不产生混叠，则N点的离散傅里叶反变换得到的序列就是原序列项x(n)。 对比上面叙述的频域采样定理，得到一个有用的结论：频域采样，时域信号周期延拓。 2) 频域采样定理的验证 给定信号如下： 编写程序分别对频谱函数在区间0，2上等间隔采样32点和16点，得到：4.10 4.10 频率采样定理频率采样定理第4章 MATLAB在数字信号处理中的应用 3仿真思路 直接调用MATLAB函数fft计算再分别对进行32点和16点IFFT，得到：分别画出的幅度谱，并绘图显示的波形，进行对比和分析。 就得到频谱函数在0，2上的32点和16点频率采样。4.10 4.10

31、 频率采样定理频率采样定理第4章 MATLAB在数字信号处理中的应用 4程序代码 M=26;N=32;n=0:26; n1=0:M/2; n2=M/2+1:M; xa=n1+1; xb=27-n2; x=xa,xb; subplot(3,2,2); stem(n,x); title(三角波序列x(n) ) w=linspace(0,pi,1000);%设定频率向量X=x*exp(-j*n*w);subplot(3,2,1);plot(w/pi,abs(X);title(序列x(n)的连续幅度谱)X32k=fft(x,32);k1=0:31;subplot(3,2,3);stem(k1,abs(

32、X32k);axis(0,15,0,200);title(32 点频率采样幅度谱)x32n=ifft(X32k);nx32=0:31;subplot(3,2,4);stem(nx32,x32n);title(32 点IDFT)第4章 MATLAB在数字信号处理中的应用 X16k=fft(x,16); k2=0:15; subplot(3,2,5); stem(k2,abs(X16k); axis(0,8,0,200); title(16 点频率采样幅度谱) x16n=ifft(X16k); nx16=0:15; subplot(3,2,6); stem(nx16,x16n); axis(0,2

33、0,0,28); title(16 点IDFT)5运行结果与分析 由图4.23可知，在一定的条件下，可以由频域离第4章 MATLAB在数字信号处理中的应用 散采样恢复原来的信号，这条件就是：如果序列x(n)的长度为M，则只有当频域采样点数NM时，才可由频域采样X(k)恢复原序列x(n)，否则将产生时域混叠现象。 6思考题 (1) 频率采样定理和时域采样定理有什么不同？ (2) 频率采样定理在什么情况下能恢复原来的信号？4.10 4.10 频率采样定理频率采样定理第4章 MATLAB在数字信号处理中的应用第4章 MATLAB在数字信号处理中的应用图 4.23 对不同信号的采样 1实验目的 (1)

34、了解维纳滤波器的原理。 (2)了解维纳滤波器的MATLAB实现方法。 2实验原理 滤波技术是信号分析、处理技术的重要分支，无论是信号的获取、传输，还是信号的处理和交换都离不开滤波技术，它对信号安全可靠和有效灵活地传递是至关重要的。信号分析检测与处理的一个十分重要的内容就是从噪声中提取信号，实现这种功能的有效手段之一是设计具有最佳线性过滤特性的滤波器，当伴有噪声的信号通过这种滤波器时，它可以将信号尽4.11 4.11 维纳滤波算法维纳滤波算法第4章 MATLAB在数字信号处理中的应用 可能精确地重现或对信号做出尽可能精确的估计，而对所伴随噪声进行最大限度的抑制。维纳滤波器就是这种滤波器的典型代表

35、之一。 1) 维纳滤波概述 维纳(Wiener)是用来解决从噪声中提取信号的一种滤波方法。这种线性滤波问题可以看成是一种估计问题或一种线性估计问题。 如果一个线性系统的单位样本响应为h(n)，当输入一个随机信号x(n)，且4.11 4.11 维纳滤波算法维纳滤波算法第4章 MATLAB在数字信号处理中的应用其中x(n)表示信号，v(n)表示噪声，则输出y(n)为： 希望x(n)通过线性系统h(n)后得到的y(n)尽量接近于s(n)，因此称y(n)为s(n)的估计值，即： 则维纳滤波器的输入-输出关系可用图4.24表示。4.11 4.11 维纳滤波算法维纳滤波算法第4章 MATLAB在数字信号处

36、理中的应用 图4.24 维纳滤波器的输入-输出关系 用h(n)进行过滤问题实际上是一种统计估计问题。一般地，从当前的和过去的观察值x(n)，x(n 1)， x(n 2)，估计当前的信号值称为过滤或滤波；从过去的观察值，估计当前的或者将来的信号值称为外推或预测；从过去的观察值，估计过去的 信号值称为平滑或内插。因此维纳滤波器又常常被称为最佳线性过滤与预测或线性最优估计。这里所谓的最佳与最优是以最小均方误差为准则的。 用e(n)表示信号的真实值和估计值之间的误差，显然e(n)可能是正值，也可能是负值，并且它是一个随机变量。因此用它的均方误差来表达误差是合理的，均方误差最小即它的平方的统计期望最小。

37、采用最小均方误差准则作为最佳过滤准则的原因还在于它的理论分析比较简单，不要求对概率的描述。 2) 维纳-霍夫方程的求解4.11 4.11 维纳滤波算法维纳滤波算法第4章 MATLAB在数字信号处理中的应用 为了按最小均方误差准则来确定维纳滤波器的冲激响应h(n)，令(n) 对h(j)的导数等于零，即可得维纳-霍夫方程：4.11 4.11 维纳滤波算法维纳滤波算法第4章 MATLAB在数字信号处理中的应用其中： 维纳-霍夫方程右端的求和范围没有具体标明，实际上有3种情况。 有限冲激响应(FIR)维纳滤波器，i从0到N1取得有限个整数值。 非因果无限冲激响应(非因果IIR)维纳滤波器，i从到+取所

38、有整数值。 因果无限冲激响应(因果IIR)维纳滤波器，i从0到+取正整数值。 上述3种情况下维纳-霍夫方程的解法不同，这里只描述FIR维纳滤波器的求解。 设滤波器冲激响应序列的长度为N，冲激响应矢量为：4.11 4.11 维纳滤波算法维纳滤波算法第4章 MATLAB在数字信号处理中的应用滤波器输入数据矢量为：则滤波器的输出为：这样，上述所示的维纳滤波器的标准方程可写成： 其中，P是s(n)与x(n)的互相关函数，它是一个N维列矢量；R是x(n)的自相关函数，是N阶方阵。 利用求逆矩阵的方法直接求解上式，得：4.11 4.11 维纳滤波算法维纳滤波算法第4章 MATLAB在数字信号处理中的应用o

39、pt表示“最佳”，这就是FIR维纳滤波器的冲激响应。 3仿真思路 假设一个点目标在 x、y平面上绕单位圆做圆周运动，由于外界干扰，其运动轨迹发生了偏移。其中，x方向的干扰为均值为0、方差为0.05的高斯噪声；y方向的干扰为均值为0、方差为0.06的高斯噪声。 (1)产生满足要求的x方向和y方向随机噪声样本500 个。 (2)明确期望信号和观测信号。 (3)试设计一个FIR维纳滤波器，确定最佳传递函数并用该滤波器处理观测信号，得到其最佳估计。 (4)分别绘制出x方向和y 方向的期望信号、噪声信号、观测信号、滤波后信号、最小均方误差信号的曲线图。 (5)在同一幅图中绘制出期望信号、观测信号和滤波后

40、点目标的运动轨迹。 4程序代码4.11 4.11 维纳滤波算法维纳滤波算法第4章 MATLAB在数字信号处理中的应用 clear; clf; sita=0:pi/249.5:2*pi; xnoise=sqrt(0.05)*randn(1,500);%产生x轴方向噪声 ynoise=sqrt(0.06)*randn(1,500);%产生y轴方向噪声 x=cos(sita)+xnoise;%产生x轴方向观测信号 y=sin(sita)+ynoise;%产生y轴方向观测信号 %产生维纳滤波中x方向上观测信号的自相关矩阵 rxx=xcorr(x); for i=1:100 for j=1:100 mr

41、xx(i,j)=rxx(500-i+j); endend第4章 MATLAB在数字信号处理中的应用 clear; clf; sita=0:pi/249.5:2*pi; xnoise=sqrt(0.05)*randn(1,500);%产生x轴方向噪声 ynoise=sqrt(0.06)*randn(1,500);%产生y轴方向噪声 x=cos(sita)+xnoise;%产生x轴方向观测信号 y=sin(sita)+ynoise;%产生y轴方向观测信号 %产生维纳滤波中x方向上观测信号的自相关矩阵 rxx=xcorr(x); for i=1:100 for j=1:100 mrxx(i,j)=r

42、xx(500-i+j); endend第4章 MATLAB在数字信号处理中的应用 xd=cos(sita);%产生维纳滤波中x 方向上观测信号与期望信号的互相关矩阵 rxd=xcorr(x,xd); for i=1:100 mrxd(i)=rxd(499+i); end hoptx=inv(mrxx)*mrxd; %由维纳-霍夫方程得到的x 、 方向上的滤波器最优解 fx=conv(x,hoptx); %滤波后x方向上的输出 nx=sum(abs(xd).2); eminx=nx-mrxd*hoptx; %x方向上最小均方误差%产生维纳滤波中y 方向上观测信号的自相关矩阵 ryy=xcorr(

43、y);第4章 MATLAB在数字信号处理中的应用 for i=1:100 for j=1:100 mryy(i,j)=ryy(500-i+j); end end yd=sin(sita);%产生维纳滤波中y方向上观测信号与期望信号的互相关矩阵 ryd=xcorr(y,yd); for i=1:100 mryd(i)=ryd(499+i); end hopty=inv(mryy)*mryd; %由维纳-霍夫方程得到的y方向上的滤波器最优解 第4章 MATLAB在数字信号处理中的应用fy=conv(y,hopty); %滤波后y方向上的输出ny=sum(abs(yd).2);eminy=ny-mr

44、yd*hopty; %y方向上最小均方误差subplot(2,4,1)plot(xd);title(x 方向期望信号);subplot(2,4,2)plot(xnoise);title(x 方向噪声信号);subplot(2,4,3)plot(x);title(x 方向观测信号);subplot(2,4,4)n=0:500;plot(n,eminx);第4章 MATLAB在数字信号处理中的应用title(x 方向最小均方误差);subplot(2,4,5)plot(yd);title(y 方向期望信号);subplot(2,4,6)plot(ynoise);title(y 方向噪声信号);su

45、bplot(2,4,7)plot(y);title(y 方向观测信号);subplot(2,4,8)plot(n,eminy);title(y 方向最小均方误差);figure;plot(xd,yd,k);第4章 MATLAB在数字信号处理中的应用 hold on; plot(x,y,b:); hold on; plot(fx,fy,g-); title(最终结果); 5运行结果与分析 x轴和y轴信号、噪声及最小均方误差如图4.25 所示。滤波后得到的信号与原始信号和噪声信号的对比如图4.26所示，滤波后的结果与期望信号还是很接近的，整体上达到了最优滤波的效果。 6思考题 (1)维纳滤波原理其

46、实是一个最优估计的问题。简述上述程序实现维纳滤波的过程。第4章 MATLAB在数字信号处理中的应用第4章 MATLAB在数字信号处理中的应用 图4.25 信号、噪声及最小均方误差滤 图 4.26 最终结果 (2)上例中，估计信号与原始信号之间还有一种事实上 的误差存在，为什么？如何使这种误差最小化？第4章 MATLAB在数字信号处理中的应用 1实验目的 (1)了解卡尔曼滤波器的原理。 (2)学习卡尔曼滤波器的仿真实现方法。 2实验原理 卡尔曼滤波器是一个“ 最优化自回归数据处理算法。对于解决大部分的问题，它是最优、效率最高甚至是最有用的。其广泛应用已经超过30年，包括机器人导航、控制、传感器数

47、据融合甚至在军事方面的雷达系统以及导弹追踪等。近年来更被应用于计算机图像处理，例如，面部识别、图像分割、图像边缘检测等方面。4.12 4.12 卡尔曼滤波算法卡尔曼滤波算法第4章 MATLAB在数字信号处理中的应用 卡尔曼滤波原理卡尔曼滤波原理 首先要引入一个离散控制过程的系统，该系统可用一个线性随机微分方程来描述： X(k)=AX(k1)+BU(k)+W(k) 再加上系统的测量值： Z(k)=H X(k)+V(k) 上两式子中，X(k)是k时刻的系统状态，U(k)是k 时刻对系统的控制量。A和B是系统参数，对于多模型系统，它们为矩阵。Z(k)是k时刻的测量值， H是测量系统的参数，对于多测量

48、系统，H为矩阵。W(k)和V(k)分别表示过程和测量的噪声。4.12 4.12 卡尔曼滤波算法卡尔曼滤波算法第4章 MATLAB在数字信号处理中的应用 它们被假设成高斯白噪声，其协方差分别是Q、R这 里假设它们不随系统状态变化而变化。 由于满足上面的条件(线性随机微分系统，过程和测量都是高斯白噪声)，卡尔曼滤波器是最优的信息处理器。下面来估算系统的最优化输出。 首先利用系统的过程模型预测下一个状态的系统。假设现在的系统状态是k，根据系统的模型，可以基于系统的上一状态而预测出现在状态： X(k|k-1)=AX(k-1|k-1)+BU(k) 上式中，X(k|k-1)是利用上一个状态预测的结果，X(k-1|k-1)是上一个状态最优的结果，U(k)为现在状4.12 4.12 卡尔曼滤波算法卡尔曼滤波算法第4章 MATLAB在数字信号处理中的应用 态的控制量，如果没有控制量，它可以为0。 到现在为止，系统结果已经更新了，可是对应于X(k|k-1)的协方差还没更新。用P表示协方差： P(k|k1)=A P(k1|k1) A+Q 上式中，P(k|k1)是X(k|k1)对应的协方差，P(k1|k1)是X(k1|k1)对应的协方差，A表示A的转置矩阵，Q是系统过程