完整的维纳滤波器Matlab源程序

完整的维纳滤波器Matlab源程序 学术 2009-11-20 20:31:58 阅读308 评论0 字号：大中小 完整的维纳滤波器Matlab源程序 clear;clc;%输入信号A=1; %信号的幅值f=1000; %信号的频率fs=105; %采样频率t=(0:999); %采样点Mlag=100; %相关函数长度变量 x=A*cos(2*pi*f*t/fs); %输入正弦波信号xmean=mean(x); %正弦波信号均值xvar=var(x,1); %正弦波信号方差xn=awgn(x,5); %给正弦波信号加入信噪比为20dB的高斯白噪声figure(1)plot(t,xn) %绘制输入信号图像title(输入信号图像) xlabel(x轴单位：t/s,color,b)ylabel(y轴单位：f/HZ,color,b)xnmean=mean(xn) %计算输入信号均值xnms=mean(xn.2) %计算输入信号均方值xnvar=var(xn,1) %计算输入信号方差Rxn=xcorr(xn,Mlag,biased); %计算输入信号自相关函数figure(2)subplot(221)plot(-Mlag:Mlag),Rxn) %绘制自相关函数图像title(输入信号自相关函数图像)f,xi=ksdensity(xn); %计算输入信号的概率密度，f为样本点xi处的概率密度subplot(222)plot(xi,f) %绘制概率密度图像title(输入信号概率密度图像)X=fft(xn); %计算输入信号序列的快速离散傅里叶变换Px=X.*conj(X)/600; %计算信号频谱subplot(223)semilogy(t,Px) %绘制在半对数坐标系下频谱图像title(输入信号在半对数坐标系下频谱图像)xlabel(x轴单位：w/rad,color,b)ylabel(y轴单位：w/HZ,color,b)pxx=periodogram(xn); %计算输入信号的功率谱密度subplot(224)semilogy(pxx) %绘制在半对数坐标系下功率谱密度图像title(输入信号在半对数坐标系下功率谱密度图像)xlabel(x轴单位：w/rad,color,b)ylabel(y轴单位：w/HZ,color,b)%fir滤波wp=0.4*pi; %通带截止频率ws=0.6*pi; %阻带截止频率DB=ws-wp; %过渡带宽度N0=ceil(6.6*pi/DB); M=N0+mod(N0+1,2); %计算fir滤波器阶数wc=(wp+ws)/2/pi; %计算理想低通滤波器通带截止频率（关于归一化）hn=fir1(M,wc); %调用fir1计算FIRDF的h(n)y1n=filter(hn,1,xn); %将输入信号通过fir滤波器figure(3)plot(y1n) %绘制经过fir滤波器后信号图像title(经过fir滤波器后信号图像)xlabel(x轴单位：f/HZ,color,b)ylabel(y轴单位：A/V,color,b)y1nmean=mean(y1n) %计算经过fir滤波器后信号均值y1nms=mean(y1n.2) %计算经过fir滤波器后信号均方值y1nvar=var(y1n,1) %计算经过fir滤波器后信号方差Ry1n=xcorr(y1n,Mlag,biased); %计算经过fir滤波器后信号自相关函数figure(4)subplot(221)plot(-Mlag:Mlag),Ry1n) %绘制自相关函数图像title(经过fir滤波器后信号自相关函数图像)f,y1i=ksdensity(y1n); %计算经过fir滤波器后信号的概率密度，f为样本点xi处的概率密度subplot(222)plot(y1i,f) %绘制概率密度图像title(经过fir滤波器后信号概率密度图像)Y1=fft(y1n); %计算经过fir滤波器后信号序列的快速离散傅里叶变换Py1=Y1.*conj(Y1)/600; %计算信号频谱subplot(223)semilogy(t,Py1) %绘制在半对数坐标系下频谱图像title(经过fir滤波器后信号在半对数坐标系下频谱图像)xlabel(x轴单位：w/rad,color,b)ylabel(y轴单位：w/HZ,color,b)py1n=periodogram(y1n); %计算经过fir滤波器后信号的功率谱密度subplot(224)semilogy(py1n) %绘制在半对数坐标系下功率谱密度图像title(经过fir滤波器后信号在半对数坐标系下功率谱密度图像)xlabel(x轴单位：w/rad,color,b)ylabel(y轴单位：w/HZ,color,b)%维纳滤波N=100; %维纳滤波器长度Rxnx=xcorr(xn,x,Mlag,biased); %产生输入信号与原始信号的互相关函数rxnx=zeros(N,1); rxnx(:)=Rxnx(101:101+N-1);Rxx=zeros(N,N); %产生输入信号自相关矩阵Rxx=diag(Rxn(101)*ones(1,N);for i=2:N c=Rxn(101+i)*ones(1,N+1-i); Rxx=Rxx+diag(c,i-1)+diag(c,-i+1);endRxx;h=zeros(N,1);h=inv(Rxx)*rxnx; %计算维纳滤波器的h(n)yn=filter(h,1,xn); %将输入信号通过维纳滤波器figure(5)plot(yn) %绘制经过维纳滤波器后信号图像title(经过维纳滤波器后信号信号图像)xlabel(x轴单位：f/HZ,color,b)ylabel(y轴单位：A/V,color,b)ynmean=mean(yn) %计算经过维纳滤波器后信号均值ynms=mean(yn.2) %计算经过维纳滤波器后信号均方值ynvar=var(yn,1) %计算经过维纳滤波器后信号方差Ryn=xcorr(yn,Mlag,biased); %计算经过维纳滤波器后信号自相关函数figure(6)subplot(221)plot(-Mlag:Mlag),Ryn) %绘制自相关函数图像title(经过维纳滤波器后信号自相关函数图像)f,yi=ksdensity(yn); %计算经过维纳滤波器后信号的概率密度，f为样本点xi处的概率密度subplot(222)plot(yi,f) %绘制概率密度图像title(经过维纳滤波器后信号概率密度图像)Y=fft(yn); %计算经过维纳滤波器后信号序列的快速离散傅里叶变换Py=Y.*conj(Y)/600; %计算信号频谱subplot(223)semilogy(t,Py) %绘制在半对数坐标系下频谱图像title(经过维纳滤波器后信号在半对数坐标系下频谱图像)xlabel(x轴单位：w/rad,color,b)ylabel(y轴单位：w/HZ,color,b)pyn=periodogram(yn); %计算经过维纳滤波器后信号的功率谱密度subplot(224)semilogy(pyn) %绘制在半对数坐标系下功率谱密度图像title(经过维纳滤波器后信号在半对数坐标系下功率谱密度图像)xlabel(x轴单位：w/rad,color,b)ylabel(y轴单位：w/HZ,color,b)-Matlab维纳滤波器最小均方误差可以为零吗v(n)=sin(2*pi/4*n)s(n)=sin(2*pi/16*n)x(n)=s(n)+v(n)y(n)=s(n)e(n)=y(n)-x(n)*h(n)h=rxx-1*ryx rxx,ryx分别是x的自相关函数和y、x的互相关函数=(rvv+rss)-1*rss计算的r00、r11前8个数值如下：rss=0.0313,0.0289,0.0227,0.0141,0.0047,-0.0039,-0.0105,-0.0144;rvv=0.5,0,-0.25,0,0.5,0,-0.25,0;fir取6阶时，对应最小均方误差为：0.0089fir取8阶时，对应最小均方误差为：0.0074均没有达到最小误差为零的要求，应该可以理解为x(n)经过维纳滤波器h(n)的输出没有得到所期望的y(n)=s(n).但若采用一般的滤波器，由于s(n)、v(n)的频谱分离，应该能够无失真的恢复出信号s(n).若照这样理解的话是不是可以认为一般的滤波器在某些情况下要比用维纳滤波器好一些，还是我理解有错误，欢迎大家讨论。-matlab 中哪个命令可以实现维纳滤波器的功能急用维纳滤波器,我需要它来仿真信号的消噪 本帖最后由 eight 于 2008-3-24 17:18 编辑 搜索更多相关主题的帖子: matlab 维纳滤波器 信号 UID65964 帖子10 精华0 积分1 体能10 点 威望0 点 管理积分0 点 阅读权限10 性别男 来自湖南 在线时间2 小时 注册时间2007-6-5 最后登录2007-6-22 查看详细资料 TOP 花如月 翰林学士 个人空间 发短消息 加为好友 当前离线 沙发 大 中 小 发表于 2007-6-6 14:39 只看该作者 I,map=imread(eight.tif);J1=imnoise(I,gaussian,0,0.02); % 受高斯噪声干扰，结果如图3-17(b)所示K noise=wiener2(J1, 5 5);%自适应维纳滤滤subplot(131),imshow(I);title(原图l)subplot(132),imshow(J1);title(噪声干扰图)subplot(133),imshow(K);title(滤波后图像)自适应维纳滤滤波示例，仅供参考。更具体的用法参wiener2帮助文档.设计方法很多：冲击响应不变法，双线性Z变换法等具体步骤：（1） 将数字滤波器性能指标变换成模拟滤波器指标 （2） 归一化到低通滤波器的设计指标 （3） 设计模拟低通滤波器 （4） 转换成所需的模拟滤波器转移函数 （5） 将所得的转移函数数字化得到需要的数字滤波器说明：此方法适合与带通、带阻、高通滤波器的设计以上每个步骤都有相应的实现函数，不再一一列出函数butter（巴特沃数字斯滤波器）函数chey1（切比雪夫1型数字斯滤波器）可以实现上边的5个步骤；更详细的请参看文献数字信号处理 清华大学出版社 胡广书=这里给出一个陷波器的例子:% 50Hz notch filter% sample frequency=400%clear all;clc;b=1 -sqrt(2) 1;a=1 -sqrt(2)*0.999 0.999;db, mag, pha, grd, w=freqz_m(b, a); subplot(221); plot(w*200/pi, db); title( Magnitude Response );xlabel(frequency in Hz); ylabel(dB); axis(0, 100, -200, 5);set(gca, XTickMode, manual, XTick, 0, 50, 100); set(gca, YTickmode, manual, YTick, -200, -100, 0); gridtitle(Notch filter response);t0=1:8000;t=1:256;t1=1:100;t2=1:128;x=sin(2*pi*50*t0/400)+0.5*sin(2*pi*100*t0/400);x1=x(t);y=filter(b,a,x1);subplot(222); plot(x1);title(Original waveform);X=fft(x1);subplot(223); plot(t2*400/256,abs(X(t2);xlabel(frequency in Hz); ylabel(|H|); axis(0, 200, 0, 150);title(Spectrum for original);set(gca, XTickMode, manual, XTick, 0, 50, 100, 150); set(gca, YTickmode, manual, YTick, 50, 100); gri