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

电气工程学院信号处理课程设计报告设计题目：维纳滤波器专业：生物医学工程指导教师:叶立夏学生姓名:叶立夏学号:20104320144起迄日期:2013年12月20日—2014年1月15日如果有做课程设计的同学不懂的，可以联系我！Qq：745581242目录前言………………………………………………31设计任务及指标…………………………………31.1课程设计的内容和要求………………41.2对课程设计成果的要求………………42设计思想…………………………………………42.1概述……………………………………42.2主要仪器设备及耗材…………………53课程设计的组成部分……………………………53.1具体操作…………………………………54实验分析…………………………………………94.1原始图像显示…………………………94.2噪声的强度对维纳滤波器的影响……134.3阶数对滤波效果的影响………………174.4数据长度对维纳滤波的影响…………215设计总结………………………………………265.1思考题…………………………………265.2实验心得………………………………266主要参考文献…………………………………27附录……………………………………………27前言去除信号中的噪声影响是信号处理中的一个重要内容，而滤波则是实现这一功能的重要手段之一。滤波器可以分为两类，及经典滤波器和现代滤波器经典滤波器是假定输入信号中嘚瑟有用成分和希望去除的成分各自占有不同的频带。当输入信号通过一个滤波器是可将欲去除的成分有效的去除，如果信号和噪声的频谱相互重叠，那么经典滤波器将无能为力。现代滤波器理论研究的主要内容是从含有噪声的数据记录中估计出信号的某些特征灬信号本身。一旦信号被估计出，那么估计出的信号的信噪比将比原信号的高。现代滤波器把信号和噪声都视为随机信号，利用它们的统计特征导出一套最佳的估值算法，然后用硬件或软件予以实现。现代滤波器理论源于维纳在20世纪40年代及其以后的工作，因此维纳滤波器便是这一类滤波器的典型代表。维纳滤波器，也是最小平方滤波器，其基本思路为：设计一个滤波器，使其与输入信号滤波后的输出与期望输出在最小平方意义下的最佳逼近。寻求最小均方误差的实质就是解维纳-霍夫方程。1.设计任务及指标1.1课程设计的内容和要求（包括原始数据、技术要求、工作要求等）：本设计的目的是产生用于信号滤波的维纳霍夫方程。并且要求调节该滤波器的参数使该滤波器能够最好的还原原始波形，以适应不同原始信号都能够被提取出来。设计要求：1.已知信号的自相关函数和噪声的能量，编写程序求解维纳-霍夫方程，寻找最优滤波器。2.编写程序仿真信号，噪声和观察波形，然后把观察信号通过滤波器得到的信号估计与原始信号比较，观察是否达到了去噪的目的。3.选择不同信号（仿真信号，实际采集的心电，脑电信号），人工添加噪声，调整噪声的相对强度，观察滤波效果。1.2.对课程设计成果的要求：利用MATLAB软件编写程序。根据维纳滤波的原理，通过调节参数获得效果最佳的滤波器，用正弦衰减信号作为测试信号观察滤波效果，再将心电信号和脑电信号作为输入，观察滤波结果；调节滤波器的参数（噪声强度，阈值，数据长度）观察参数对滤波器的影响。2.设计思想2.1概述根据正交原理可以推导出维纳-霍夫方程，满足该方程的滤波器输出信号的估计值与信号在最小均方误差意义下最接近。,,0,,optjmjRmhjRmxxxs（2-1）根据滤波器的形式，维纳滤波器可以分为三种情况：非因果IIR型，因果IIR型，FIR型，对于实时性有要求的情况下用后两种形式。图2.1维纳滤波器对于FIR型维纳滤波器，维纳-霍夫方程的形式为1,,1,010optNjmjRmhjRNmxxxs（2-2）或者写成矩阵形式nwnsnxnhnsnyˆxsoptxxRHR（2-3）其中ToptNhhh110H（2-4）TxsxsxsxsNRRR110R（2-5）021201110xxxxxxxxxxxxxxxxxxxxRNRNRNRRRNRRRR（2-6）这样，如果信号和噪声的二阶统计特性已知，则易求解xsxxRRH1（2-7）维纳滤波的均方误差是1020NmxsoptssmRmhRneE（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滤波器合适的阶数%判据是当阶数增加而均方误差没有明显下降时，则认为阶数足够whileabs(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的全零矩阵作为仿真信号的自相关函数forj=1:Nforn=1:NR_xx(j,n)=Rxx(abs(j-n)+1);endendh=inv(R_xx)*Rxs';%求解互相关函数转置的逆矩阵%xsxxRRH1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);%仿真信号，可以自己生成，任意形式%loadecgdata;%实际心电信号%s=ecgdata(1:M)';%loadeedata;%实际脑电信号%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.调节输入信号为心电信号，心电信号的获取是从专门的心电数据网站提取的，调节信号源以及心电信号的输出方式显示为数据txt的方式，将数据复制存档成txt格式文件，然后复制txt文件到Excel文件中，变成3列，运用列变化提取一列就成为一维的心电信号文件最后Excel截图显示为：图3.1网站心电数据显示图图3.2存档为txt格式的文件图3.3存为Excel截图图3.4存为txt一维心电3.1.4.脑电数据的提取方法跟心电数据的基本一样来源网站信号源的选取为EEG信号，其他部分与心电相似。图3.5网站脑电信号来源注：心电信号和脑电信号都是以txt的格式存取的，其是通过load函数与主程序的输入函数相连接的。具体格式为q=load('C:\Users\tang\Desktop\eeg.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