维纳滤波器-维纳霍夫方程课程设计.doc

电气工程学院信号处理课程设计报告设计题目： 维纳滤波器 专 业： 生物医学工程 指导教师: 叶立夏 学生姓名: 叶立夏 学 号: 20104320144 起迄日期: 2013年12月20日2014年1月15日如果有做课程设计的同学不懂的，可以联系我！Qq：745581242目 录前言3 1 设计任务及指标3 1.1 课程设计的内容和要求4 1.2 对课程设计成果的要求42 设计思想42.1 概述42.2 主要仪器设备及耗材53 课程设计的组成部分53.1具体操作54 实验分析94.1 原始图像显示94.2 噪声的强度对维纳滤波器的影响13 4.3 阶数对滤波效果的影响174.4 数据长度对维纳滤波的影响21 5 设计总结265.1思考题265.2实验心得266 主要参考文献 27 附录27前言去除信号中的噪声影响是信号处理中的一个重要内容，而滤波则是实现这一功能的重要手段之一。滤波器可以分为两类，及经典滤波器和现代滤波器经典滤波器是假定输入信号中嘚瑟有用成分和希望去除的成分各自占有不同的频带。当输入信号通过一个滤波器是可将欲去除的成分有效的去除，如果信号和噪声的频谱相互重叠，那么经典滤波器将无能为力。现代滤波器理论研究的主要内容是从含有噪声的数据记录中估计出信号的某些特征灬信号本身。一旦信号被估计出，那么估计出的信号的信噪比将比原信号的高。现代滤波器把信号和噪声都视为随机信号，利用它们的统计特征导出一套最佳的估值算法，然后用硬件或软件予以实现。现代滤波器理论源于维纳在20世纪40年代及其以后的工作，因此维纳滤波器便是这一类滤波器的典型代表。维纳滤波器，也是最小平方滤波器，其基本思路为：设计一个滤波器，使其与输入信号滤波后的输出与期望输出在最小平方意义下的最佳逼近。寻求最小均方误差的实质就是解维纳-霍夫方程。1.设计任务及指标1.1课程设计的内容和要求（包括原始数据、技术要求、工作要求等）：本设计的目的是产生用于信号滤波的维纳霍夫方程。并且要求调节该滤波器的参数使该滤波器能够最好的还原原始波形，以适应不同原始信号都能够被提取出来。设计要求：1.已知信号的自相关函数和噪声的能量，编写程序求解维纳-霍夫方程，寻找最优滤波器。2.编写程序仿真信号，噪声和观察波形，然后把观察信号通过滤波器得到的信号估计与原始信号比较，观察是否达到了去噪的目的。3.选择不同信号（仿真信号，实际采集的心电，脑电信号），人工添加噪声，调整噪声的相对强度，观察滤波效果。1.2. 对课程设计成果的要求：利用MATLAB软件编写程序。根据维纳滤波的原理，通过调节参数获得效果最佳的滤波器，用正弦衰减信号作为测试信号观察滤波效果，再将心电信号和脑电信号作为输入，观察滤波结果；调节滤波器的参数（噪声强度，阈值，数据长度）观察参数对滤波器的影响。2. 设计思想2.1概述根据正交原理可以推导出维纳-霍夫方程，满足该方程的滤波器输出信号的估计值与信号在最小均方误差意义下最接近。 （2-1）根据滤波器的形式，维纳滤波器可以分为三种情况：非因果IIR型，因果IIR型，FIR型，对于实时性有要求的情况下用后两种形式。图2.1 维纳滤波器对于FIR型维纳滤波器，维纳-霍夫方程的形式为 （2-2）或者写成矩阵形式 （2-3）其中 （2-4） （2-5） （2-6）这样，如果信号和噪声的二阶统计特性已知，则易求解 （2-7）维纳滤波的均方误差是 （2-8）2.2 主要仪器设备及耗材MATLAB软件，计算机3.课程设计的组成部分3.1具体操作3.1.1.建立一个以WH命名的函数（源程序上半部分）function h,e = WH(Rss,Rww,M)% 求解维纳-霍夫方程的函数，其中M为信号的长度e1 = 10; %给e1赋初值e0 = 0; %给e0赋初值N = 0; % 给N 赋初值% 以下循环的目的是找出FIR滤波器合适的阶数% 判据是当阶数增加而均方误差没有明显下降时，则认为阶数足够while abs(e0-e1)1e-6 % e1和e0不够接近则循环，确定阈值 N = N+1; %使N递增，使其逐步靠近需要满足的条件 e0 = e1; Rxs = Rss(M:(M+N-1); Rxx = Rww(M:(M+N-1)+Rss(M:(M+N-1);%输入信号总的自相关函数 R_xx = zeros(N); %产生一个N*N的全零矩阵作为仿真信号的自相关函数 for j = 1:N for n = 1:N R_xx(j,n) = Rxx(abs(j-n)+1); end end h = inv (R_xx)*Rxs; %求解互相关函数转置的逆矩阵 %e1 = Rss(M)-h*Rxs; %endN % 显示N的最终值e = e1;用于利用线性最小均方误差求解维纳霍夫方程。3.1.2.建立一个主程序，运行主程序将输入信号设置为测试信号即正弦衰减信号（源程序的下半部分）噪声为强度为0.4的白噪声。% 主程序clear; clc;M = input(信号的长度 M = );n = 1:M;s = exp(-0.002*n).*sin(pi*n/50); % 仿真信号，可以自己生成，任意形式% load ecgdata; % 实际心电信号% s = ecgdata (1:M);% load eedata; % 实际脑电信号% s = eegdata (1:M);w = 0.4*randn(1,M); % 白噪声，系数代表噪声相对强度x = s+w; % 仿真信号Rss = xcorr(s,s); %估计信号自相关函数Rww = xcorr(w,w); % 估计噪声自相关函数h,e = WH(Rss,Rww,M);ss = filter(h,1,x); %用维纳滤波器滤波figure;subplot(2,2,1);plot(n,s);title(信号);subplot(2,2,2);plot(n,w);title(噪声);subplot(2,2,3);plot(n,x);title(观测值);subplot(2,2,4);plot(n,ss);title(信号估计);figure;plot(n,ss-s);title(估计误差);error=mean(s-ss).2)3.1.3.调节输入信号为心电信号，心电信号的获取是从专门的心电数据网站/cgi-bin/atm/ATM提取的，调节信号源以及心电信号的输出方式显示为数据txt的方式，将数据复制存档成txt格式文件，然后复制txt文件到Excel文件中，变成3列，运用列变化提取一列就成为一维的心电信号文件最后Excel截图显示为：图3.1 网站心电数据显示图图3. 2 存档为txt格式的文件图3. 3 存为Excel截图 图3. 4 存为txt一维心电3.1.4.脑电数据的提取方法跟心电数据的基本一样来源网站/cgi-bin/atm/ATM信号源的选取为EEG信号，其他部分与心电相似。图3. 5 网站脑电信号来源注：心电信号和脑电信号都是以txt的格式存取的，其是通过load函数与主程序的输入函数相连接的。具体格式为q=load(C:UserstangDesktopeeg.txt);s = q(1:M);4 实验分析改变主程序的参数，调节参数使其符合最优滤波器的要求，分析各个参数对滤波器的影响，运行结果分析分为如下几个：4.1、原始图像显示：w为强度为0.4白噪声 M=1000的信号图4.1.1、S为指数正弦衰减信号s = exp(-0.002*n).*sin(pi*n/50)，w为强度为0.4白噪声 M=1000 源程序如下：clear; clc;M = input(信号的长度 M = );n = 1:M;t=10*n;s = exp(-0.002*n).*sin(pi*n/50); w = 0.1*randn(1,M); % 白噪声，系数代表噪声相对强度x = s+w; % 仿真信号Rss = xcorr(s,s); %估计信号自相关函数Rww = xcorr(w,w); % 估计噪声自相关函数h,e = WH(Rss,Rww,M);ss = filter(h,1,x); %用维纳滤波器滤波figure;subplot(2,2,1);plot(t,s);title(信号);xlabel(x轴单位：t/ms,color,b) ylabel(y轴单位：mV,color,b) subplot(2,2,2);plot(t,w);title(噪声);xlabel(x轴单位：t/ms,color,b) ylabel(y轴单位：mV,color,b) subplot(2,2,3);plot(t,x);title(观测值);xlabel(x轴单位：t/ms,color,b) ylabel(y轴单位：mV,color,b) subplot(2,2,4);plot(t,ss);title(信号估计);xlabel(x轴单位：t/ms,color,b) ylabel(y轴单位：mV,color,b) figure;plot(t,ss-s);title(估计误差);error=mean(s-ss).2) w = 0.1*randn(1,M); % 白噪声，系数代表噪声相对强度x = s+w; % 仿真信号Rss = xcorr(s,s); %估计信号自相关函数Rww = xcorr(w,w); % 估计噪声自相关函数h,e = WH(Rss,Rww,M);ss = filter(h,1,x); %用维纳滤波器滤波figure;subplot(2,2,1);plot(t,s);title(信号);xlabel(x轴单位：t/ms,color,b) ylabel(y轴单位：mV,color,b) subplot(2,2,2);plot(t,w);title(噪声);xlabel(x轴单位：t/ms,color,b) ylabel(y轴单位：mV,color,b) subplot(2,2,3);plot(t,x);title(观测值);xlabel(x轴单位：t/ms,color,b) ylabel(y轴单位：mV,color,b) subplot(2,2,4);plot(t,ss);title(信号估计);xlabel(x轴单位：t/ms,color,b) ylabel(y轴单位：mV,color,b) figure;plot(t,ss-s);title(估计误差);error=mean(s-ss).2)实验输出图像图4. 1 测试信号输出图形显示信号的长度 M = 1000 阶数 N =38 最小均方误差error =0.016结果分析：观察上图，比较信号与估计信号，我们可以看出噪声对信号有很大的影响，虽然还看得出有用信号的轮廓，但是噪声的影响还是没能消除，滤波效果还算可以。4.1.2、S为心电信号，w为强度为0.4的白噪声 ，M=500,采样频率125Hz.（心电和脑电信号的横坐标均为时间，纵坐标均为幅值，且都有单位）实验输出图像图4. 2 心电信号滤波结果显示信号的长度 M = 500 阶数N =46 最小均方误差error =0.0448 简要分析：从上面的图可以看出滤波效果不好，比起指数正弦衰减信号而言，估计信号直接很难辨别，但是虽然滤波效果不好但却估计出了峰值(R波)，比较估计误差图和估计信号，我们不难看出在估计信号产生波峰的时候，估计误差图中对应的产生了与其相反的信号，这就是噪声对心电信号的影响。4.1.3、S为脑电信号，w为强度为0.1的白噪声, M=1024,采样频率250Hz.实验输出图像图4. 3 脑电信号滤波结果显示信号的长度 M = 1024 阶数N =5 最小均方误差error = 0.0100结果分析：在此次运行的结果图中可以看出脑电信号的滤波效果是特别差的了，几乎与原信号和噪声信号相似，这主要是因为维纳滤波主要处理的事有规律的信号，对周期信号处理的效果比较好。而在此次运行的信号中，信号与噪声信号相似，所以信号估计滤波效果不好。换句话说就是：区分度越大，滤波效果越好。4.2、噪声的强度对维纳滤波器的影响：M=1024，改变噪声的强度为0.2、1、2的各信号分析及图4.2.1、S为指数正弦衰减信号，w为强度为0.2、1、2白噪声 M=1024当W为强度为0.2的白噪声时，实验结果图4.4 当W为强度为0.2的白噪声时，实验结果信号的长度 M = 1024 阶数N =34 最小均方误差error =0.0017当W为强度为1的白噪声时，实验结果 图4.5 当W为强度为1的白噪声时，实验结果信号的长度 M = 1024 阶数N = 236 最小均方误差error =0.0291当W为强度为2的白噪声时，实验结果图4.6 当W为强度为2的白噪声时，实验结果信号的长度 M = 1024 阶数N = 74 最小均方误差error = 0.0015结果分析：对上图进行观察分析可以看出：相比之下均值为0.2的与原始信号最相似，均值为2的滤波效果最差，此外我们还可以看出N的值和最小均方差都在增大，所以我们大致可以得出以下结论：均值越小，最小均方误差越小，估计信号越准确。4.2.2、S为心电信号，w为强度为0.2、1、2白噪声 M=500当S为心电信号，W为强度为0.2的白噪声时，实验结果图4.7 当S为心电信号，W为强度为0.2的白噪声时，实验结果信号的长度 M = 500 阶数N =290 最小均方误差error = 0.0144当S为心电信号，W为强度为1的白噪声时，实验结果图4.8 当S为心电信号，W为强度为1的白噪声时，实验结果信号的长度 M = 500 阶数N = 11 最小均方误差error = 0.0686当S为心电信号，W为强度为2的白噪声时，实验结果图4.9 当S为心电信号，W为强度为2的白噪声时，实验结果信号的长度 M = 500 阶数N = 45 最小均方误差error =0.0682结果分析：在心电信号改变噪声强度的情况下，观察其变化规律，在均值增大的同时，噪声对心电信号的R,Q,P,S和T波的影响度是很大的，使得滤波效果不明显。并且均值，最小均方差都在影响着估计信号的好坏。4.2.3、S为脑电信号，w为强度为20、50、100白噪声 M=1024当S为脑电信号，W为强度为20的白噪声时，实验结果图4.10 当S为脑电信号，W为强度为20的白噪声时，实验结果信号的长度 M = 1024 阶数N = 1011 最小均方误差error = 191.2949当S为脑电信号，W为强度为50白噪声时，实验结果图4.11 当S为脑电信号，W为强度为50白噪声时，实验结果信号的长度 M = 1024 阶数N =526 最小均方误差error =520.1239当S为脑电信号，W为强度为100白噪声时，实验结果图4.12 当S为脑电信号，W为强度为100白噪声时，实验结果信号的长度 M = 1024 阶数N = 526 最小均方误差error = 520.1639结果分析：与之前的滤波结果大同小异，脑电信号的滤波结果着实不好，着主要是因为脑电信号紊乱，不过这也反映出了实际情况，我们的大脑是时时在工作的，他的工作精细而复杂。噪声信号对脑电信号的影响非常大。4.3、阶数对滤波效果的影响：w为强度为0.4白噪声 M=1024，分别指数衰减信号，心电信号为例 【注：阶数的修改是通过修改WH功能函数实现的】4.3.1指数衰减信号，阈值为 、：S为指数衰减信号，W为强度0.4白噪声，阈值为时实验结果图4.13 S为指数衰减信号，阈值为时实验结果信号的长度 M = 1024 阶数N = 136 最小均方误差error = 0.0117S为指数衰减信号，W为强度0.4白噪声，阈值为时实验结果图4.14 S为指数衰减信号，阈值为时实验结果信号的长度 M = 1024 阶数N =141 最小均方误差error =0.0094S为指数衰减信号，W为强度0.4白噪声，阈值为时实验结果图4.15 S为指数衰减信号，阈值为时实验结果信号的长度 M = 1024 阶数N=65 最小均方误差error =0.0101结果分析：改变阈值的大小，得到以上三幅图，由此可以看出随着阈值的减小，图像滤波的效果也随之增强。4.3.2、心电信号，数据长度M=500，阈值为 、：S为心电信号，W为强度0.4白噪声，阈值为时图4.16 S为心电信号，阈值为时信号的长度 M = 500 阶数N =110 最小均方误差error =0.0363S为心电信号，W为强度0.4白噪声，阈值为时实验结果图4.17 S为心电信号，阈值为时实验结果信号的长度 M = 500 阶数N =226 最小均方误差error = 0.0365结果分析：由误差信号可知，数据的长度一定时，可以通过改变滤波器的阶数来减小最小误差，从而改善滤波器的滤波效果，除此之外，经过多次运行程序后可知，当阶数达到某个值时，误差的改善不再明显，因此滤波器的阶数对实验结果有很大影响，增加滤波器的阶数可以提高滤波器的性能。4.3.3脑电信号，数据长度M=1024，阈值为 、：S为脑电信号，W为强度50白噪声，阈值为时实验结果图4.18 S为脑电信号，阈值为时实验结果信号的长度 M = 1024 阶数N =970 最小均方误差error =533.2368S为脑电信号，W为强度50白噪声，阈值为时实验结果图4.19 S为脑电信号，阈值为时实验结果信号的长度 M = 1024 阶数N =523 最小均方误差error = 521.34154.4、数据长度对维纳滤波的影响：改变数据长度，观察噪声对信号的影响4.4.1指数衰减信号w为强度为0.4白噪声 ，阈值为M分别为1000、3000、4000指数衰减信号 信号长度M=1000实验结果图4.20 指数衰减信号 信号长度M=1000实验结果信号的长度 M = 1000 阶数N =166 最小均方误差error =0.0078指数衰减信号 信号长度M=3000实验结果图4.21 指数衰减信号 信号长度M=3000实验结果信号的长度 M = 3000 阶数N = 190 最小均方误差error =0.0060指数衰减信号 信号长度M=4000实验结果图4.22 指数衰减信号 信号长度M=4000实验结果信号的长度 M = 4000 阶数N =848 最小均方误差error =0.0045结果分析：对比观察发现，信号长度越大，N值也逐渐在增大，最小均方差也逐渐缓慢减小，因此信号估计图就越接近原图像，着主要是因为M越大，样本点数也就越大，原始信号采集量就越大，出来的结果自然也比较准确。4.4.2心电信号w为强度为0.4白噪声 阈值为M分别为1000、3000、4000(在源程序中加load ecgdata; % 实际心电信号s=eegdata,eegdata,eegdata;s = s (1:M);实现周期延拓)心电信号，数据长度M=1000时的实验结果图4.23 心电信号，数据长度M=1000时的实验结果信号的长度 M = 1000 阶数N =95 最小均方误差error = 0.0432心电信号，数据长度M=3000时的实验结果 图4.24 心电信号，数据长度M=3000时的实验结果信号的长度 M = 3000 阶数N =134 最小均方误差error = 0.0414心电信号，数据长度M=4000时的实验结果图4.25 心电信号，数据长度M=4000时的实验结果信号的长度 M = 4000 阶数N =38 最小均方误差error =0.0565结果分析：因为心电信号及脑电信号有局限所以要通过周期延拓才能取到长度比1024大的，观察上图，不难发现样本点数越大越好，这样就可以更精确的显示结果。4.4.3脑电信号w为强度为50白噪声 M分别为1000、3000、4000脑电信号，数据长度M=1000时实验结果图4.26 脑电信号，数据长度M=1000时实验结果信号的长度 M = 1000 阶数N =221 最小均方误差error = 623.2513脑电信号，数据长度M=3000时实验结果图4.27 脑电信号，数据长度M=3000时实验结果信号的长度 M = 3000 阶数N =482 最小均方误差error =684.2442脑电信号，数据长度M=4000时实验结果图4.28 脑电信号，数据长度M=4000时实验结果信号的长度 M = 4000 阶数N =612 最小均方误差error = 873.7261结果分析：仿心电信号，脑电信号也具有相似的结果，只是脑电信号的原信号没心电信号清晰，才导致观察起来有障碍，但还是可以看出一定的规律性。对改变信号长度的小结：由对指数衰减信号，心电信号，脑电信号改变长度的观察，当滤波器的阶数一定，观测数据的长度增加时，可以减小估计信号与期望信号之间的差值，因此，观测数据的长度对实验有着重要的影响，即增加信号样本个数可以提高维纳滤波的性能。5设计总结5.1思考题观察实验结果，对于几种不同的信号，维纳滤波是否都取得了较好的效果？如果效果不好，试分析原因。答：观察实验结果，对于指数衰减信号，心电信号，脑电信号这三种信号而言，维纳滤波仅对指数衰减信号的滤波效果更好一些，对心电信号次之，对脑电信号的滤波效果不好，这些主要是受输入信号（输入信号的区分度越大，滤波效果越好）；噪声强度（噪声强度越小，估计信号越准确）；最小均方误差（最小均方误差越小，滤波效果越好）；阶数（提高阶数，滤波效果越好）；信号的长度（增大信号长度使得滤波效果更精确）和维纳滤波器本身的局限性（只适用于平稳随机信号）等的影响。5.2实验心得在做本次实验之前，我对这次试验的要求以及原理可以说是一点都了解，茫然的很，知道这次试验要用matlab软件实现，自己就慌了，当初学习这个软件的时候三心二意，学得乱七八糟。对于维纳滤波，也了解甚少，当初学得也忘记的一干二净了。多亏了同学的帮忙耐心教导，讲解了浅显的matlab软件的应用，和试验步骤与基本方法，自己按照网上的方法一步一步将软件安装好，总算有了个开头。然后再去图书馆借了相关的书籍，翻阅出原来学过的专业书生物医学信号处理，上网查询matlab相关的实现维纳滤波的程序。从了解维纳滤波的原理，到以及matlab的实现方法，细细琢磨，才将此次设计的原理看懂一点。通过此次试验设计，了解要用严谨的态度来对待事情，不能马马虎虎，讲究就好。实验中，学到不少的新的知识，各个方面要亲历亲为，不懂的东西查资料或者问同学，知道弄懂为止。维纳滤波是一种最基本的提取信号，去除噪声的方法，适用于从噪声中分离出有用的信号波形。其优点是适应面较广，无论是平稳随机是连续的还是离散的，是标量的还是向量的，都可应用。缺点是在要求得到半无限区间内的全部观察数据的条件很难满足，同时也不能用于噪声为非平稳的随机过程的情况，对于向量情况应用也不方便。因此，维纳滤波在实际问题中应用不多。因为没有在学校的原因，没有能与邱老师您交流，给您说声对不起，自己也感到遗憾错过了您的指导。6.主要参考文献1刘浩.MATLAB R2001a完全自学一点通M.北京：电子工业出版社，20132张志涌。掌握与精通matlab。北京：北京航空大学出版社，19983 唐向宏，岳恒立，郑雪峰。基于matlab及在电子信息类课程中的应用。北京：电子工业出版社，2006附录WH功能函数：function h,e = WH(Rss,Rww,M)% 求解维纳-霍夫方程的函数，其中M为信号的长度e1 = 10; %给e1赋初值e0 = 0; %给e0赋初值N = 0; % 给N 赋初值% 以下循环的目的是找出FIR滤波器合适的阶数% 判据是当阶数增加而均方误差没有明显下降时，则认为阶数足够while abs(e0-e1)1e-6 % e1和e0不够接近则循环，确定阈值 N = N+1; %使N递增，使其逐步靠近需要满足的条件 e0 = e1; Rxs = Rss(M:(M+N-1); Rxx = Rww(M:(M+N-1)+Rss(M:(M+N-1);%输入信号总的自相关函数 R_xx = zeros(N); %产生一个N*N的全零矩阵作为仿真信号的自相关函数 for j = 1:N for n = 1:N R_xx(j,n) = Rxx(abs(j-n)+1); end end h = inv (R_xx)*Rxs; %求解互相关函数转置的逆矩阵 %e1 = Rss(M)-h*Rxs; %endN % 显示N的最终值e = e1;测试信号霍夫方程源程序：% 主程序clear; clc;M = input(信号的长度 M = );n = 1:M; %确定数据的输入是1到Ms = exp(-0.002*n).*sin(pi*n/50); % 仿真信号，可以自己生成，任意形式%q=load(D:eeg.txt); %s = q(1:M);% load eedata; % 实际脑电信号% s = eegdata (1:M);w = 0.8*randn(1,M); % 白噪声，系数代表噪声相对强度x = s+w; % 仿真信号Rss = xcorr(s,s); %估计信号自相关函数Rww = xcorr(w,w); % 估计噪声自相关函数h,e = WH(Rss,Rww,M);%自定义的功能函数，用来求解维纳霍夫ss = filter(h,1,x); %用维纳滤波器滤波,输入x为滤波前序列% ss为滤波结果序列，h/1提供滤波器系数，h为分子，1为分母figure;subplot(2,2,1);plot(n,s);title(信号);subplot(2,2,2);plot(n,w);title(噪声);subplot(2,2,3);plot(n,s+w);title(观测值);subplot(2,2,4);plot(n,ss);title(信号估计);figure; % 控制窗口数量，在下一个窗口里面plot(n,ss-s);title(估计误差);error=mean(s-ss).2) %求最小均方误差脑电图霍夫方程源程序：% 主程序clear; clc;M = input(信号的长度 M = );n = 1:M;%s = exp(-0.002*n).*sin(pi*n/50); % 仿真信号，可以自己生成，任意形式% load ecgdata; % 实际脑电信号%s =ecgdata(1:M);q=load(D:eeg.txt); s = q(1:M);w = 0.4*randn(1,M); % 白噪声，系数代表噪声相对强度x = s+w; % 仿真信号Rss =