[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');