參考文章：

形象易懂講解算法I——小波變換

信號分析的主要目的是尋找一種簡單有效的信號變換方法，使信號包含的重要特征能顯示出來。
在小波變換興起之前，Fourier級數(shù)展開和Fourier分析是刻畫函數(shù)空間、求解微分方程、進行數(shù)學(xué)計算的有效方法和有效的數(shù)學(xué)工具。
從物理直觀上看，一個周期性振動的量可以看成是具有簡單頻率的簡諧振動的疊加。Fourier級數(shù)展開則是這一物理過程的數(shù)學(xué)描述。

一.Fourier變換

Fourier變換基本原理

傅里葉變換把無限長的三角函數(shù)作為基函數(shù)：

這個基函數(shù)會伸縮、會平移（其實本質(zhì)并非平移，而是兩個正交基的分解）。
縮得窄，對應(yīng)高頻；伸得寬，對應(yīng)低頻。然后這個基函數(shù)不斷和信號做相乘。某一個尺度（寬窄）下乘出來的結(jié)果，就可以理解成信號所包含的當(dāng)前尺度對應(yīng)頻率成分有多少。于是，基函數(shù)會在某些尺度下，與信號相乘得到一個很大的值，因為此時二者有一種重合關(guān)系。那么我們就知道信號包含該頻率的成分的多少。

仔細體會可以發(fā)現(xiàn)，這一步其實是在計算信號和三角函數(shù)的相關(guān)性。


這兩種尺度能乘出一個大的值（相關(guān)度高），所以信號包含較多的這兩個頻率成分，在頻譜上這兩個頻率會出現(xiàn)兩個峰。
以上，就是粗淺意義上傅里葉變換的原理。

Fourier變換實例

舉例如下：

fun = @(t) cos(2*pi*10*t)+cos(2*pi*25*t)+cos(2*pi*50*t)+cos(2*pi*100*t); x = 0:0.001:1; y = fun(x); Dfy = fft(y); % 離散Fourier變換 Dfy_shift = fftshift(Dfy); % 對稱變換得到對稱的Fourier頻譜 figure plot(y); axis([0,1000,-3,4]),xlabel('x'),ylabel('y'); % 原始平穩(wěn)信號figure plot(abs(Dfy)),axis([0,250,0,600]),xlabel('n'),ylabel('|Dfy|'); % |Dft|為離散頻譜幅值信息

相應(yīng)具體圖像如下：

做完FFT（快速傅里葉變換）后，可以在頻譜上看到清晰的四條線，信號包含四個頻率成分。

Fourier變換對非平穩(wěn)過程，具有局限性。
下面簡單介紹：信號雖具有相同頻譜，但具體序列并不相同。

如上圖，最上邊的是頻率始終不變的平穩(wěn)信號。而下邊兩個則是頻率隨著時間改變的非平穩(wěn)信號，它們同樣包含和最上信號相同頻率的四個成分。

做FFT后，我們發(fā)現(xiàn)這三個時域上有巨大差異的信號，頻譜（幅值譜）卻非常一致。尤其是下邊兩個非平穩(wěn)信號，我們從頻譜上無法區(qū)分它們，因為它們包含的四個頻率的信號的成分確實是一樣的，只是出現(xiàn)的先后順序不同。

可見，傅里葉變換處理非平穩(wěn)信號有天生缺陷。它只能獲取一段信號總體上包含哪些頻率的成分，但是對各成分出現(xiàn)的時刻并無所知。因此時域相差很大的兩個信號，可能頻譜圖一樣。

Fourier變換是平穩(wěn)信號分析的最重要的工具。然而在實際運用中，所遇到的信號大多數(shù)并不平穩(wěn)，而是時變頻率信號，這時人們需要知道信號在突變時刻所對應(yīng)的頻率成分。

二.短時傅里葉變換（Short-time Fourier Transform, STFT）

1946年諾貝爾獎獲得者D.Gabor引入了短時傅里葉變換(Short-time Fourier Transform)，雖然短時傅里葉變換隨著窗口在時間軸t上滑動可以摳取頻譜上的所有信息，但從本質(zhì)上講，短時傅里葉變換是一種單一分辨率的信號分析方法，因為它使用一個固定的短時窗函數(shù)（常用Gauss函數(shù)），窗口大小缺乏自適應(yīng)功能，窗口位置可以移動，但是形狀不能改變。

通俗理解，“把整個時域過程分解成無數(shù)個等長的小過程，每個小過程近似平穩(wěn)，再進行傅里葉變換，就知道在哪個時間點上出現(xiàn)了什么頻率。”這就是短時傅里葉變換。

使用STFT存在一個問題，我們應(yīng)該用多寬的窗函數(shù)？窗太寬太窄都有問題：

窗太窄，窗內(nèi)的信號太短，會導(dǎo)致頻率分析不夠精準，頻率分辨率差。

窗太寬，時域上又不夠精細，時間分辨率低。

三.小波分析（Wavelet Analysis）

從事石油信號處理的法國工程師J.Morlet在研究地下巖石油層分布時，于1974年首先提出了小波變換的概念，并且J.Morlet構(gòu)建了光滑柔性的、時頻域局部性能都較好的Morlet小波函數(shù)。數(shù)學(xué)家Meyer，Mallat，Daubechies，K.Chui 等人的工作為小波分析的誕生和發(fā)展奠定了基礎(chǔ)。

小波變換簡介

小波變換做出的改變主要是將無限長的三角函數(shù)基換成了有限長的會衰減的小波基。

從公式可以看出，不同于傅里葉變換，變量只有頻率ω，小波變換有兩個變量：尺度a（scale）和平移量 τ（translation）。尺度a控制小波函數(shù)的伸縮，對應(yīng)于頻率（反比）；平移量 τ控制小波函數(shù)的平移，對應(yīng)于時間。

當(dāng)伸縮、平移到這么一種重合情況時，也會相乘得到一個大的值。這時候和傅里葉變換不同的是，這不僅可以知道信號有這樣頻率的成分，而且知道它在時域上存在的具體位置。

而當(dāng)我們在每個尺度下都平移著和信號乘過一遍后，我們就知道信號在每個位置都包含哪些頻率成分。

小波變換基本原理

兩類：
1.連續(xù)小波變換CWT
2.離散小波變換DWT

區(qū)別： DWT是數(shù)據(jù)的緊湊表示，對于降噪和數(shù)據(jù)壓縮特別有用，而CWT換對于特征提取更好。

CWT函數(shù)

matlab自帶的有兩種實現(xiàn)方式，一種是2006年版本推出的函數(shù)cwt，一種是2016年版本推出的函數(shù)cwt。兩個函數(shù)名稱相同，用法不同。

【舊版本】cwt用法為：

coefs = cwt(x,scales,'wname')

用于獲取實的或者復(fù)的連續(xù)小波變換系數(shù)，其中“wname”為選取的小波，可以用waveinfo查看小波的相關(guān)信息。
舉例：

%針對同一個分段函數(shù)，采用“Haar“小波分解，得到小波的頻譜 fun = @(x)0.*(x<0 |x>5)+x.*(x>=0 & x<1)+1.*(x>=1 & x<4) +(5-x).*(x>=4 & x<5);%分段函數(shù) x=-5:0.2:10;%x自變量 y = fun(x); figure(1) Cab = cwt(y,1:2:32,'haar','3Dplot');%連續(xù)小波變換，并進行頻譜幅值的三維顯示 figure(2) Cab = cwt(y,1:2:32,'haar','plot');%連續(xù)小波變換，并進行頻譜幅值的圖像顯示 %這段程序也可以按照Gabor變換的方式顯示，即 a = 1:2:32; b = 1:length(y); Cab = cwt(y,1:2:32,'haar');%連續(xù)小波變換，并進行頻譜幅值的三維顯示

【新版本】cwt用法為：

[wt,f] = cwt(x,wname,fs)

最明顯的區(qū)別在于，新版本取消了自定義尺度函數(shù)scales的功能。同時新版本還更新替換了一部分小波函數(shù)。
舊版本支持 ‘haar’，‘db’，‘sym’，‘cmor’，‘mexh’，‘gaus’，‘bior’等小波，新版本支持’morse’， 'amor’和 'bump’小波。

常見小波基函數(shù)

目前我們主要是通過用小波分析方法處理信號的結(jié)果與理論結(jié)果的誤差來判定小波基的好壞，由此決定小波基。
常用小波基有Haar小波、Daubechies(dbN)小波、Mexican Hat(mexh)小波、Morlet小波、Meyer小波等。
1.Haar小波
Haar函數(shù)是小波分析中最早用到的一個具有緊支撐的正交小波函數(shù)，也是最
簡單的一個小波函數(shù)，它是支撐域在[0,1]∈t范圍內(nèi)的單個矩形波。
Haar小波的時域和頻域波形如下：

[phi,g1,xval]=wavefun('haar',20);subplot(2,1,1); plot(xval,g1,'LineWidth',2); xlabel('t') title('haar 時域'); g2=fft(g1); g3=abs(g2);subplot(2,1,2); plot(g3,'LineWidth',2); xlabel('f') title('haar 頻域')

2.Morlet小波
連續(xù)復(fù)小波變換
Morlet小波的時域和頻域波形圖：

d=-6; h=6; n=100; [g1,x]=morlet(d,h,n);subplot(2,1,1); plot(x,g1,'LineWidth',2); xlabel('t'); title('morlet 時域'); g2=fft(g1); g3=(abs(g2));subplot(2,1,2); plot(g3,'LineWidth',2); xlabel('f'); title('morlet 頻域');

數(shù)據(jù)處理

1.小波系數(shù)實部圖

2.小波模部等值線圖

3.小波方差圖
小波方差隨尺度a的變化過程，稱為小波方差圖，它能反映信號波動的能量隨尺度a的分布。因此，小波方差圖可用來確定信號中不同種尺度擾動的相對強度和存在的主要時間尺度，即主周期。

4.主周期趨勢圖

根據(jù)小波方差檢驗的結(jié)果，繪制演變的第一和第二主周期小波系數(shù)圖。

具體案例如下：

%進行連續(xù)小波變換得到小波系數(shù)矩陣,選擇復(fù)morlet小波函數(shù)Cab = cwt(Data,1:5:40,'cmor1-1','plot');% 求得系數(shù)的實部realPart=real(Cab);figure(2)contourf(realPart,10,'-');colormap('HSV');colorbar;xlabel('年份');ylabel('時間尺度/年');set(gca,'XTick',[1:10*12:60*12],'XTickLabel',{'1961','1971','1981','1991','2001','2011'}) %更新XTickLabel% 小波方差是模的平方的算數(shù)平均Module=abs(Cab);Module2=Module.^2;fangcha=mean(Module2,1);figure(3);plot(fangcha/12,'k-','linewidth',1);xlabel('時間尺度/年');ylabel('小波方差');set(gca,'XTick',[0:10*12:60*12],'XTickLabel',{'0','10','20','30','40','50','60'}) %更新XTickLabel

得到的相應(yīng)圖像如下所示：

相應(yīng)小波方差圖如下圖所示：

總結(jié)

以上是生活随笔為你收集整理的MATLAB利用小波分析提取周期的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

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