《MATLAB輔助現(xiàn)代工程數(shù)字信號(hào)處理》課件第8章

第8章非平穩(wěn)信號(hào)分析與處理8.1短時(shí)傅里葉變換

8.2維格納時(shí)頻分布

8.3小波變換

8.4小結(jié)

8.1短時(shí)傅里葉變換

8.1.1時(shí)域窗法

對(duì)于時(shí)變信號(hào)，需了解它的局部傅里葉變換，最簡(jiǎn)單的方法就是對(duì)信號(hào)乘上一個(gè)滑動(dòng)窗來研究它的傅里葉變換。

設(shè)離散信號(hào)為x(n)，在n0時(shí)刻局部的傅里葉變換為(8.1)當(dāng)n0變化時(shí)，可以得到信號(hào)x(n)的短時(shí)傅里葉變換為(8.2)式中，w(n)為窗函數(shù)，其長(zhǎng)度為Nw，通?？蛇x為中心對(duì)稱。定義n時(shí)刻的短時(shí)段信號(hào)為

xn(m)=x(m)w(n－m)

(8.3)

因此，時(shí)域窗法就是用一個(gè)中心對(duì)稱的滑動(dòng)窗函數(shù)截取觀測(cè)信號(hào)，對(duì)不同時(shí)刻的短時(shí)段信號(hào)進(jìn)行傅里葉變換，最后得到各段信號(hào)構(gòu)成的時(shí)變譜陣。短時(shí)傅里葉變換對(duì)信號(hào)分析的示意圖如圖8.1所示。圖8.1短時(shí)傅里葉變換示意圖通常對(duì)于短時(shí)傅里葉變換X(n，ω)進(jìn)行時(shí)域抽取，值得注意的是抽取的間隔L必須滿足L≤Nw/2，這樣才能保證信號(hào)不丟失，即可以從時(shí)域信號(hào)恢復(fù)原信號(hào)。抽取后的短時(shí)傅里葉變換為(8.4)8.1.2頻域窗法

如果用W(ω)表示窗函數(shù)w(n)的傅里葉變換，則信號(hào)的STFT還可表示為

(8.5)對(duì)于固定點(diǎn)頻率ω0其實(shí)現(xiàn)結(jié)構(gòu)如圖8.2(a)所示。由圖可以看出，對(duì)每個(gè)頻率值，STFT是時(shí)序列x(n)與窗函數(shù)w(n)ejω0n的卷積，然后與序列e－jω0n相乘。從頻率解釋就是，先將信號(hào)頻譜X(ω)經(jīng)過帶通濾波器，然后下變頻。對(duì)于固定點(diǎn)頻率ω0，還有另一種實(shí)現(xiàn)方法，其結(jié)構(gòu)如圖8.2(b)所示。STFT表示為(8.6)圖8.2短時(shí)傅里葉變換頻域窗法原理圖由圖可以看出，對(duì)每個(gè)頻率值，STFT是序列x(n)

e－jω0n與窗函數(shù)w(n)的卷積，從頻率解釋就是先將信號(hào)頻譜X(ω)下變頻為X(ω+ω0)，然后通過濾波器濾波。

因此，頻域窗法就是將信號(hào)通過具有不同中心頻率、帶寬為Δω的窄帶濾波器W(ω)進(jìn)行濾波。為使其具有較高的頻率分辨率，W(ω)帶外衰減越大越好。8.1.3窗函數(shù)的選取

1.不同窗的選擇

矩形窗具有最窄的主瓣寬度，但同時(shí)也具有最大的旁瓣峰值和最慢的旁瓣峰值衰減速度，矩形窗存在的邊界效應(yīng)會(huì)造成頻譜泄漏，大大影響分析性能;漢寧窗和漢明窗具有較小的旁瓣峰值和較大的旁瓣峰值衰減速度，但主瓣稍寬于矩形窗，在計(jì)算機(jī)短時(shí)傅里葉變換時(shí)，將大大減少頻譜泄漏，具有較優(yōu)的分析性能。

理想的最佳分析窗函數(shù)應(yīng)具有最窄的主瓣寬度、最小的旁瓣峰值和最大的衰減速度，因此，應(yīng)在不同條件下通過比較選擇最合適的分析窗函數(shù)。

2.不同的窗長(zhǎng)度選擇

在選定了分析窗函數(shù)之后，合理地選擇分析窗長(zhǎng)度也有利于分析精度的提高。若窗口過長(zhǎng)，則信號(hào)不能近似看做平穩(wěn)信號(hào)，傅里葉變換將失去作用，也會(huì)加大運(yùn)算量;若窗口過小，又會(huì)丟失信息。因此，慎重選取窗函數(shù)的長(zhǎng)度也是很有必要的。

MATLAB信號(hào)處理工具箱提供了函數(shù)specgram，用于計(jì)算短時(shí)傅里葉變換的頻譜圖，其調(diào)用格式為

B=specgram(a)

B=specgram(a，nfft，F(xiàn)s)

[B，f]=specgram(a，nfft，F(xiàn)s)

[B，f，t]=specgram(a，nfft，F(xiàn)s)

B=specgram(a，nfft，F(xiàn)s，windows)

B=specgram(a，nfft，F(xiàn)s，windows，noverlap)

specgram(a)

B=specgram(a，f，F(xiàn)s，windows，noverlap)

其中，B=specgram(a)為返回信號(hào)向量a的譜圖，該調(diào)用缺省時(shí)，參數(shù)配置如下:

(1)FFT計(jì)算點(diǎn)數(shù)nfft＝min(256，length(a));

(2)采樣頻率Fs＝2;

(3)窗函數(shù)windows＝hanning;

(4)分段重疊的點(diǎn)數(shù)noverlap＝length(windows)/2;當(dāng)a為實(shí)數(shù)時(shí)，specgram只在正頻域內(nèi)計(jì)算離散傅里葉變換。當(dāng)a為偶數(shù)時(shí)，B的行數(shù)為(nfft/2+1);當(dāng)a為奇數(shù)時(shí)，B的行數(shù)為nfft/2，B的列數(shù)為fix((n－noverlap)/(length(windows)－noverlap)。當(dāng)a為復(fù)數(shù)時(shí)，specgram函數(shù)要計(jì)算正、負(fù)頻域的離散傅里葉變換，這時(shí)B為nfft行的復(fù)矩陣，從B中第一點(diǎn)開始，時(shí)間隨列數(shù)線性增加，而頻率從0開始隨行數(shù)線性下降。

B＝specgram(a，nfft)采用指定的nfft點(diǎn)FFT，點(diǎn)數(shù)為2的冪次時(shí)，運(yùn)算速度較快。

[B，f]＝specgram(a，nfft，F(xiàn)s)除了得到譜圖B外，還返回離散傅里葉變換向量f，F(xiàn)s決定了f的頻率分辨率，并不影響B(tài)的結(jié)果。

[B，f，t]=specgram(a，nfft，F(xiàn)s)附加地得到記錄時(shí)域窗函數(shù)與a相交時(shí)刻的向量t。

B=specgram(a，nfft，F(xiàn)s，windows)中windows指定所采用的窗函數(shù)和窗的寬度。窗函數(shù)的寬度應(yīng)該小于或等于nfft，否則補(bǔ)零。

B=specgram(a，nfft，F(xiàn)s，windows，noverlap)中noverlap指定各分段重疊的點(diǎn)數(shù)。

specgram(a)直接用圖像的方式顯示出譜圖。

【例8.1】

利用MATLAB函數(shù)chirp生成兩個(gè)chirp信號(hào)，一個(gè)頻率由小變大，一個(gè)由大變小，再將兩信號(hào)合成得到信號(hào)y(t)，求y(t)的頻譜圖。

MATLAB程序如下:

%MATLABPROGRAM8-1

clc;

t=0:0.001:1.024-.001;

N=1024;

%得到兩個(gè)chirp信號(hào)并相加

y1=chirp(t，0，1，350);

y2=chirp(t，350，1，0);

y=y1+y2;

subplot(211);plot(t，y1);

ylabel(′Chirpsignaly1′)

%求兩個(gè)chirp信號(hào)和的短時(shí)傅里葉變換

[S，F(xiàn)，T]=specgram(y，127，1，hanning(127)，126);

subplot(212);

surf(T/1000，F(xiàn)，abs(S).^2)

view(－80，30);

shadingflat;

colormap(cool);

xlabel('Time');ylabel('Frequency');zlabel('spectrogram');

程序運(yùn)行結(jié)果如圖8.3所示。該圖表明了信號(hào)的能量隨時(shí)間和頻率的分布。圖8.3短時(shí)傅里葉變換的頻譜圖8.1.4時(shí)頻分析的MATLAB實(shí)現(xiàn)

【例8.2】

對(duì)于解析信號(hào)x(t)=2ejπkt2，k=4，0≤T≤5，信號(hào)帶寬fc=kT=20，采樣頻率fs=5fc，采樣點(diǎn)數(shù)N=Tfs，試用STFT分析其特性。

MATLAB程序如下:

%MATLABPROGRAM8-2

clf;

k=4;T=5;

fc=k*T;fs=5*fc;

Ts=1/fs;N=T*fs;

t=[0:N-1];

%解析信號(hào)x(t)

x=2*exp(j*k*pi*(t*Ts).^2);

subplot(221);

%信號(hào)時(shí)域特性

plot(t*Ts，real(x));

xlabel(′t′);ylabel(′x(t)′);title(′x(t)′);

X=fft(x);

X=fftshift(X);

subplot(222);

%信號(hào)頻域特性

plot((t-N/2)*fs/N，abs(X));

xlabel(′f′);ylabel(′|X(w)|′);title(′X(w)′);

Nw=20;

L=Nw/2;

Tn=(N－Nw)/L+1;

nfft=32;

TF=zeros(Tn，nfft);

%短時(shí)傅里葉變換

fori=[1:Tn]

xw=y((i－1)*L+1:i*L+L);

temp=fft(xw，nfft);

temp=fftshift(temp);

TF(i，:)=temp;

end

subplot(223);

fnew=((1:nfft)-nfft/2)*fs/nfft;

tnew=(1:Tn)*L*Ts;

%時(shí)頻分布三維圖繪制

[F，T]=meshgrid(fnew，tnew);

mesh(F，T，abs(TF));

subplot(224);

%繪制曲面的等高線圖

contour(F，T，abs(TF));

程序運(yùn)行結(jié)果如圖8.4所示。圖8.4解析信號(hào)的短時(shí)傅里葉變換頻譜分析

【例8.3】

對(duì)于實(shí)信號(hào)x(t)=2sin(kπt2)，k=4，0≤T≤5，信號(hào)帶寬fc=kT=20，采樣頻率fs=5fc，采樣點(diǎn)數(shù)N=Tfs，試用STFT分析其特性。

MATLAB程序如下:

%MATLABPROGRAM8-3

clf;

k=4;T=5;

fc=k*T;

fs=5*fc;

Ts=1/fs;

N=T*fs;

t=[0:N-1];

%實(shí)信號(hào)x(t)

x=2*sin(k*pi*(t*Ts).^2);

subplot(221);

%信號(hào)時(shí)域特性

plot(t*Ts，real(x));

xlabel(′t′);ylabel(′x(t)′);title(′x(t)′);

X=fft(x);

X=fftshift(X);

subplot(222);

%信號(hào)頻域特性

plot((t-N/2)*fs/N，abs(X));

xlabel(′f′);ylabel(′|X(w)|′);title(′X(w)′);

Nw=20;

L=Nw/2;

Tn=(N－Nw)/L+1;

nfft=32;

TF=zeros(Tn，nfft);

%短時(shí)傅里葉變換

fori=[1:Tn]

xw=x((i－1)*10+1:i*10+10);

temp=fft(xw，nfft);

temp=fftshift(temp);

TF(i，:)=temp;

end

subplot(223);

fnew=((1:nfft)－nfft/2)*fs/nfft;

tnew=(1:Tn)*L*Ts;

%時(shí)頻分布三維圖繪制

[F，T]=meshgrid(fnew，tnew);

mesh(F，T，abs(TF));

subplot(224);

%繪制曲面的等高線圖

contour(F，T，abs(TF));程序運(yùn)行結(jié)果如圖8.5所示。圖8.5實(shí)信號(hào)的短時(shí)傅里葉變換頻譜分析由例8.2和例8.3比較可得，對(duì)于解析信號(hào)只有正頻率，但實(shí)信號(hào)還有負(fù)頻率，時(shí)頻平面的坐標(biāo)為0≤T≤5，頻率范圍為0≤f≤20。

【例8.4】

對(duì)于解析信號(hào)x(t)=2ejπkt2，k=4，0≤T≤5，信號(hào)帶寬fc=kT=20，采樣頻率fs=5fc，采樣點(diǎn)數(shù)N=Tfs。將復(fù)高斯白噪聲信號(hào)加入x(t)后，得到合成信號(hào)y(t)，試用STFT分析y(t)的特性。MATLAB程序如下:

%MATLABPROGRAM8-4

clf;

k=4;T=5;

fc=k*T;fs=5*fc;

Ts=1/fs;N=T*fs;

t=[0:N－1];

%解析信號(hào)x(t)

x=2*exp(j*k*pi*(t*Ts).^2);

%復(fù)高斯白噪聲

noise=randn(1，N)+j*randn(1，N);

%合成輸出信號(hào)y(t)

y=x+0.2*noise;

subplot(221);

%信號(hào)時(shí)域特性

plot(t*Ts，real(y));

xlabel(′t′);ylabel(′y(t)′);title(′y(t)′);

Y=fft(y);

Y=fftshift(Y);

subplot(222);

%信號(hào)頻域特性

plot((t－N/2)*fs/N，abs(Y));

xlabel(′f′);ylabel(′|Y(w)|′);title(′Y(w)′);

Nw=20;

L=Nw/2;

Tn=(N－Nw)/L+1;

nfft=32;

TF=zeros(Tn，nfft);

%短時(shí)傅里葉變換

fori=[1:Tn]

yw=y((i－1)*10+1:i*10+10);

temp=fft(yw，nfft);

temp=fftshift(temp);

TF(i，:)=temp;

end

subplot(223);

fnew=((1:nfft)-nfft/2)*fs/nfft;

tnew=(1:Tn)*L*Ts;

%時(shí)頻分布三維圖繪制

[F，T]=meshgrid(fnew，tnew);

mesh(F，T，abs(TF));

subplot(224);

%繪制曲面的等高線圖

contour(F，T，abs(TF));程序運(yùn)行結(jié)果如圖8.6所示。該結(jié)果表明:當(dāng)有噪聲存在時(shí)，視頻譜變差。圖8.6STFT分析含噪聲的復(fù)合時(shí)變信號(hào)

【例8.5】

有一個(gè)時(shí)變信號(hào)x(n)，0≤n≤N－1，N=400。在n∈[40，90]時(shí)有正弦信號(hào)sin(2πf1t);在n∈[300，360]時(shí)有正弦信號(hào)sin(2πf2t);其余各點(diǎn)信號(hào)值為零。其中，歸一化頻率f1=0.2Hz，f2=0.3Hz。試用STFT分析其特性。

MATLAB程序如下:

%MATLABPROGRAM8-5

clf;

N=400;t=[0:N－1];

x=zeros(1，N);%數(shù)字信號(hào)x(n)

f1=0.2;f2=0.3;

%復(fù)合信號(hào)

x(40:90)=sin(2*pi*f1*(t(40:90)－40));

x(300:360)=sin(2*pi*f2*(t(300:360)－300));

subplot(221);

%信號(hào)時(shí)域特性

plot(t，x);

xlabel(′t′);ylabel(′x(t)′);title(′x(t)時(shí)域特性′);

X=fft(x);

X=fftshift(X);

subplot(222);

%信號(hào)頻域特性

f=(t-N/2)/N;

plot(f，abs(X));

xlabel(′f′);ylabel(′|X(w)|′);title(′X(w)頻域特性′);

Nw=20;

L=Nw/2;

Tn=(N－Nw)/L+1;

nfft=32;

TF=zeros(Tn，nfft);

fori=[1:Tn]

xw=x((i-1)*L+1:i*L+L);

temp=fft(xw，nfft);

temp=fftshift(temp);

TF(i，:)=temp;

end

subplot(223);

fnew=((1:nfft)－nfft/2)/nfft;

tnew=(1:Tn)*L;

[F，T]=meshgrid(fnew，tnew);

mesh(F，T，abs(TF));

subplot(224);

contour(F，T，abs(TF));程序運(yùn)行結(jié)果如圖8.7所示。圖8.7STFT分析時(shí)變信號(hào)

8.2維格納時(shí)頻分布

8.2.1連續(xù)時(shí)間維格納分布

1.維格納分布的定義

設(shè)連續(xù)時(shí)間信號(hào)x(t)定義于整個(gè)時(shí)域，且是復(fù)值的，即x(t)∈C，t∈R。該信號(hào)的自維格納分布定義為

(8.7)是函數(shù)對(duì)τ的傅里葉變換。兩個(gè)連續(xù)時(shí)間信號(hào)x(t)與y(t)的互維格納分布定義為(8.8)當(dāng)不會(huì)引起混淆時(shí)，常將上述兩個(gè)函數(shù)稱為維格納分布。

x(t)與y(t)兩個(gè)信號(hào)的傅里葉譜X(ω)與Y(ω)的維格納分布定義為(8.9)

Wx,y(t，ω)與WX,Y(ω，t)的關(guān)系為

WX,Y(ω，t)=Wx,y(t，ω)

(8.10)式(8.10)表明，兩個(gè)信號(hào)傅里葉譜的維格納分布可由其時(shí)間信號(hào)的維格納分布中頻率和時(shí)間變量互相交換而得到，即維格納分布的時(shí)域和頻域間存在對(duì)稱性。

2.維格納分布的性質(zhì)

1)對(duì)稱性

復(fù)值信號(hào)x(t)的維格納分布是實(shí)函數(shù)，即

Wx(t，ω)=W*x(t，ω)

(8.11)

進(jìn)而可得實(shí)信號(hào)x(t)的維格納分布是頻率的偶函數(shù)，即

Wx(t，ω)=Wx(t，－ω)

(8.12)

2)位移性

(1)x(t)時(shí)移τ，其Wx(t，ω)也時(shí)移τ，即

(8.13)

(2)x(t)為ejΩt調(diào)制，其維格納分布頻移Ω，即

x(t)ejΩt

Wx(t，ω－Ω)

(8.14)

3)定義域的同一性

(1)若x(t)限制在t1≤t≤t2中，則其Wx(t，ω)也被限制在同一時(shí)域范圍中;

(2)若X(ω)限制在ω1≤ω≤ω2中，則其Wx(t，ω)也被限制在同一時(shí)域范圍中;

(3)若x(t)是因果信號(hào)，則其Wx(t，ω)也是因果信號(hào)的，即

Wx(t，ω)=0，t<0

(8.15)

解析信號(hào)xa(t)的Wx(t，ω)限制在上半平面，即

Wx(t，ω)=0，ω<0

(8.16)

4)反演性

(1)某一時(shí)刻t的值x(t)，可以通過時(shí)刻等于t/2處將Wx(t/2，ω)對(duì)頻率ω作傅里葉反變換得到，其中只差一個(gè)比例函數(shù)x*(0)，即

(8.17)

(2)頻率ω的值X(ω)，可以通過頻率等于ω/2處將Wx(t，ω/2)對(duì)時(shí)間t作傅里葉變換得到，其中只差一個(gè)比例系數(shù)X*(0)，即(8.18)這一性質(zhì)可用來由維格納分布恢復(fù)原信號(hào)x(t)。

5)積分性

(1)在固定時(shí)刻t，Wx(t，ω)沿全頻軸的積分等于該時(shí)刻x(t)的瞬時(shí)功率|x(t)|2，即

(8.19)

(2)在固定頻率ω，Wx(t，ω)沿全時(shí)軸的積分等于該頻率的能譜密度|X(ω)|2，即(8.20)

(3)由式(8.19)和式(8.20)可推出:Wx(t，ω)沿時(shí)、頻兩個(gè)軸的雙重積分等于信號(hào)的能量E，即(8.21)

6)基本運(yùn)算

(1)加法運(yùn)算，若x(t)=x1(t)+x2(t)，則

(8.22)由于維格納分布不是線性運(yùn)算，引入交叉項(xiàng)干擾

，從而使得對(duì)各分量信號(hào)維格納分布的直觀解釋發(fā)生困難。

(2)卷積運(yùn)算，若

，則與其對(duì)應(yīng)的維格納分布也在時(shí)軸上褶積，即(8.23)

(3)乘法運(yùn)算，若y(t)=x(t)m(t)，且對(duì)應(yīng)的頻譜X(ω)與M(ω)在頻域上褶積，則和它們對(duì)應(yīng)的維格納分布也將在頻域上褶積，即(8.24)

3.偽維格納分布

維格納分布可作為能量分布來表示信號(hào)的瞬時(shí)特征，但它是在全時(shí)軸(－∞<t<∞)上定義的，因此不便于實(shí)時(shí)處理。

實(shí)際工作中只能取有限長(zhǎng)的數(shù)據(jù)來進(jìn)行分析，這相當(dāng)于對(duì)原始信號(hào)施加一個(gè)隨時(shí)間滑動(dòng)的窗函數(shù)。如果用t表示隨研究時(shí)刻而滑動(dòng)的窗位置，并將時(shí)間變量改用τ表示，則被處理的信號(hào)為

xt(τ)=x(τ)·w(τ－t)

(8.25)

式中，w(t－τ)通常是以t為中心而且對(duì)稱的函數(shù)。

根據(jù)維格納分布的乘法運(yùn)算公式有

(8.26)其中，Ww(t，ω)為窗函數(shù)w(t)的維格納分布，即(8.27)通常不用考慮式(8.25)中的全部，而只需了解窗函數(shù)滿足t=τ處的維格納分布，稱為信號(hào)x(t)的偽維格納分布，簡(jiǎn)記作PWx(t，ω)，即(8.28)因?yàn)樗m然在形式上很像維格納分布，但實(shí)際已不是信號(hào)x(t)在原始意義下的維格納分布了。式(8.28)說明，PWx(t，ω)是Wx(t，ω)在頻域上被Ww(0，ω)褶積而得到的，故加滑動(dòng)窗的結(jié)果是:在時(shí)域上把數(shù)據(jù)截短，而在頻域上對(duì)Wx(t，ω)起平滑作用。頻域的平滑會(huì)降低在頻率軸方向上的分辨率。因此，在前述維格納分布的諸性質(zhì)中，有關(guān)時(shí)間方面的各項(xiàng)性質(zhì)仍然保持，但有關(guān)頻率方面的各項(xiàng)性質(zhì)則受影響，需重新考慮。

4.維格納分布與模糊函數(shù)的關(guān)系

模糊函數(shù)(AF)也是一種常用的時(shí)頻表示，廣泛應(yīng)用于雷達(dá)、聲納等領(lǐng)域，并在很多書中已有詳述。本節(jié)側(cè)重討論它與維格納分布的不同點(diǎn)及相互間的關(guān)系，這有助于了解后面將要介紹的信號(hào)時(shí)頻分布的統(tǒng)一表示。

為了研究的需要，定義信號(hào)x(t)在時(shí)域的瞬時(shí)自相關(guān)函數(shù)為(8.29)設(shè)信號(hào)x(t)的頻譜為X(ω)，定義其在頻域的瞬時(shí)譜相關(guān)函數(shù)為(8.30)其維格納分布為(8.31)其模糊函數(shù)為(8.32)

5.信號(hào)時(shí)頻分布統(tǒng)一表述

在眾多的時(shí)頻表示中，維格納分布和模糊函數(shù)是實(shí)際中用得較廣的兩種時(shí)頻表示法。盡管它們的函數(shù)形式和性質(zhì)各不相同，但具有一定的聯(lián)系和共性:具有特定性質(zhì)的時(shí)頻分布是對(duì)核函數(shù)施以一定約束條件得到的。因而信號(hào)時(shí)頻分布可以統(tǒng)一表述。

許多常用的二次型能量化時(shí)頻表示都可以由包含核函數(shù)的廣義雙線性時(shí)頻表示得出，即(8.33)式中，j(ξ，τ)表示核函數(shù)，它決定時(shí)頻分布Px(t，f)的特性。采用不同的核函數(shù)，將得到不同的時(shí)頻分布。8.2.2離散時(shí)間維格納分布

1.離散時(shí)間的維格納分布定義

該定義是連續(xù)時(shí)間信號(hào)維格納分布定義的直接引申。將連續(xù)時(shí)間信號(hào)維格納分布定義寫成離散形式即得到離散時(shí)間信號(hào)的維格納分布變換(8.34)令k′=2k，得(8.35)對(duì)應(yīng)的頻域定義為

(8.36)

由上式可見，頻域的重復(fù)周期是π，即(8.37)

2.離散時(shí)間的維格納分布性質(zhì)

離散時(shí)間維格納分布定義可保持連續(xù)時(shí)間信號(hào)維格納分布定義中有關(guān)時(shí)間上的一些特性，但頻率上的一些特性卻被破壞了。

(1)頻域總和為

(8.38)這說明Wx(n，ω)的周期為π，但也反映了頻域存在混淆現(xiàn)象。

(2)頻域的反演公式為(8.39)

3.離散時(shí)間的維格納分布計(jì)算

維格納分布的計(jì)算量大，可用FFT來計(jì)算，但目前的各種快速算法還不能從根不上解決其計(jì)算量大的問題，距實(shí)際應(yīng)用尚有一定距離。本節(jié)主要介紹FFT實(shí)現(xiàn)的方法。

計(jì)算離散時(shí)間信號(hào)x(n)的維格納分布，需對(duì)定義式(8.35)作截?cái)嗉哟疤幚?，加窗函?shù)后的維格納分布為(8.40)式中，w(l)為時(shí)寬(2L－1)的窗函數(shù)(當(dāng)|l|>L時(shí)，w(l)=0)。為了能用基2FFT來計(jì)算維格納分布，需在頻域?qū)x(n，ω)計(jì)算采樣值，由于Wx(n，ω)在頻域的重復(fù)周期為π，若一周期內(nèi)的采樣點(diǎn)為N，則采樣間隔Δω=π/N。此外，還應(yīng)采用補(bǔ)零的方法，使得N＝2L，以方便計(jì)算。

(8.41)由此可得(8.42)由于FFT通常是在(0~L－1)范圍內(nèi)計(jì)算的，故可用如下的重排序列來實(shí)現(xiàn)，即令(8.43)代入式(8.42)得(8.44)對(duì)于MATLAB6.x、MATLAB7.x以及MATLABR2008版來說，此函數(shù)不存在，wig2函數(shù)的源代碼如下:

function[wx，waxis]=wig2(x0，nfft，flag)

%---------------parameterchecks-------------

[m，n]=size(x0);

if(min(m，n)~=1)

disp([′wig2:inputargumentxisa′，int2str(m)，′by′，int2str(n)，...′array′])

error(′Inputargumentxmustbeavector′);

end

if(exist(′flag′)~=1)flag=1;end

if(all(imag(x0)==0)&flag~=0)x0=hilbert(x0);end

%----------findpoweroftwoforFFT-------------

lx=length(x0);

lfft=2^nextpow2(2*lx);%minimumFFTlength

if(exist(′nfft′)~=1)nfft=lfft;end

if(isempty(nfft))nfft=lfft;end

if(nfft<2*lx)

disp([′WIG2:FFTlengthmustexceedtwicethesignallength′])

disp([′resettingFFTlengthto′，int2str(lfft)])

nfft=lfft;

end

x=zeros(nfft，1);x(1:lx)=x0(:);cx=conj(x);

wx=zeros(nfft，lx);

L1=lx-1;

%-------computer(tau，t)=cx(t-tau/2)*x(t+tau/2)------

forn=0:L1

indm=max(－n，－L1+n):min(n，L1－n);

indy=indm+(indm<0)*nfft;%outputindicesy(m;n)

y=zeros(nfft，1);

y(indy+1)=x(n+indm+1).*cx(n-indm+1);

wx(:，n+1)=y;

end

%---------WD(f，t)=FT(tau-->f)r(tau，t)--------

wx=fft(wx);

wx=real(wx.′);%forceittobereal

%--------displaytheWD-------------

nfftby2=nfft/2;

wx=wx(:，[nfftby2+1:nfft，1:nfftby2]);

waxis=[－nfftby2:nfftby2－1]/(2*nfft);

taxis=1:lx;

contour(waxis，taxis，abs(wx)，8);

xlabel(′frequency′);ylabel(′timeinsamples′);

title(′WS′);grid;

set(gcf，′Name′，′HosaWIG2′)

MATLAB5.3版的信號(hào)處理工具箱提供了函數(shù)wig2c，用于計(jì)算基于Choi-Williams核的維格納分布，調(diào)用格式為

[wx，waxis]=wig2c(x，nfft，sigma，flag)

其中，wig2c函數(shù)返回信號(hào)向量x的維格納分布wx，

waxis為相應(yīng)的頻率點(diǎn)向量。nfft為FFT的計(jì)算點(diǎn)數(shù)，為了避免混疊，nfft值應(yīng)大于等于信號(hào)x長(zhǎng)度的兩倍，否則，采用缺省值，即大于x長(zhǎng)度兩倍且為2的冪次的最小數(shù)。sigma為Choi-Williams核的參數(shù)，隨著sigma的增大，維格納分布抑制交叉項(xiàng)的能力減弱，其缺省值為0.05。如果flag為非零值，且x為實(shí)信號(hào)，則使用x的解析信號(hào)作變換，flag缺省值為1。在高版本的MATLAB中，此函數(shù)不存在。8.2.3時(shí)頻分布的MATLAB實(shí)現(xiàn)

為了分析維格納時(shí)頻分布的原理，在利用MATLAB編程時(shí)，未采用MATLAB5.3版信號(hào)處理工具箱中的內(nèi)部函數(shù)wig2和wig2c，而直接編寫仿真程序。

【例8.6】

對(duì)于解析信號(hào)x(t)=2ejπkt2，k=4，0≤T≤4，信號(hào)帶寬fc=kT=16，采樣頻率fs=4fc，采樣點(diǎn)數(shù)N=Tfs，試用維格納分布分析其特性。

MATLAB程序如下:

%MATLABPROGRAM8-6

clf;

k=4;T=4;

fc=k*T;fs=4*fc;

Ts=1/fs;N=T*fs;t=[0:N－1];

%解析信號(hào)x(t)

x=2*exp(j*k*pi*(t*Ts).^2);

subplot(221);

%信號(hào)時(shí)域特性

plot(t*Ts，real(x));

xlabel(′t′);ylabel(′x(t)′);legend(′x(t)′);

X=fft(x);

X=fftshift(X);

subplot(222);

%信號(hào)頻域特性

plot((t-N/2)*fs/N，abs(X));

xlabel(′f′);ylabel(′|X(w)|′);legend(′X(w)′);

R=zeros(N，N);

%維格納分布

forn=[0:N－1]

M=min(n，N－1－n);

fork=[0:M]

R(n+1，k+1)=x(n+k+1)*conj(x(n－k+1));

end

fork=[N－1:－1:N－M]

R(n+1，k+1)=conj(R(n+1，N－k+1));

end

end

TF=zeros(N，N);

forn=[0:N－1]

temp=fft(R(n+1，:));

temp=fftshift(temp);

TF(n+1，:)=temp;

end

fnew=(t－N/2)*fs/2/N;

tnew=(0:N－1)*Ts;

[F，T]=meshgrid(fnew，tnew);

subplot(223);

mesh(F，T，abs(TF));

subplot(224);

contour(F，T，abs(TF));

程序運(yùn)行結(jié)果如圖8.8所示。圖8.8解析信號(hào)的維格納分布圖(fs=4fc)

【例8.7】采用例8.6中的解析信號(hào)，僅將采樣頻率變?yōu)閒s=2fc，試用維格納分布分析其特性。

MATLAB程序如下:

%MATLABPROGRAM8-7

clf;

k=4;T=4;

fc=k*T;fs=2*fc;

Ts=1/fs;N=T*fs;

t=[0:N－1];

x=2*exp(j*k*pi*(t*Ts).^2);

subplot(221);

plot(t*Ts，real(x));

xlabel(′t′);ylabel(′x(t)′);legend(′x(t)′);

X=fft(x);

X=fftshift(X);

subplot(222);

plot((t－N/2)*fs/N，abs(X));

xlabel(′f′);ylabel(′|X(w)|′);legend(′X(w)′);

R=zeros(N，N);

forn=[0:N－1]

M=min(n，N－1－n);

fork=[0:M]

R(n+1，k+1)=x(n+k+1)*conj(x(n－k+1));

end

fork=[N－1:－1:N－M]

R(n+1，k+1)=conj(R(n+1，N－k+1));

end

end

TF=zeros(N，N);

forn=[0:N－1]

temp=fft(R(n+1，:));

temp=fftshift(temp);

TF(n+1，:)=temp;

end

fnew=(t－N/2)*fs/2/N;

tnew=(0:N－1)*Ts;

[F，T]=meshgrid(fnew，tnew);

subplot(223);

mesh(F，T，abs(TF));

subplot(224);

contour(F，T，abs(TF));

程序運(yùn)行結(jié)果如圖8.9所示。圖8.9解析信號(hào)的維格納分布圖(fs=2fc)從例8.6的輸出時(shí)頻平面中可以清晰地看出，頻率隨時(shí)間線性變換，其變化率k為4，可見整個(gè)分析是正確的。從例8.7的輸出結(jié)果可以看出，對(duì)于帶寬為fc的模擬信號(hào)進(jìn)行FFT分析，采樣頻率至少是截止頻率的兩倍，即fs=2fc，這就是奈奎斯特準(zhǔn)則的要求。然而對(duì)于維格納分布變換，采樣頻率至少是信號(hào)截止頻率的四倍，才能保證不會(huì)混淆，即fs=4fc。

【例8.8】

對(duì)于解析信號(hào)x(t)=2ejπkt2，k=4，0≤T≤4，信號(hào)帶寬fc=kT=16，采樣頻率fs=5fc，采樣點(diǎn)數(shù)N=Tfs。在信號(hào)中加入復(fù)高斯白噪聲，用維格納分布分析其特性。MATLAB程序如下:

%MATLABPROGRAM8-8

clf;

k=4;T=4;fc=k*T;

fs=5*fc;

Ts=1/fs;N=T*fs;

t=[0:N－1];

x=zeros(1，N);

x=2*exp(j*k*pi*(t*Ts).^2);

noise=randn(1，N)+j*randn(1，N);

x=x+0.2*noise;

subplot(221);

plot(t*Ts，real(x));

xlabel(′t′);ylabel(′x(t)′);legend(′x(t)′);

X=fft(x);

X=fftshift(X);

subplot(222);

plot((t－N/2)*fs/N，abs(X));

xlabel(′f′);ylabel(′|X(w)|′);legend(′X(w)′);

R=zeros(N，N);

%維格納分布

forn=[0:N－1]

M=min(n，N－1－n);

fork=[0:M]

R(n+1，k+1)=x(n+k+1)*conj(x(n－k+1));

end

fork=[N－1:－1:N－M]

R(n+1，k+1)=conj(R(n+1，N－k+1));

end

end

TF=zeros(N，N);

forn=[0:N－1]

temp=fft(R(n+1，:));

temp=fftshift(temp);

TF(n+1，:)=temp;

end

fnew=(t－N/2)*fs/2/N;

tnew=(0:N－1)*Ts;

[F，T]=meshgrid(fnew，tnew);

subplot(223);

mesh(F，T，abs(TF));

subplot(224);

contour(F，T，abs(TF));

程序運(yùn)行結(jié)果如圖8.10所示。該例說明了維格納分布變換對(duì)噪聲信號(hào)不太敏感，時(shí)頻變換后信噪比較高。圖8.10含噪聲信號(hào)的維格納分布

【例8.9】

數(shù)字信號(hào)x(n)，0≤n≤N－1，N=400，在n∈[40，90]時(shí)有正弦信號(hào)sin(2πf1t);在n∈[280，320]時(shí)有正弦信號(hào)sin(2πf2t);其余各點(diǎn)信號(hào)值為零。其中，歸一化頻率f1=0.2Hz，f2=0.3Hz。用維格納分布分析其特性。

MATLAB程序如下:

%MATLABPROGRAM8-9

clf;

N=400;

t=[0:N－1];

%數(shù)字信號(hào)x(n)

x=zeros(1，N);

f1=0.2;f2=0.3;

%復(fù)合信號(hào)

x(40:90)=sin(2*pi*f1*(t(40:90)－40));

x(280:320)=sin(2*pi*f2*(t(280:320)－280));

subplot(221);

plot(t，x);

xlabel(′t′);ylabel(′x(t)′);legend(′x(t)′);

X=fft(x);

X=fftshift(X);

subplot(222);

plot((t-N/2)/N，abs(X));

xlabel(′f′);ylabel(′|X(w)|′);legend(′X(w)′);

R=zeros(N，N);

forn=[0:N－1]

M=min(n，N－1－n);

fork=[0:M]

R(n+1，k+1)=x(n+k+1)*conj(x(n－k+1));

end

fork=[N-1:－1:N－M]

R(n+1，k+1)=conj(R(n+1，N－k+1));

end

end

TF=zeros(N，N);

forn=[0:N－1]

temp=fft(R(n+1，:));

temp=fftshift(temp);

TF(n+1，:)=temp;

end

fnew=(t－N/2)/N;

tnew=[0:N－1];

[F，T]=meshgrid(fnew，tnew);

subplot(223);

mesh(F，T，abs(TF));

subplot(224);

contour(F，T，abs(TF));

程序運(yùn)行結(jié)果如圖8.11所示。圖8.11數(shù)字信號(hào)的維格納分布

8.3小波變換

小波變換提出了變化的時(shí)間窗，當(dāng)需要精確的低頻信息時(shí)，采用長(zhǎng)的時(shí)間窗;當(dāng)需要精確的高頻信息時(shí)，采用短的時(shí)間窗。小波變換用的不是時(shí)間-頻率域，而是時(shí)間-尺度域。尺度越大，采用越長(zhǎng)的時(shí)間窗;尺度越小，采用越短的時(shí)間窗，即尺度與時(shí)間成反比。

假設(shè)Ψ(t)為平方可積函數(shù)，即Ψ(t)∈L2(R)，若其傅里葉變換為，并滿足完全重構(gòu)條件

(8.45)則稱Ψ(t)為一個(gè)基小波或母小波。母函數(shù)Ψ(t)經(jīng)伸縮和平移后為(8.46)則稱Ψa,b(t)為一個(gè)小波序列。其中，a為伸縮因子，b為平移因子。

Ψa,b(t)的傅里葉變換為(8.47)8.3.1連續(xù)小波變換

連續(xù)小波變換(CWT)又稱積分小波變換(IWT)。任意函數(shù)f(t)的連續(xù)小波變換定義為

(8.48)其重構(gòu)公式(逆變換ICWT)為(8.49)由于基小波Ψ(t)生成的小波Ψa,b(t)，在小波變換中對(duì)被分析的信號(hào)起著觀測(cè)窗的作用，所以，Ψ(t)還應(yīng)該滿足一般函數(shù)的約束條件(8.50)故是一個(gè)連續(xù)函數(shù)。為滿足完全重構(gòu)條件，在原點(diǎn)的值必須為0。為了使信號(hào)重構(gòu)的實(shí)現(xiàn)在數(shù)值上是穩(wěn)定的，小波Ψ(t)的傅里葉變化需滿足穩(wěn)定性條件(8.51)8.3.2離散小波變換

在實(shí)際運(yùn)用中，尤其是在計(jì)算機(jī)上實(shí)現(xiàn)時(shí)，連續(xù)小波必須加以離散化。因此，下面討論連續(xù)小波Ψa,b(t)和連續(xù)小波變換Wf(a，b)的離散化。需要強(qiáng)調(diào)的是，離散化都是針對(duì)連續(xù)的尺度參數(shù)a和連續(xù)平移參數(shù)b的，而不是針對(duì)時(shí)間變量t的。這與時(shí)間離散化不同。

在連續(xù)小波中，函數(shù)Ψa,b(t)為

(8.52)式中，b∈R，a∈R+，且a≠0，Ψ是容許的。在離散化中，限制a只取正值，這樣相容性條件就變?yōu)?8.53)將連續(xù)小波變換中尺度參數(shù)a和平移參數(shù)b的離散公式，分取為a=aj0，b=kaj0b0，j∈Z，擴(kuò)展步長(zhǎng)a0≠1，是固定值,且一般假定a0>1。因此，對(duì)應(yīng)的離散小波函數(shù)Ψj,k(t)為(8.54)

離散化小波變換系數(shù)表示為(8.55)重構(gòu)公式為(8.56)式中，C是一個(gè)與信號(hào)無關(guān)的常數(shù)。然而，為了保證重構(gòu)信號(hào)的精度，在選擇a0和b0時(shí)，網(wǎng)格點(diǎn)應(yīng)盡可能密(即a0和b0盡可能小)。因?yàn)榫W(wǎng)格點(diǎn)越稀疏，小波函數(shù)Ψj,k(t)和離散小波系數(shù)Cj,k就越少，信號(hào)重構(gòu)的精確度也就會(huì)越低。8.3.3小波變換在突變信號(hào)檢測(cè)中的應(yīng)用

設(shè)j(t)為具有低通性質(zhì)的平滑函數(shù)，以它的一階和二階導(dǎo)數(shù)作為小波，對(duì)f(t)作小波變換，可以得到(8.57)(8.58)式中，“”表示卷積，j1(t)和j2(t)表示j(t)的一階和二階導(dǎo)數(shù)，對(duì)f(t)作小波變換，相當(dāng)于f(t)被j(t)平滑后，再對(duì)t求一階或二階導(dǎo)數(shù)。對(duì)某一固定a值，f(t)j

(t)的拐點(diǎn)既是Wj1f(a，b)的極值點(diǎn)，又是Wj1f(a，b)的過零點(diǎn)，由此可檢測(cè)出信號(hào)的急劇變化。選用具有低通性質(zhì)的高斯函數(shù)作為平滑函數(shù)，以其一階導(dǎo)數(shù)、二階導(dǎo)數(shù)(墨西哥草帽小波函數(shù))作為小波基函數(shù)進(jìn)行突變點(diǎn)分析。各小波函數(shù)的表達(dá)式為

(8.59)上述小波函數(shù)具有對(duì)稱、可微及可積的性質(zhì)，在時(shí)頻兩域都為高斯型且呈平方指數(shù)衰減特性，具有很好的時(shí)頻局域性，用作小波變換可準(zhǔn)確地識(shí)別信號(hào)突變點(diǎn)。

【例8.10】

利用MATLAB編寫程序，繪制500s長(zhǎng)的小波函數(shù)j(t，a)。小波基位于250s處，尺度步長(zhǎng)為10。MATLAB程序如下:

%MATLABPROGRAM8-10

Ls=500;t0=250;a=10;

t=(1:1:Ls);%給出時(shí)間序列

x=exp(－((t－t0)/a).^2/2)/a;

%高斯小波基函數(shù)

y=(t-t0).*exp(－((t－t0)/a).^2/2)/a;%一階導(dǎo)數(shù)

%小波變換二階導(dǎo)數(shù)墨西哥草帽(MexscroHat)

z=(1－((t－t0)/a).^2).*exp(－(