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

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

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

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

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

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

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

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

MATLAB語(yǔ)法要點(diǎn)

0．２MATLAB語(yǔ)法要點(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矩陣對(duì)應(yīng)元素運(yùn)算:乘右除左除冪

.*./.\.^

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

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

語(yǔ)句體

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

表達(dá)式2為步長(zhǎng),省略時(shí)默認(rèn)步長(zhǎng)為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語(yǔ)句格式:while表達(dá)式語(yǔ)句體

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

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

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

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

語(yǔ)句體1

elseif

邏輯表達(dá)式2

語(yǔ)句體2

elseif

邏輯表達(dá)式3

語(yǔ)句體3…else

語(yǔ)句體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．３幾個(gè)MATLAB常用函數(shù)繪圖函數(shù)

plot(t,x)繪制函數(shù)x(t)的曲線(xiàn)圖;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)傅立葉反變換第一章信號(hào)的時(shí)頻分析1.1引言1.2Hilbert變換與解析信號(hào)1.3時(shí)頻分布及其性質(zhì)1.4二次型時(shí)頻分布的交叉項(xiàng)1.5Wigner-Ville分布及其應(yīng)用1.1

引言

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

s(t)的Fourier變換為

(1.1)

1.1.1Fourier變換回顧

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

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

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

S(ω)的Fourier逆變換為

Fourier變換可以用來(lái)對(duì)信號(hào)進(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;%向量對(duì)應(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時(shí)間信號(hào)及其振幅譜由s(t)的Fourier變換

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

S(ω)的模|S(ω)

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

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

1.1.2Fourier變換的局限性圖1.2時(shí)變信號(hào)及其振幅譜

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

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

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

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

*本章討論時(shí)-頻分布理論

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

即式中

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

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

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

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

1.2.6解析信號(hào)的快速數(shù)值計(jì)算

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

注：取N為2的整數(shù)次冪，補(bǔ)零個(gè)數(shù)大于原數(shù)據(jù)個(gè)數(shù)。例1.3用MATLAB編程計(jì)算實(shí)信號(hào)的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