维纳-霍夫方程 实验报告.doc

生物医学信号处理评分大理大学实验报告 课程名称 生物医学信号处理 实验名称 维纳-霍夫方程 专业班级 姓 名 羽卒兰cl 学 号 实验日期 2016年5月20日 实验地点 20152016学年度第 3 学期一、 实验目的学习求解维纳-霍夫方程，寻找最小均方误差意义下的最优滤波器。二、实验环境 1、硬件配置：Intel（R） Core(TM) i5-4210U CPU 1.7GHz 1.7GHz 安装内存（RAM）：4.00GB 系统类型：64位操作系统 2、软件环境：MATLAB R2013b软件三、实验内容(包括本实验要完成的实验问题及需要的相关知识简单概述)已知信号的自相关函数和噪声的能量，编写程序求解维纳-霍夫方程，寻找最优滤波器。编写程序仿真信号，噪声和观察波形，然后把观察信号通过滤波器得到的信号估计与原始信号比较，观察是否达到了去噪的目的。选择不同信号（仿真信号，实际采集的心电，脑电信号），人工添加噪声，调整噪声的相对强度，观察滤波效果。四、实验结果与分析 （包括实验原理、数据的准备、运行过程分析、源程序（代码）、图形图象界面等）实验原理根据正交原理可以推导出维纳-霍夫方程，满足该方程的滤波器输出信号的估计值与信号在最小均方误差意义下最接近。根据滤波器的形式，维纳滤波器可以分为三种情况：非因果IIR型，因果IIR型，FIR型，对于实时性有要求的情况下用后两种形式。图1 维纳滤波器对于FIR型维纳滤波器，维纳-霍夫方程的形式为：或者写成矩阵形式：其中：这样，如果信号和噪声的二阶统计特性已知，则易求解。维纳滤波的均方误差是：源程序：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的取值加1e0 = e1; %把e1赋给e0Rxs = Rss(M:(M+N-1); %信号和噪声互不相关Rxx = Rww(M:(M+N-1)+Rss(M:(M+N-1);%信号和噪声的自相关之和R_xx = zeros(N); %生成N阶的零矩阵for j = 1:N for n = 1:N R_xx(j,n) = Rxx(abs(j-n)+1); %双for为了求出对称矩阵 endendh = inv(R_xx)*Rxs; %求rxs的逆矩阵e1 = Rss(M)-h*Rxs; %转置endN % 显示N的最终值e = e1; %真实值与估计值之间的估计误差%主程序clear; clc; %清空M = input(信号的长度 M = ); %输入n = 1:M; %n的取值范围从1到Ms = exp(-0.002*n).*sin(pi*n/50); %指数衰减正弦,仿真信号，可以自己生成，任意形式% load ecgdata; %导入实际心电信号的数据%load eegdata; %导入脑电信号的数据%load icpdata; %导入颅内压信号的数据%load respdata; %导入个呼吸信号的数据%s = ecgdata(1:M); %取生成心电信号(1:M)的逆矩阵%s = eegdata(1:M); %取生成脑电信号(1:M)的逆矩阵%s = icpdata(1:M); %取生成颅内压信号(1:M)的逆矩阵%s = respdata(1:M); %取生成呼吸信号(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(ss-s).2)%求最小均方误差思考题：观察实验结果，对于几种不同的信号，维纳滤波是否都取得了较好的效果？如果效果不好，试分析原因。1.只改变数据长度，观察其对信号的影响1.1仿真信号，噪声强度w=0.4，阈值取10-4，M分别取1024，2000，3000，4000 信号的长度M=1024 N=103 error=0.0092 信号的长度M=2000 N=114 error=0.0079 信号的长度M=3000 N=200 error=0.0044 信号的长度M=4000 N=204 error=0.0050图1 仿真信号的维纳滤波图（数据长度M分别为1024，2000,3000,4000）答：从图1可以看出，随着仿真信号长度M分别为1024，2000,3000,4000的增大，N值也逐渐增大，最小均方误差error有下降的趋势，因此信号估计图越接近原信号，主要是因为M越大，样本点数越大，原始信号的采集量越大，结果比较准确。1.2心电信号，噪声强度w=0.4，阈值取10-4，M分别取128，256，512，1024 信号的长度M=128 N=46 error=0.0326 信号的长度M=256 N=76 error=0.0189 信号的长度M=512 N=89 error=0.0234 信号的长度M=1024 N=176 error=0.0229 图2 心电信号的维纳滤波图（数据长度M分别为128,256,512,1024） 答：从图2可以看出，随着心电信号长度M分别为128,256,512，1024的增大，N值也逐渐增大，最小均方误差error有下降的趋势，因此信号估计图越接近原信号,但是心电信号有局限性，所以取值的时候要注意，要想取到1024以上，可以进行周期延拓。1.3脑电信号，噪声强度w=0.4，阈值取10-4，M分别取128，256，512，1024 信号的长度M=128 N=28 error=0.1261 信号的长度M=256 N=77 error=0.1497 信号的长度M=512 N=162 error=0.1149 信号的长度M=1024 N=444 error=0.1317 图3 脑电信号的维纳滤波图（数据长度M分别为128,256,512,1024）答：从图3可以看出，随着脑电信号长度M分别为128,256,512，1024的增大，N值逐渐增大，最小均方误差error有上升的趋势，那是因为脑电信号是非平稳随机信号，不适用，脑电信号有局限性，所以取值的时候要注意，要想取到1024以上，可以进行周期延拓。1.4颅内压信号，噪声强度w=0.4，阈值取10-4，M分别取128，256，512，1024 信号的长度M=128 N=81 error=0.1266 信号的长度M=256 N=188 error=0.0744 信号的长度M=512 N=197 error=0.0819 信号的长度M=1024 N=264 error=0.0747 图4 颅内压信号的维纳滤波图（数据长度M分别为128,256,512,1024）答：从图4可以看出，随着颅内压信号长度M分别为128,256,512，1024的增大，N值逐渐增大，最小均方误差error有下降的趋势，颅内压信号有局限性，所以取值的时候要注意，要想取到1024以上，可以进行周期延拓。1.5呼吸信号，噪声强度w=0.4，阈值取10-4，M分别取128，256，512，1024 信号的长度M=128 N=3 error=0.1622 信号的长度M=256 N=54 error=0.1404 信号的长度M=512 N=113 error=0.1581 信号的长度M=1024 N=264 error=0.1513 图5 呼吸信号的维纳滤波图（数据长度M分别为128,256,512,1024）答：从图5可以看出，随着呼吸信号长度M分别为128,256,512，1024的增大，N值逐渐增大，最小均方误差error有下降的趋势，呼吸信号有局限性，所以取值的时候要注意，要想取到1024以上，可以进行周期延拓。2.只改变噪声强度，观察噪声强度对信号的影响2.1仿真信号，阈值取10-6，M=1024，W分别取0.1,0.2,0.4,0.8 信号的长度M=1024 N=37 error=0.0014（W=0.1） 信号的长度M=1024 N=92 error=0.0042（W=0.2） 信号的长度M=1024 N=88 error=0.0118（W=0.4） 信号的长度M=1024 N=183 error=0.0234（W=0.8） 图6 仿真信号的维纳滤波图（噪声强度W分别取0.1,0.2,0.4,0.8）答：从图6可以看出，随着仿真信号的噪声强度0.1,0.2,0.4,0.8的增大，N有逐渐增大的趋势，最小均方误差error越大，维纳滤波效果越差。2.2心电信号，阈值取10-6，M=1024，W分别取0.1,0.2,0.4,0.8 信号的长度M=1024 N=46 error=0.0042（W=0.1） 信号的长度M=1024 N=87 error=0.0119（W=0.2） 信号的长度M=1024 N=120 error=0.0220（W=0.4） 信号的长度M=1024 N=322 error=0.0319（W=0.8） 图7 心电信号的维纳滤波图（噪声强度W分别取0.1,0.2,0.4,0.8）答：从图7可以看出，随着心电信号的噪声强度0.1,0.2,0.4,0.8的增大，N逐渐增大，最小均方误差error越大，维纳滤波效果越差。2.3脑电信号，阈值取10-6，M=1024，W分别取0.1,0.2,0.4,0.8 信号的长度M=1024 N=3 error =0.0100（W=0.1） 信号的长度M=1024 N=13，error =0.0404（W=0.2） 信号的长度M=1024 N=39 error =0.1316（W=0.4） 信号的长度M=1024 N =121 error =0.4509（W=0.8）图8 脑电信号的维纳滤波图（噪声强度W分别取0.1,0.2,0.4,0.8）答：从图8可以看出，随着脑电信号的噪声强度增大，N逐渐增大，最小均方误差error越大，维纳滤波效果越差。2.4颅内压信号，阈值取10-6，M=1024，W分别取0.1,0.2,0.4,0.8 信号的长度M=1024 N=11 error =0.0091（W=0.1） 信号的长度M=1024 N=20，error =0.0301（W=0.2） 信号的长度M=1024 N=59 error =0.0692（W=0.4） 信号的长度M=1024 N =125 error =0.1760（W=0.8）图9 颅内压信号的维纳滤波图（噪声强度W分别取0.1,0.2,0.4,0.8）答：从图9可以看出，随着颅内压信号的噪声强度增大，N逐渐增大，最小均方误差error越大，维纳滤波效果越差。2.5呼吸信号，阈值取10-6，M=1024，W分别取0.1,0.2,0.4,0.8 信号的长度M=1024 N=2 error=0.0095（W=0.1） 信号的长度M=1024 N=47 error=0.0365（W=0.2） 信号的长度M=1024 N=110 error=0.1503（W=0.4） 信号的长度M=1024 N=228 error=0.6752（W=0.8） 图10 呼吸信号的维纳滤波图（噪声强度W分别取0.1,0.2,0.4,0.8）答：从图10可以看出，随着呼吸信号的噪声强度0.1,0.2,0.4,0.8的增大，N逐渐增大，最小均方误差error越大，维纳滤波效果越差。3.只改变阶数，观察阶数对信号的影响3.1仿真信号，W=0.4，M=1024，阈值分别取10-1，10-3，10-6，10-9 信号的长度M=1024 N=58 error=0.0094（10-1） 信号的长度M=1024 N=85 error=0.0091（10-3） 信号的长度M=1024 N=87 error=0.0087（10-6） 信号的长度M=1024 N=128 error=0.0083（10-9） 图11 仿真信号的维纳滤波图（阈值分别取10-1，10-3，10-6，10-9）答：从图11中可以看出，随着仿真信号的阈值10-1，10-3，10-6，10-9的减小，N逐渐增大，最小均方误差error有下降的趋势，维纳滤波效果越好。3.2心电信号，W=0.4，M=1024，阈值分别取10-1，10-3，10-6，10-9 信号的长度M=1024 N=58 error=0.0257（10-1） 信号的长度M=1024 N=102 error=0.0197（10-3） 信号的长度M=1024 N=109 error=0.0211（10-6） 信号的长度M=1024 N=99 error=0.0223（10-9） 图12 心电信号的维纳滤波图（阈值分别取10-1，10-3，10-6，10-9）答：从图12中可以看出，随着心电信号的阈值10-1，10-3，10-6，10-9的减小，N有增大的趋势，最小均方误差error有下降的趋势，维纳滤波效果越好。3.3脑电信号，W=0.4，M=1024，阈值分别取10-1，10-3，10-6，10-9 信号的长度M=1024 N=27 error=0.1386（10-1） 信号的长度M=1024 N=114 error=0.1311（10-3） 信号的长度M=1024 N=34 error=0.1349（10-6） 信号的长度M=1024 N=80 error=0.1380（10-9） 图13 脑电信号的维纳滤波图（阈值分别取10-1，10-3，10-6，10-9）答：从图13中可以看出，随着脑电信号的阈值10-1，10-3，10-6，10-9的减小，N有增大的趋势，最小均方误差error有下降的趋势，维纳滤波效果越好。3.4颅内压信号，W=0.4，M=1024，阈值分别取10-1，10-3，10-6，10-9 信号的长度M=1024 N=25 error=0.0861（10-1） 信号的长度M=1024 N=56 error=0.0783（10-3） 信号的长度M=1024 N=113 error=0.0802（10-6） 信号的长度M=1024 N=228 error=0.0752（10-9） 图14 颅内压信号的维纳滤波图（阈值分别取10-1，10-3，10-6，10-9）答：从图14中可以看出，随着颅内压信号的阈值10-1，10-3，10-6，10-9的减小，N逐渐增大，最小均方误差error越小，维纳滤波效果越好。3.5呼吸信号，W=0.4，M=1024，阈值分别取10-1，10-3，10-6，10-9 信号的长度M=1024 N=3 error=0.1612（10-1） 信号的长度M=1024 N=8 error=0.1633（10-3） 信号的长度M=1024 N=12 error=0.1580（10-6） 信号的长度M=1024 N=3 error=0.1534（10-9） 图15 呼吸信号的维纳滤波图（阈值分别取10-1，10-3，10-6，10-9）答：从图15中可以看出，随着呼吸信号的阈值10-1，10-3，10-6，10-9的减小，N有增大的趋势，最小均方误差error有下降的趋势，维纳滤波效果越好。4.只改变信号的类型，观察信号类型对信号的影响（五种信号）4.1 W=0.4，M=1024，阈值=10-6，信号为：仿真信号 图16 信号的长度M=1024 N=57 error=0.0076（仿真信号维纳滤波图）答：观察图16，比较仿真的信号与估计信号，可以看出噪声对信号有很大影响，虽然可以看出有用信号的轮廓，但没能消除噪声的影响，滤波效果还可以。4.2 W=0.4，M=1024，阈值=10-6，信号为：心电信号 图17 信号的长度M=1024 N=66 error=0.0212（心电信号维纳滤波图）答：观察图17，比较心电的信号与估计信号，可以看出噪声对信号有很大影响，估计信号很难辨别，估计误差图中产生了与原信号相反的波峰，滤波效果差。4.3 W=0.4，M=1024，阈值=10-6，信号为：脑电信号 图18 信号的长度M=1024 N=40 error=0.1468（脑电信号维纳滤波图）答：观察图18，比较脑电信号与估计信号，可以看出估计信号和原信号噪声相似，脑电信号对噪声很敏感，最小均方误差大。因为维纳滤波处理的是平稳随机信号，对脑电信号这种非平稳随机信号的处理效果不好。4.4 W=0.4，M=1024，阈值=10-6，信号为：颅内压信号 图19 信号的长度M=1024 N=45 error=0.0797（颅内压信号维纳滤波图）答：观察图19，比较颅内压信号与