談?wù)凪atlab短時(shí)傅里葉變換和小波變換的時(shí)頻

簡介

本文主要給定一小段音頻，通過短時(shí)傅里葉變換和小波變換制作時(shí)頻圖。音頻的采樣率為44100，

短時(shí)傅里葉變換

在matlab中，短時(shí)傅里葉變換的分析函數(shù)為spectrogram，其使用情況如下：

功能 ：使用短時(shí)傅里葉變換得到信號的頻譜圖。

語法 ：

[S,F,T,P]=spectrogram(x,window,noverlap,nfft,fs)

[S,F,T,P]=spectrogram(x,window,noverlap,F,fs)

說明：當(dāng)使用時(shí)無輸出參數(shù)，會自動繪制頻譜圖；有輸出參數(shù)，則會返回輸入信號的短時(shí)傅里葉變換。當(dāng)然也可以從函數(shù)的返回值S,F,T,P繪制頻譜圖，具體參見例子。

參數(shù) ：

x---輸入信號的向量。默認(rèn)情況下，即沒有后續(xù)輸入?yún)?shù)，x將被分成8段分別做變換處理，如果x不能被平分成8段，則會做截?cái)嗵幚?。默認(rèn)情況下，其他參數(shù)的默認(rèn)值為：window---窗函數(shù)，默認(rèn)為nfft長度的海明窗Hamming；noverlap---每一段的重疊樣本數(shù)，默認(rèn)值是在各段之間產(chǎn)生50%的重疊；nfft---做FFT變換的長度，默認(rèn)為256和大于每段長度的最小2次冪之間的最大值。另外，此參數(shù)除了使用一個(gè)常量外，還可以指定一個(gè)頻率向量F；fs---采樣頻率，默認(rèn)值歸一化頻率。

Window---窗函數(shù)，如果window為一個(gè)整數(shù)，x將被分成window段，每段使用Hamming窗函數(shù)加窗。如果window是一個(gè)向量，x將被分成length(window)段，每一段使用window向量指定的窗函數(shù)加窗。所以如果想獲取specgram函數(shù)的功能，只需指定一個(gè)256長度的Hann窗。

Noverlap---各段之間重疊的采樣點(diǎn)數(shù)。它必須為一個(gè)小于window或length(window)的整數(shù)。其意思為兩個(gè)相鄰窗不是尾接著頭的，而是兩個(gè)窗有交集，有重疊的部分。

Nfft---計(jì)算離散傅里葉變換的點(diǎn)數(shù)。它需要為標(biāo)量。

Fs---采樣頻率Hz，如果指定為[]，默認(rèn)為1Hz。

S---輸入信號x的短時(shí)傅里葉變換。它的每一列包含一個(gè)短期局部時(shí)間的頻率成分估計(jì)，時(shí)間沿列增加，頻率沿行增加。如果x是長度為Nx的復(fù)信號，則S為nfft行k列的復(fù)矩陣，其中k取決于window，如果window為一個(gè)標(biāo)量，則k = fix((Nx-noverlap)/(window-noverlap))；如果window為向量，則k = fix((Nx-noverlap)/(length(window)-noverlap))。對于實(shí)信號x，如果nfft為偶數(shù)，則S的行數(shù)為(nfft/2+1)，如果nfft為奇數(shù)，則行數(shù)為(nfft+1)/2，列數(shù)同上。

F---在輸入變量中使用F頻率向量，函數(shù)會使用Goertzel方法計(jì)算在F指定的頻率處計(jì)算頻譜圖。指定的頻率被四舍五入到與信號分辨率相關(guān)的最近的DFT容器(bin)中。而在其他的使用nfft語法中，短時(shí)傅里葉變換方法將被使用。對于返回值中的F向量，為四舍五入的頻率，其長度等于S的行數(shù)。

T---頻譜圖計(jì)算的時(shí)刻點(diǎn)，其長度等于上面定義的k，值為所分各段的中點(diǎn)。

P---能量譜密度PSD(Power Spectral Density)，對于實(shí)信號，P是各段PSD的單邊周期估計(jì)；對于復(fù)信號，當(dāng)指定F頻率向量時(shí)，P為雙邊PSD。P矩陣的元素計(jì)算公式如下P(I,j)=k|S(I,j)|2，其中的的k是實(shí)值標(biāo)量，定義如下對于單邊PSD，計(jì)算公式如下，其中w(n)表示窗函數(shù)，F(xiàn)s為采樣頻率，在0頻率和奈奎斯特頻率處，分子上的因子2改為1；

MATLAB程序：

复制[Au, Fs]=audioread('audio.mp3');   % Fs 采樣率 44100
[B, F, T, P] = spectrogram(Au(:,1),1024,512,1024,Fs);   % B是F大小行T大小列的頻率峰值，P是對應(yīng)的能量譜密度
figure
imagesc(T,F,10*log10(abs(P)));
set(gca,'YDir','normal')
colorbar;
xlabel('時(shí)間 t/s');
ylabel('頻率 f/Hz');
title('短時(shí)傅里葉時(shí)頻圖');

注意：

• nfft越大，頻域的分辨率就越高（分辨率=fs/nfft），但離瞬時(shí)頻率就越遠(yuǎn)；
• noverlap影響時(shí)間軸的分辨率，越接近nfft，分辨率越高，相應(yīng)的冗余就越多，計(jì)算量越大，但計(jì)算機(jī)只要能承受，問題不大。

小波變換

首先，在matlab中，小波變換的分析函數(shù)為cwt，其使用情況如下：

功能 ：實(shí)現(xiàn)一維連續(xù)小波變換的函數(shù)。

語法 ：

COEFS=cwt(S, SCALES, 'wname')

COEFS=cwt(S, SCALES, 'wname', 'plot')

COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE')

COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE', XLIM)

參數(shù) ：

COEFS=cwt(S, SCALES, 'wname') 采用'wname'小波，在正、實(shí)尺度SCALES下計(jì)算向量一維小波系數(shù)。

COEFS=cwt(S, SCALES, 'wname', 'plot') 除了計(jì)算小波系數(shù)外，還加以圖形顯示。

COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE') 計(jì)算并畫出連續(xù)小波變換的系數(shù)，并使用PLOTMODE對圖形著色。

COEFS=cwt(S, SCALES, 'wname', 'plot') 相當(dāng)于 格式 COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE') 中的語法 COEFS=cwt(S, SCALES, 'wname', 'absglb')

COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE', XLIM) 能夠計(jì)算并畫出連續(xù)小波變換的系數(shù)。系數(shù)使用PLOTMODE和XLIM進(jìn)行著色。其中：XLIM=[x1,x2]，并且有如下關(guān)系：1<=x1<=x2<=length(S)。

MODE值含義：

'lvl' scale-by-scale著色模式

'glb' 考慮所有尺度的著色模式

'abslvl'或'lvlabs' 使用系數(shù)絕對值的scale-by-scale著色模式

'absglb'或'glbabs' 使用系數(shù)絕對值并考慮所有尺度的著色模式

COEFS行的大小等于SCALES尺度的長度，COEFS列的大小等于信號S的長度。

MATLAB程序：

复制totalscal=1024*16;
wavename='cmor3-3';
Fc=centfrq(wavename); % 小波的中心頻率
c=2*Fc*totalscal;    
scals=c./(1:totalscal);
f=scal2frq(scals,wavename,1/Fs); % 將尺度轉(zhuǎn)換為頻率   頻率在0-500Hz取1024
coefs = cwt(Au(totalscal,1),scals,wavename); % 求連續(xù)小波系數(shù)
t=0:1/Fs:(totalscal-1)/Fs;
figure
imagesc(t,f,abs(coefs));
set(gca,'YDir','normal')
colorbar;
xlabel('時(shí)間 t/s');
ylabel('頻率 f/Hz');
title('小波時(shí)頻圖');