To

learn, to share, to debate, then comes

progress.

1.算法背景：

信號濾波的實(shí)質(zhì)為從觀測信號中提取有效信號，隨著數(shù)學(xué)理論的發(fā)展與實(shí)際應(yīng)用的需求，基于不同原理的濾波方法被不斷地提出來，雖然依據(jù)的準(zhǔn)則，推導(dǎo)的過程各有差異，但最終的目的均是減小信號估計(jì)的誤差，使濾波系統(tǒng)的輸出信號盡可能地接近實(shí)際信號。

Wiener濾波是第二次世界大戰(zhàn)中，為了解決火力控制系統(tǒng)精確跟蹤問題，Wiener相繼提出了平穩(wěn)隨機(jī)過程的最優(yōu)線性濾波理論，首次將數(shù)理統(tǒng)計(jì)知識和線性系統(tǒng)理論聯(lián)系起來，形成了對隨機(jī)信號作平滑，濾波和預(yù)測的最新估計(jì)理論。在此后的發(fā)展中，Wiener濾波被應(yīng)用于更多的領(lǐng)域，并沿用至今。

2.算法原理：

（1）有限長濾波器

對于一列輸入信號x，一般的無限長線性濾波器輸出為：

y(n)=

Σh(m)x(n-m) m=0…∞

實(shí)際中，濾波器的長度，即階數(shù)是有限長的，設(shè)為M，則有：

y(n)=

Σh(m)x(n-m) m=0…M

即濾波器的當(dāng)前時刻輸出為前M個時刻的值經(jīng)過加權(quán)之后得到的。

為便于書寫與理解，上式可以寫為矩陣形式：

y(n)=H(m)*X(n)

如果期望信號d已知，則可以計(jì)算輸出與期望信號之間的誤差：

e(n)=d(n)-y(n)= d(n)- H(m)*X(n) m=0…M

Wiener濾波的目標(biāo)就是，如何確定一個長為M的系數(shù)序列H，使得上述誤差值最小。

（2）最小均方誤差濾波

根據(jù)目標(biāo)函數(shù)的不同，又可以將濾波算法細(xì)分為不同的類別，一般來說有最小均方誤差，最小二乘誤差等等，這里只討論最小均方誤差。

令目標(biāo)函數(shù)為：

Min

E[e(n)^2]= E[(d(n)- H(m)*X(n))^2]

當(dāng)濾波器的系數(shù)最優(yōu)時，目標(biāo)函數(shù)對系數(shù)的倒數(shù)應(yīng)該為0，即：

dE[e(n)^2]/dH=0

2 E[ (d(n)-

H(m)*X(n))]* X(n)=0

E[(d(n)

X(n))- H(m)E[X(n)X(n)]=0

根據(jù)隨機(jī)過程的知識，上式可以表達(dá)為：

Rxd-H*Rxx=0

其中Rxd與Rxx分別為輸入信號與期望信號的相關(guān)矩陣與輸入信號的自相關(guān)矩陣。

從而有：

H=Rxx-1*Rxd

至此，便得到了Wiener濾波的基本原理與公式推導(dǎo)。

3.算法應(yīng)用與實(shí)現(xiàn)

理解了算法的原理之后，下邊舉一個小的例子來考察如何應(yīng)用Wienar濾波處理實(shí)際問題。

問題背景：一個點(diǎn)目標(biāo)在x，y平面上繞單位圓做圓周運(yùn)動，由于外界干擾，其運(yùn)動軌跡發(fā)生了偏移。其中，x方向的干擾為均值為0，方差為0.05的高斯噪聲；y方向干擾為均值為0，方差為0.06的高斯噪聲。

問題分析與思路：

將物體的運(yùn)動軌跡分解為X方向和Y方向，并假設(shè)兩個方向上運(yùn)動相互獨(dú)立。分別將運(yùn)動軌跡離散為一系列點(diǎn)，作為濾波器的輸入，分別在兩個方向上進(jìn)行濾波，最終再合成運(yùn)動軌跡。

程序設(shè)計(jì)思路：

生成期望信號-添加噪聲-計(jì)算相關(guān)矩陣-求解最佳濾波器系數(shù)-濾波運(yùn)算-輸出信號-合成軌跡

4.結(jié)果與分析

5.源代碼

%***********************************************

%該程序使用Wiener濾波方法對圓周運(yùn)動軌跡進(jìn)行控制

%信號模型：d=s+no

觀測信號=期望信號+噪聲信號

%進(jìn)行一次Wiener濾波，得到最佳濾波器系數(shù)

17.4 by

Howie

clear

close all

N=500;

theta=linspace(0,2*pi,N); %極坐標(biāo)參數(shù)

s_x=cos(theta); %x，y方向上的期望信號

s_y=sin(theta);

no_x=normrnd(0,sqrt(0.05),1,N); %高斯白噪聲

no_y=normrnd(0,sqrt(0.06),1,N);

d_x=s_x+no_x; %觀測信號

d_y=s_y+no_y;

M=500;%M為濾波器的階數(shù)

%% 對x方向上數(shù)據(jù)進(jìn)行濾波

rxx=xcorr(d_x);

Rxx=zeros(N);

% temp=toeplitz(rxx);

for

i=1:N %觀測信號的相關(guān)矩陣

for

j=1:N

Rxx(i,j)=rxx(N+i-j);

end

end

rxd=xcorr(s_x,d_x); %觀測信號與期望信號的相關(guān)矩陣

Rxd=rxd(N:N+M-1); %向量而非矩陣

hopt_x=Rxx\Rxd';

% de_x=conv(hopt_x,d_x);

de_x=zeros(1,N);

for n=1:N

for

i=1:n-1

de_x(n)=de_x(n)+hopt_x(i)*d_x(n-i);

end

end

de_x(1:2)=d_x(1:2);

ems_x=sum(d_x.^2)-Rxd*hopt_x;

e_x=de_x-s_x;

% de_x(N-1:N)=d_x(N-1:N);

%% 對y方向上數(shù)據(jù)進(jìn)行濾波 處理思路同x方向

ryy=xcorr(d_y);

Ryy=zeros(N);

for i=1:N

for

j=1:N

Ryy(i,j)=ryy(N+i-j);

end

end

% temp=toeplitz(ryy);

% Ryy=temp(1:M,N:N+M-1);

ryd=xcorr(s_y,d_y);

% temp=toeplitz(ryd);

%

Ryd=temp(1:N,N:length(temp));

Ryd=ryd(N:N+M-1);

hopt_y=Ryy\Ryd';

% de_y=conv(hopt_y,d_y);

de_y=zeros(1,N);

for n=1:N

for

i=1:n-1

de_y(n)=de_y(n)+hopt_y(i)*d_y(n-i);

end

end

de_y(1:2)=d_y(1:2);

ems_y=sum(d_y.^2)-Ryd*hopt_y;

e_y=de_y-s_y;

% de_y(N-1:N)=d_y(N-1:N);

%% plot

figure

plot(s_x,s_y,'r','linewidth',2)

hold on

plot(d_x,d_y,'b')

hold on

plot(de_x,de_y,'k-')

title('維納濾波預(yù)測軌跡')

legend('期望軌跡','觀測軌跡','濾波軌跡')

%% %% x方向上繪圖

figure

suptitle('X方向上維納濾波效果')

subplot(321)

plot(s_x)

title('期望信號')

subplot(322)

plot(no_x)

title('噪聲信號')

subplot(323)

plot(d_x)

title('觀測信號')

subplot(324)

plot(de_x)

title('濾波后信號')

subplot(325)

plot(ems_x,'o')

title('最小均方誤差')

subplot(326)

plot(e_x)

title('絕對誤差')

%% y方向上繪圖

figure

suptitle('Y方向上維納濾波效果')

subplot(321)

plot(s_y)

title('期望信號')

subplot(322)

plot(no_y)

title('噪聲信號')

subplot(323)

plot(d_y)

title('觀測信號')

subplot(324)

plot(de_y)

title('濾波后信號')

subplot(325)

plot(ems_y,'o')

title('最小均方誤差')

subplot(326)

plot(e_y)

title('絕對誤差')

總結(jié)

以上是生活随笔為你收集整理的维纳滤波python 函数_Wiener维纳滤波基本原理及其算法实现的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

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