文章目錄

• • 維納（Wiener）濾波
  • 模型結(jié)構(gòu)
  • 使用條件
  • 原理公式推導(dǎo)
  • 仿真分析——Matlab代碼
  • • 一、參考信號(hào)$d(n)$為原始信號(hào)$s(n)$
    • 二、參考信號(hào)$d(n)$為加性高斯白噪聲$v(n)$
    • 三、參考信號(hào)$d(n)$為輸入信號(hào)自身$x(n)$
  • 總結(jié)

維納（Wiener）濾波

Wiener濾波的核心目標(biāo)就是使均方誤差最小，從而可以推導(dǎo)出維納-霍夫方程。

在信號(hào)處理中，濾波器可以分為FIR和IIR兩種，在這里主要介紹維納FIR濾波器，如下圖1所示。換句話說，就是要想方設(shè)法求出最優(yōu)濾波器的系數(shù)。

圖1 Wiener濾波的橫向?yàn)V波器

模型結(jié)構(gòu)

首先，先來看一下wiener濾波器的一般結(jié)構(gòu)，如下圖2所示。

圖2 wiener濾波的原理框圖

其中$s(n)$是輸入的原始信號(hào)，經(jīng)過噪聲信道$v(n)$后，變成了$x(n)$，濾波器后輸出得到$s^{\prime }(n)$，期望響應(yīng)是$d(n)$ ，那誤差$e(n)=d(n)?s^{\prime }(n)$ ，濾波器的系數(shù)為$h(n)$。

對(duì)于信號(hào)$s(n)$和噪聲$v(n)$的混合體$x(n)=s(n)+v(n)$，按照均方誤差最小的準(zhǔn)則，從$x(t)$中分離出信號(hào)$s(t)$的理論，稱為維納濾波理論。

這里要著重說明的一點(diǎn)是，期望響應(yīng)$d(n)$一般上來說是對(duì)原始信號(hào)的$s(n)$估計(jì)，后面會(huì)講到別的情況。如果你想要對(duì)一個(gè)未知的信號(hào)進(jìn)行濾波，用wiener濾波的方法是不行的。因?yàn)?#xff0c;在設(shè)計(jì)濾波器系數(shù)的時(shí)候，必須用到期望信號(hào)$d(n)$。所以，必須要知道$d(n)$，或者$d(n)$的一些其他的特征。

使用條件

1、輸入信號(hào)是寬平穩(wěn)信號(hào)。那么寬平穩(wěn)是什么呢？通俗的來說就是與時(shí)間無關(guān)的信號(hào)，或者與時(shí)間的起點(diǎn)無關(guān)，只與時(shí)間間隔有關(guān)。

2、必須知道期望信號(hào)，或者它的一些信號(hào)特征才行（具體看下面的公式推導(dǎo)部分）。

原理公式推導(dǎo)

為了簡(jiǎn)便運(yùn)算，假設(shè)所使用的信號(hào)是實(shí)信號(hào)。

輸出信號(hào)：$s^{\prime }(n)=h(n)?x(n)={\sum }_{k=?\infty }^{+\infty }h(n)x(n?k)$

誤差：$e(n)=d(n)?s^{\prime }(n)==d(n)?{\sum }_{k=?\infty }^{+\infty }h(n)x(n?k)$

均方誤差：$J(h)=E[e^2(n)]$

為了使均方誤差最小，需要對(duì)$J$進(jìn)行求導(dǎo)，并讓導(dǎo)數(shù)為0，即可得到最佳濾波器的系數(shù)了。

$${?}_{h}J=?2E[x(n?k)e(n)]$$

所以，
$$E[x(n?k)e(n)]=0(1)$$

$$E\left[x(n?k)\left(d(n)?\sum _{i=?\infty }^{\infty }h(i)x(n?i)\right)\right]=0$$

$$\sum _{i=?\infty }^{\infty }h(i)E\left[x(n?k)x(n?i)\right]=E\left[x(n?k)d(n)\right]$$

$$\sum _{i=?\infty }^{\infty }h(i)r_x(i?k)=r_{xd}(?k)$$

針對(duì)$M$階FIR濾波器，維納-霍夫方程（Wiener-Hopf）為${\sum }_{i=0}^{M?1}h(i)r_x(i?k)=r_{xd}(?k)$，那么，寫成矩陣的形式就是
$$\begin{array}{l}Rh=r_{xd}\\ h=R^{?1}{r}_{xd}\end{array}$$
其中$R$是信號(hào)$x(n)$的自相關(guān)矩陣，$r_{xd}$是信號(hào)$x(n)和d(n)$的互相關(guān)矩陣。

仿真分析——Matlab代碼

由于Wiener濾波器的參數(shù)求解過程中必須要用到參考信號(hào)，所以這里仿真采用的信號(hào)Signal為
$s=A?sin\left(2\pi f_{1}t\right)+B?sin\left(2\pi f_{2}t\right)$ ，
即為兩個(gè)疊加的正弦信號(hào)。其中$A=10,B=15,f_1=1000,f_2=2000$ 。
根據(jù)采樣定理，這里假定采樣頻率$f_s=100000$，采樣間隔$T=1/f_s$，則$s_a(t)|{}_{t=nT}=s_a(nT)(0\le n\le 999)$。
對(duì)于不同的$n$值，$s(n)$是一個(gè)有序的數(shù)字序列:$s(n)=\left\{s_a(0),s_a(T),s_a(2T),?\right\}$。信號(hào)傳輸過程中，信道中的噪聲為加性高斯白噪聲。原始信噪比定義為SNR=3。

上面提到期望信號(hào)$d(n)$是必不可少的，所以，對(duì)期望信號(hào)的設(shè)定也會(huì)決定結(jié)果的好壞。所以，$d(n)$也可以稱作參考信號(hào)。有以下幾種不同的情況：

• 參考信號(hào)$d(n)$為原始信號(hào)$s(n)$
• 參考信號(hào)$d(n)$為加性高斯白噪聲$v(n)$
• 參考信號(hào)$d(n)$為輸入信號(hào)自身$x(n)$

下面對(duì)三種不同的情形進(jìn)行仿真與對(duì)比分析, 其中濾波器系數(shù)長(zhǎng)度N均為100。
（全部三種完整的MATLAB代碼整合版在我的資源“維納（Wiener）濾波及Matlab代碼”中）

一、參考信號(hào)$d(n)$為原始信號(hào)$s(n)$

%% wiener filter for different d(n) %% StarryHuang 2021.1.9 close all; clear; clc; %% 信號(hào)產(chǎn)生，對(duì)原始信號(hào)進(jìn)行采樣 A=10; % 信號(hào)的幅值 B=15; % 信號(hào)的幅值 f1=1000; % 信號(hào)1的頻率 f2=2000; % 信號(hào)2的頻率 fs=10^5; % 采樣頻率 t=0:999; % 采樣點(diǎn)t = [0,999],長(zhǎng)度1000 M = length(t); % 信號(hào)長(zhǎng)度 s=A*sin(2*pi*f1*t/fs) + B*sin(2*pi*f2*t/fs); % 采樣正弦波信號(hào)為Signal SNR = 3; % 初始信噪比 x=awgn(s,SNR,'measured'); %觀測(cè)信號(hào) x=s+v.，給正弦波信號(hào)加入信噪比為-3dB的高斯白噪聲 v=x - s; % 高斯白噪聲，誤差信號(hào)，每次運(yùn)行都不一樣 %% 第一種情況——期望信號(hào)d(n)為感興趣的原信號(hào)Signal d = s; %% 第二種情況——期望信號(hào)d(n)為Noise v % d = v; %% 第三種情況—— 期望信號(hào)d(n)為x(n-1) % d = [x(1),x(1:end-1),]; % d(n)=x(n-1)%% 維納濾波 N = floor(length(x)*0.1); % 濾波器的階數(shù)，向下取整 % N=100; % 維納濾波器階數(shù) Rxx=xcorr(x,N-1,'biased'); % 自相關(guān)函數(shù)1*(2N-1)維度，返回一個(gè)延遲范圍在[-N，N]的互相關(guān)函數(shù)序列,對(duì)稱的 % 變成矩陣 N*N維度 for i=1:Nfor j=1:NmRxx(i,j)=Rxx(N-i+j); % N*N維度end end%產(chǎn)生維納濾波中x 方向上觀測(cè)信號(hào)與期望信號(hào)d的互相關(guān)矩陣 Rxd=xcorr(x,d,N-1,'biased'); % 互相關(guān)函數(shù)1*(2N-1)維度% 變成矩陣1*N維度 for i=1:NmRxd(i)=Rxd(N-1+i); % 1*N維度 endh = inv(mRxx)*mRxd'; % 由wiener-Hopf方程得到濾波器最優(yōu)解, h是N*1維度%% 檢驗(yàn)wiener濾波效果 y = conv(x,h); % 濾波后的輸出,長(zhǎng)度為M+N-1，要截取前M個(gè)。 y = y(1:M); % yy = filter(h,1,x); % 用卷積或者直接用filter都可以 Py=sum(power(y,2))/M; %濾波后信號(hào)y的功率 e = d - y; % 輸出減去期望等于濾波誤差 J = sum(power(e,2))/M % 濾波后噪聲功率 SNR_after = 10*log10((Py-J)/J) % 濾波后信噪比 db單位 imv = 10*log10((Py-J)/J/power(10,SNR/10)) % 濾波后較濾波前信噪比提高了imv dB。%% 畫圖 % 原始信號(hào)x，噪聲v，觀測(cè)波形d figure(1), subplot(311) plot(t,s) % 輸入函數(shù) title(' Signal原信號(hào)')subplot(312) plot(t,v) % 輸入函數(shù) title(' Noise噪聲')subplot(313) plot(t,x) % 輸入函數(shù) title(' X(n)觀測(cè)波形')%% d = s % 期望和濾波后的信號(hào)對(duì)比 figure(2), subplot(211) plot(t, d, 'r:', t, y, 'b-','LineWidth',1); legend('期望信號(hào)','濾波后結(jié)果'); title('期望信號(hào)與濾波結(jié)果對(duì)比'); xlabel('觀測(cè)點(diǎn)數(shù)');ylabel('信號(hào)幅度'); axis([0 1000 -50 50])subplot(212), plot(t, e); title('輸出誤差'); xlabel('觀測(cè)點(diǎn)數(shù)');ylabel('誤差幅度'); axis([0 1000 -50 50])% 濾波前后對(duì)比 figure(3), subplot(211) plot(t, x); title('維納濾波前'); xlabel('觀測(cè)點(diǎn)數(shù)');ylabel('信號(hào)幅度'); axis([0 1000 -50 50])subplot(212), plot(t, y); title('維納濾波后'); xlabel('觀測(cè)點(diǎn)數(shù)');ylabel('誤差幅度'); axis([0 1000 -50 50])

濾波后信噪比SNR_after = 13.187dB；濾波后較濾波前信噪比提高了10.4dB。

二、參考信號(hào)$d(n)$為加性高斯白噪聲$v(n)$

只要將代碼的開頭d換成v，畫圖函數(shù)修改為%% d=v部分即可。

% %% d = v % figure(2) % plot(t, d, 'r:', t, y, 'b-','LineWidth',1); % legend('期望信號(hào)','濾波后結(jié)果'); title('期望信號(hào)與濾波結(jié)果對(duì)比'); % xlabel('觀測(cè)點(diǎn)數(shù)');ylabel('信號(hào)幅度'); % axis([0 1000 -50 50]) % % figure(3) % plot(t, s, 'r:', t, x - y, 'b-','LineWidth',1); % legend('Signal原始信號(hào)','噪聲抵消后結(jié)果'); title('Signal原始信號(hào)與噪聲抵消后結(jié)果對(duì)比'); % xlabel('觀測(cè)點(diǎn)數(shù)');ylabel('信號(hào)幅度'); % axis([0 1000 -50 50])

濾波后信噪比SNR_after = 10.023dB；濾波后較濾波前信噪比提高了7dB。

三、參考信號(hào)$d(n)$為輸入信號(hào)自身$x(n)$

只要將代碼的開頭d換成x即可。


濾波后信噪比SNR_after = -0.812dB；濾波后較濾波前信噪比提高了-3.812dB。

總結(jié)

參考信號(hào)選為原始信號(hào)時(shí)的濾波效果最好。雖然在第三種方案中，參考信號(hào)選為自身信號(hào)時(shí)的濾波效果不好，第三種Wiener濾波器模型對(duì)于許多應(yīng)用仍然具有很大的實(shí)用價(jià)值, 因?yàn)樵谠S多實(shí)際應(yīng)用中, 既無法提前獲知系統(tǒng)的期望響應(yīng), 也不具備獲得與輸入信號(hào)高相關(guān)噪聲的條件。

總結(jié)

以上是生活随笔為你收集整理的维纳（Wiener）滤波及Matlab代码的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。