現代信號處理方法主講：

王忠仁2012年2月授課內容〇、MATLAB語法要點

一、信號的時頻分析二、Radon-Wigner變換三、分數階Fourier變換四、高階譜估計參考書[1]

張賢達.現代信號處理.清華大學出版社.1995[2]

張賢達,保錚.非平穩信號分析與處理.國防工業出版社.1998[3]

張賢達.時間序列分析-高階統計量方法.清華大學出版社.1996[4]

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

MATLAB語法要點

0．２MATLAB語法要點算術運算符：標量運算:

加減乘除冪+-*/^矩陣運算:

加減乘右除左除冪轉置復共軛轉置

+-*/\^.’’

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

標量情況有:2/5=0.4;2\5=2.5矩陣對應元素運算:乘右除左除冪

.*./.\.^

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

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

語句體

end其中:表達式1為循環初值;

表達式2為步長,省略時默認步長為1;

表達式3為循環終值.例

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

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

end其中:當表達式值為真時執行語句體;

當表達式值為假時終止循環.例:找出階乘超過10100的最小整數.

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

if語句格式:if邏輯表達式1

語句體1

elseif

邏輯表達式2

語句體2

elseif

邏輯表達式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常用函數繪圖函數

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

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

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

引言

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

s(t)的Fourier變換為

(1.1)

1.1.1Fourier變換回顧

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

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

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

S(ω)的Fourier逆變換為

Fourier變換可以用來對信號進行頻譜分析或頻率濾波。(1.2)例1.1用MATLAB編程形成雷克子波,并進行頻譜分析.%峰值頻率為fp的雷克子波的表達式:%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;%數乘向量,pi是MATLAB定義的特殊變量πtt2=tt.*tt;%向量對應元素相乘得新向量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)%傅立葉變換的前一半數據取模

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)所有時間加權累積得到的，

S(ω)的模|S(ω)

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

在現代信號處理中常常要了解時變信號的時頻局部變化特征，Fourier變換無能為力。

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

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

1.1.3聯合時-頻表示時間域頻率域時頻平面圖1.3聯合時-頻表示平面

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

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

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

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

即式中

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

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

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

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

1.2.6解析信號的快速數值計算

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

注：取N為2的整數次冪，補零個數大于原數據個數。例1.3用MATLAB編程計算實信號的Hilbert變換與包絡.%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