現(xiàn)代信號處理方法1-1

現(xiàn)代信號處理方法主講：

王忠仁2012年2月授課內(nèi)容〇、MATLAB語法要點(diǎn)

一、信號的時頻分析二、Radon-Wigner變換三、分?jǐn)?shù)階Fourier變換四、高階譜估計參考書[1]

張賢達(dá).現(xiàn)代信號處理.清華大學(xué)出版社.1995[2]

張賢達(dá),保錚.非平穩(wěn)信號分析與處理.國防工業(yè)出版社.1998[3]

張賢達(dá).時間序列分析-高階統(tǒng)計量方法.清華大學(xué)出版社.1996[4]

王忠仁,張靜.復(fù)變函數(shù)與積分變換.高等教育出版社.20060．１MATLAB語言特點(diǎn)向量與矩陣的運(yùn)算非常方便庫函數(shù)極其豐富語言簡潔，程序簡短圖形功能強(qiáng)大比C、FORTRAN語言更高級的第四代計算機(jī)語言：機(jī)器語言匯編語言C-語言MATLAB第〇章

MATLAB語法要點(diǎn)

0．２MATLAB語法要點(diǎn)算術(shù)運(yùn)算符：標(biāo)量運(yùn)算:

加減乘除冪+-*/^矩陣運(yùn)算:

加減乘右除左除冪轉(zhuǎn)置復(fù)共軛轉(zhuǎn)置

+-*/\^.’’

若Ax=b則x=A\b(常用);若xA=b則x=A/b(不常用)

標(biāo)量情況有:2/5=0.4;2\5=2.5矩陣對應(yīng)元素運(yùn)算:乘右除左除冪

.*./.\.^

數(shù)乘矩陣:c*A冒號運(yùn)算符:a:h:b是以a為初值，h為步長,終值不超過b（但再增加一個步長就超過b）的行向量關(guān)系運(yùn)算符:小于小于等于大于大于等于等于不等于<<=>>===~=邏輯運(yùn)算符:與或非&|~注釋與語句分隔:%開始注釋,直到行末;;分隔語句,阻止前面變量的結(jié)果顯示在命令窗口；,分隔語句,前面變量的結(jié)果可顯示在命令窗口。幫助：在命令窗口鍵入help函數(shù)名循環(huán)結(jié)構(gòu):

(1)for語句格式:fort=表達(dá)式1:表達(dá)式2:表達(dá)式3

語句體

end其中:表達(dá)式1為循環(huán)初值;

表達(dá)式2為步長,省略時默認(rèn)步長為1;

表達(dá)式3為循環(huán)終值.例

forn=1:5form=1:nr(n,m)=m*n;endendr=1000024000369004812160510152025特別注意：MATLAB數(shù)組下標(biāo)不能為0或負(fù)數(shù)！

(2)while語句格式:while表達(dá)式語句體

end其中:當(dāng)表達(dá)式值為真時執(zhí)行語句體;

當(dāng)表達(dá)式值為假時終止循環(huán).例:找出階乘超過10100的最小整數(shù).

n=1;whileprod(1:n)<1e100n=n+1;end注:函數(shù)prod(x)計算向量x的各元素的積,1:n—冒號運(yùn)算省略步長1形成常向量.結(jié)果:n=70階乘1.1979e+100選擇(分支)結(jié)構(gòu):

if語句格式:if邏輯表達(dá)式1

語句體1

elseif

邏輯表達(dá)式2

語句體2

elseif

邏輯表達(dá)式3

語句體3…else

語句體elseendt=0:0.001:0.18;fori=1:length(t)

ift(i)<=0.06x(i)=sin(2*pi*60*t(i));elseift(i)>0.06&t(i)<=0.12x(i)=sin(2*pi*800*t(i)*t(i));elsex(i)=sin(2*pi*100*t(i));endend

0．３幾個MATLAB常用函數(shù)繪圖函數(shù)

plot(t,x)繪制函數(shù)x(t)的曲線圖;plot(t,x,’.’)繪制函數(shù)x(t)的散點(diǎn)圖(點(diǎn)圖還有‘o’,‘*’等);stem(t,x)繪制函數(shù)x(t)的火柴桿圖(空心);stem(t,x,’.’)繪制函數(shù)x(t)的火柴桿圖(實(shí)心).FFT

fft(x(t),n)傅立葉變換

Ifft(X(ω),n)傅立葉反變換第一章信號的時頻分析1.1引言1.2Hilbert變換與解析信號1.3時頻分布及其性質(zhì)1.4二次型時頻分布的交叉項(xiàng)1.5Wigner-Ville分布及其應(yīng)用1.1

引言

設(shè)信號s(t)為能量有限信號，即s(t)滿足數(shù)學(xué)上表明s(t)為平方可積函數(shù)，記。

s(t)的Fourier變換為

(1.1)

1.1.1Fourier變換回顧

S(ω)又稱為s(t)的頻譜，

S(ω)的模|S(ω)|稱為振幅譜，

S(ω)的輻角argS(ω)稱為相位譜。

S(ω)的Fourier逆變換為

Fourier變換可以用來對信號進(jìn)行頻譜分析或頻率濾波。(1.2)例1.1用MATLAB編程形成雷克子波,并進(jìn)行頻譜分析.%峰值頻率為fp的雷克子波的表達(dá)式:%Rc(t)=[1-2(πfpt)2]exp[-(πfpt)2]clccleartm=6.4;dt=0.1;t=-tm+dt:dt:tm;n=length(t);fp=0.5;tt=pi*fp*t;%數(shù)乘向量,pi是MATLAB定義的特殊變量πtt2=tt.*tt;%向量對應(yīng)元素相乘得新向量tt2，或tt.^2rc=(1-2*tt2).*exp(-tt2);xf=fft(rc,length(t));%傅立葉變換df=1.0/(n*dt);fend=n*df/2;w2=0:df:fend;fori=1:length(w2)%傅立葉變換的前一半數(shù)據(jù)取模

xfw(i)=abs(xf(i));endfigure,subplot(2,1,1)plot(t,rc)subplot(2,1,2)plot(w2,xfw)圖1.1時間信號及其振幅譜由s(t)的Fourier變換

可知，S(ω)是由s(t)所有時間加權(quán)累積得到的，

S(ω)的模|S(ω)

|只能整體上反映s(t)主要含有哪些頻率成分（或頻帶）。

在現(xiàn)代信號處理中常常要了解時變信號的時頻局部變化特征，F(xiàn)ourier變換無能為力。

1.1.2Fourier變換的局限性圖1.2時變信號及其振幅譜

信號的時域表示或頻域表示構(gòu)成了觀察一個信號的兩種方式，但它并不能告訴我們某種頻率分量發(fā)生在哪些時間內(nèi)。時－頻分析法是處理非平穩(wěn)信號的有力工具，它是把一維信號映射到二維時－頻平面上，這樣可以在時－頻域里反映信號的非平穩(wěn)特性。時－頻分析基本思想是：設(shè)計時間和頻率的聯(lián)合函數(shù)，用它同時描述信號在不同時間和頻率的能量密度或強(qiáng)度。

1.1.3聯(lián)合時-頻表示時間域頻率域時頻平面圖1.3聯(lián)合時-頻表示平面

圖中左邊的是時間域信號，由它可知強(qiáng)度隨時間的變化規(guī)律；最下面的是能量密度頻譜，由它可知信號存在那些頻率及其相對強(qiáng)度（大約175到340Hz，320Hz最強(qiáng)）；右圖是時間信號的聯(lián)合時－頻分布圖，通過它可以明顯地看出信號的某一頻率在什么時刻存在。如通過右圖可以知道頻率大約從175Hz開始，大至在0.8秒左右的時間內(nèi)線性地增加到大約340Hz，然后在340Hz的地方持續(xù)了大約0.2秒，最后又在約0.4秒的時間內(nèi)線性地減小至約240Hz，再比如還可以知道300Hz的頻率分別在0.8秒時刻和1.6秒時刻兩次出現(xiàn)。

比較著名的時-頻分析方法有：短時Fourier變換法Gabor變換法小波變換法時-頻分布理論

*本章討論時-頻分布理論

1.1.4時-頻分析方法1.2Hilbert變換與解析信號1.2.1Hilbert變換的定義設(shè)有實(shí)信號，它的Hilbert變換記作或并定義為

即式中

表示取積分的主值。(1.3)由兩函數(shù)f(t)與g(t)卷積定義可將Hilbert變換(1.3)式寫成(1.4)注：卷積滿足交換律。例1.2求的Hilbert變換.解:注:狄利克雷積分同理:1.2.2Hilbert變換的物理意義由實(shí)信號的Hilbert變換(1.4)式

可見的Hilbert變換是的一種濾波,濾波器的脈沖響應(yīng)為(1.5)濾波器的

頻率響應(yīng)為其中由卷積定理，Hilbert變換（1.4）式的頻域形式為(1.6)當(dāng)時，(1.7)(1.8)由（1.7）和（1.8）式可見，信號經(jīng)Hilbert變換后，其振幅譜不變，相位譜相差900，因此Hilbert變換又稱為900移相器。1.2.3Hilbert逆變換由Hilbert變換頻域形式（1.6）式又由卷積定理，則有以乘上式兩端，注意到則有于是，Hilbert逆變換為還有(1.9)其中1.2.4解析函數(shù)的定義設(shè)實(shí)信號的Hilbert變換為，稱復(fù)信號為對應(yīng)的解析信號。由（1.7）式有由(1.10)式可見，解析信號只在正頻域上擁有單邊譜。(1.10)1.2.5實(shí)信號的瞬時參數(shù)設(shè)實(shí)信號的解析復(fù)信號為

的瞬時振幅（也叫包絡(luò)）的瞬時相位的瞬時頻率(Hz)

1.2.6解析信號的快速數(shù)值計算

（包括Hilbert變換）1.取實(shí)信號s(t)的N個離散值并做FFT（快速Fourier變換）得到信號的頻域表示2.構(gòu)造3.計算,其中IFFT為FFT的反變換。這樣即得到原信號的解析信號。4.為原信號的Hilbert變換。

注：取N為2的整數(shù)次冪，補(bǔ)零個數(shù)大于原數(shù)據(jù)個數(shù)。例1.3用MATLAB編程計算實(shí)信號的Hilbert變換與包絡(luò).%clearall;clc;fp=24;a=-10;dt=0.002;TT=0.5;%t=dt:dt:TT;t=0:dt:TT;b=exp(a*t).*sin(2*pi*fp*t);nsmp=length(t);N=512;fork=nsmp+1:Nb(k)=0.0;endxf=fft(b,N);fori=1:Nifi==1

xf(i)=xf(i);

elseifi>1&i<=N/2

xf(i)=2*xf(i