[Matlab]维纳滤波器设计

[Matlab]维纳滤波器设计

​ 维纳滤波（wiener filtering) 一种基于最小均方误差准则、对平稳过程的最优估计器。这种滤波器的输出与期望输出之间的均方误差为最小，因此，它是一个最佳滤波系统。它可用于提取被平稳噪声污染的信号。

​ 从连续的（或离散的）输入数据中滤除噪声和干扰以提取有用信息的过程称为滤波，这是信号处理中经常采用的主要方法之一，具有十分重要的应用价值，而相应的装置称为滤波器。根据滤波器的输出是否为输入的线性函数，可将它分为线性滤波器和非线性滤波器两种。维纳滤波器是一种线性滤波器。

基本概念

​ 从噪声中提取信号波形的各种估计方法中，维纳（Wiener）滤波是一种最基本的方法，适用于需要从噪声中分离出的有用信号是整个信号（波形），而不只是它的几个参量。

设维纳滤波器的输入为含噪声的随机信号。期望输出与实际输出之间的差值为误差，对该误差求均方，即为均方误差。因此均方误差越小，噪声滤除效果就越好。为使均方误差最小，关键在于求冲激响应。如果能够满足维纳－霍夫方程 [3] ，就可使维纳滤波器达到最佳。根据维纳－霍夫方程，最佳维纳滤波器的冲激响应，完全由输入自相关函数以及输入与期望输出的互相关函数所决定。

维纳滤波器优缺点

维纳滤波器的优点是适应面较广，无论平稳随机过程是连续的还是离散的，是标量的还是向量的，都可应用。对某些问题，还可求出滤波器传递函数的显式解，并进而采用由简单的物理元件组成的网络构成维纳滤波器。维纳滤波器的缺点是，要求得到半无限时间区间内的全部观察数据的条件很难满足，同时它也不能用于噪声为非平稳的随机过程的情况，对于向量情况应用也不方便。因此，维纳滤波在实际问题中应用不多。

实现维纳滤波的要求是：

1.输入过程是广义平稳的

2.输入过程的统计特性是已知的。根据其他最佳准则的滤波器亦有同样要求

然而，由于输入过程取决于外界的信号、干扰环境，这种环境的统计特性常常是未知的、变化的，因而难以满足上述两个要求。这就促使人们研究自适应滤波器。

维纳滤波器原理分析：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;clear all; close all;  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%输入信号
A=1;                                                      %信号的幅值
f=1000;                                                 %信号的频率
fs=10^5;                                                %采样频率
t=(0:999);                                              %采样点
Mlag=100;                                             %相关函数长度变量   
x=A*cos(2*pi*f*t/fs);                                %输入正弦波信号
xmean=mean(x);                                    %正弦波信号均值
xvar=var(x,1);                                         %正弦波信号方差
noise=wgn(1,1000,2);%产生1行1000列的矩阵，强度为2dbw
xn=x+noise;                                         %给正弦波信号加入信噪比为20dB的高斯白噪声
plot(t,xn)    
xlabel('x轴单位：t/s','color','b')
ylabel('y轴单位：A/V','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')
 
%维纳滤波
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);
end
Rxx;
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')
subplot(4,1,1),plot(noise); title('噪声信号')
subplot(4,1,2),plot(x); title('正弦信号')
subplot(4,1,3),plot(xn); title('加噪信号')
subplot(4,1,4),plot(yn); title('维纳信号')

维纳滤波器函数设计：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%维纳滤波器函数设计
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y =wienerfilter(x,Rxx,Rxd,N) 
%进行维纳滤波 
%x是输入信号，Rxx是输入信号的自相关向量 
%Rxd是输入信号和理想信号的的互相关向量，N是维纳滤波器的长度 
%输出y是输入信号通过维纳滤波器进行维纳滤波后的输出 
h=yulewalker(Rxx,Rxd,N);								%求解维纳滤波器系数 
t=conv(x,h);											%进行滤波 
Lh=length(h);											%得到滤波器的长度 
Lx=length(x);											%得到输入信号的长度 
y=t(double(uint16(Lh/2)):Lx+double(uint16(Lh/2))-1);%输出序列y的长度和输入序列x的长度相同
%以下是维纳滤波器系数的求解 
function h=yulewalker(A,B,M)    
%求解Yule-Walker方程 
%A是接收信号的自相关向量为 Rxx(0),Rxx(1),......,Rxx(M-1) 
%B是接收信号和没有噪声干扰信号的互相关向量为 Rxd(0),Rxd(1),......,Rxd(M-1) 
%M是滤波器的长度 
%h保存滤波器的系数 
T1=zeros(1,M);%T1存放中间方程的解向量 
T2=zeros(1,M);%T2存放中间方程的解向量 
T1(1)=B(1)/A(1); 
T2(1)=A(2)/A(1); 
X=zeros(1,M); 
for i=2:M-1 
temp1=0; 
temp2=0; 
    for j=1:i-1 
        temp1=temp1+A(i-j+1)*T1(j); 
        temp2=temp2+A(i-j+1)*T2(j); 
    end 
    X(i)=(B(i)-temp1)/(A(1)-temp2); 
    for j=1:i-1 
        X(j)=T1(j)-X(i)*T2(j); 
    end 
    for j=1:i 
        T1(j)=X(j); 
    end 
temp1=0; 
temp2=0; 
    for j=1:i-1 
        temp1=temp1+A(j+1)*T2(j); 
        temp2=temp2+A(j+1)*T2(i-j); 
    end 
    X(1)=(A(i+1)-temp1)/(A(1)-temp2); 
    for j=2:i 
        X(j)=T2(j-1)-X(1)*T2(i-j+1); 
    end 
    for j=1:i 
        T2(j)=X(j); 
    end 
end 
temp1=0; 
temp2=0; 
for j=1:M-1 
	temp1=temp1+A(M-j+1)*T1(j); 
	temp2=temp2+A(M-j+1)*T2(j); 
end 
X(M)=(B(M)-temp1)/(A(1)-temp2); 
for j=1:M-1 
	X(j)=T1(j)-X(M)*T2(j); 
end 
h=X;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%维纳滤波器案例测试
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;clear all; close all;  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load handel						%加载语音信号
d=y; d=d*8;						%增强语音信号强度
d=d';
[m,n]=size(d);
T = 1/Fs; % 采样时间
t = (1:n)*T;% 时间
subplot(3,2,1);
plot(t,d);				
title('原始语音信号');
xlabel('时间/t');
ylabel('幅值/dB');
fq=fft(d,8192);						%进行傅立叶变换得到语音信号频频
subplot(3,2,2);
f=Fs*(0:4095)/8192;
plot(f,abs(fq(1:4096)));				%画出频谱图
title('原始语音信号的频域图形');
xlabel('频率 f');
ylabel('FFT');
x_noise=randn(1,n);				%（0，1）分布的高斯白噪声
x=d+x_noise;						%加入噪声后的语音信号
subplot(3,2,3);
plot(t,x);				
title('加入噪声后');
xlabel('时间/t');
ylabel('幅值/dB');
fq=fft(x,8192);						%对加入噪声后的信号进行傅立叶变换，看其频谱变化
subplot(3,2,4);
plot(f,abs(fq(1:4096)));				%画出加入噪声后信号的频谱图
title('加入噪声后语音信号的频域图形');
xlabel('频率 f');
ylabel('FFT');
%维纳滤波
yyhxcorr=xcorr(x(1:4096));			%求取信号的信号的自相关函数
size(yyhxcorr); 
A=yyhxcorr(4096:4595);
yyhdcorr=xcorr(d(1:4096),x(1:4096));			%求取信号和理想信号的互相关函数
size(yyhdcorr);
B=yyhdcorr(4096:4595);
M=500;
yyhresult=wienerfilter(x,A,B,M);				%进行维纳滤波
yyhresult=yyhresult(300:8192+299);
subplot(3,2,5);
t = (1:8192)*T;% 时间
plot(t,yyhresult);				%画出频谱图
title('进行维纳滤波');
xlabel('时间/t');
ylabel('幅值/dB');
fq=fft(yyhresult);							%对维纳滤波的结果进行傅立叶变换，看其频谱变化
subplot(3,2,6); 
f=Fs*(0:4095)/8192;
plot(f,abs(fq(1:4096)));						%画出维纳滤波后信号的频谱图
title('经过维纳滤波后语音信号的频域图形');
xlabel('频率 f');
ylabel('FFT');