一段時間沒寫公衆号，今天正好有個朋友發了一段語音，可以用來做信号分析，故分享一下MATLAB短時傅裡葉變換和小波變換的時頻分析
簡介
本文主要給定一小段音頻，通過短時傅裡葉變換和小波變換制作時頻圖。音頻的采樣率為44100，
短時傅裡葉變換
在matlab中，短時傅裡葉變換的分析函數為spectrogram，其使用情況如下：
功能：使用短時傅裡葉變換得到信号的頻譜圖。
文法：
[S,F,T,P]=spectrogram(x,window,noverlap,nfft,fs)
[S,F,T,P]=spectrogram(x,window,noverlap,F,fs)
說明：當使用時無輸出參數，會自動繪制頻譜圖；有輸出參數，則會傳回輸入信号的短時傅裡葉變換。當然也可以從函數的傳回值S,F,T,P繪制頻譜圖，具體參見例子。
參數：
x---輸入信号的向量。預設情況下，即沒有後續輸入參數，x将被分成8段分别做變換處理，如果x不能被平分成8段，則會做截斷處理。預設情況下，其他參數的預設值為：window---窗函數，預設為nfft長度的海明窗Hamming；noverlap---每一段的重疊樣本數，預設值是在各段之間産生50%的重疊；nfft---做FFT變換的長度，預設為256和大于每段長度的最小2次幂之間的最大值。另外，此參數除了使用一個常量外，還可以指定一個頻率向量F；fs---采樣頻率，預設值歸一化頻率。
Window---窗函數，如果window為一個整數，x将被分成window段，每段使用Hamming窗函數加窗。如果window是一個向量，x将被分成length(window)段，每一段使用window向量指定的窗函數加窗。是以如果想擷取specgram函數的功能，隻需指定一個256長度的Hann窗。
Noverlap---各段之間重疊的采樣點數。它必須為一個小于window或length(window)的整數。其意思為兩個相鄰窗不是尾接着頭的，而是兩個窗有交集，有重疊的部分。
Nfft---計算離散傅裡葉變換的點數。它需要為标量。
Fs---采樣頻率Hz，如果指定為[]，預設為1Hz。
S---輸入信号x的短時傅裡葉變換。它的每一列包含一個短期局部時間的頻率成分估計，時間沿列增加，頻率沿行增加。如果x是長度為Nx的複信号，則S為nfft行k列的複矩陣，其中k取決于window，如果window為一個标量，則k = fix((Nx-noverlap)/(window-noverlap))；如果window為向量，則k = fix((Nx-noverlap)/(length(window)-noverlap))。對于實信号x，如果nfft為偶數，則S的行數為(nfft/2+1)，如果nfft為奇數，則行數為(nfft+1)/2，列數同上。
F---在輸入變量中使用F頻率向量，函數會使用Goertzel方法計算在F指定的頻率處計算頻譜圖。指定的頻率被四舍五入到與信号分辨率相關的最近的DFT容器(bin)中。而在其他的使用nfft文法中，短時傅裡葉變換方法将被使用。對于傳回值中的F向量，為四舍五入的頻率，其長度等于S的行數。
T---頻譜圖計算的時刻點，其長度等于上面定義的k，值為所分各段的中點。
P---能量譜密度PSD(Power Spectral Density)，對于實信号，P是各段PSD的單邊周期估計；對于複信号，當指定F頻率向量時，P為雙邊PSD。P矩陣的元素計算公式如下P(I,j)=k|S(I,j)|2，其中的的k是實值标量，定義如下對于單邊PSD，計算公式如下，其中w(n)表示窗函數，Fs為采樣頻率，在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是對應的能量譜密度
figure
imagesc(T,F,10*log10(abs(P)));
set(gca,'YDir','normal')
colorbar;
xlabel('時間 t/s');
ylabel('頻率 f/Hz');
title('短時傅裡葉時頻圖'); 複制
注意：
- nfft越大，頻域的分辨率就越高（分辨率=fs/nfft），但離瞬時頻率就越遠；
- noverlap影響時間軸的分辨率，越接近nfft，分辨率越高，相應的備援就越多，計算量越大，但計算機隻要能承受，問題不大。
小波變換
首先，在matlab中，小波變換的分析函數為cwt，其使用情況如下：
功能：實作一維連續小波變換的函數。
文法：
COEFS=cwt(S, SCALES, 'wname')
COEFS=cwt(S, SCALES, 'wname', 'plot')
COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE')
COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE', XLIM)
參數：
COEFS=cwt(S, SCALES, 'wname') 采用'wname'小波，在正、實尺度SCALES下計算向量一維小波系數。
COEFS=cwt(S, SCALES, 'wname', 'plot') 除了計算小波系數外，還加以圖形顯示。
COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE') 計算并畫出連續小波變換的系數，并使用PLOTMODE對圖形着色。
COEFS=cwt(S, SCALES, 'wname', 'plot') 相當于 格式 COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE') 中的文法 COEFS=cwt(S, SCALES, 'wname', 'absglb')
COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE', XLIM) 能夠計算并畫出連續小波變換的系數。系數使用PLOTMODE和XLIM進行着色。其中：XLIM=[x1,x2]，并且有如下關系：1<=x1<=x2<=length(S)。
MODE值含義：
'lvl' scale-by-scale着色模式
'glb' 考慮所有尺度的着色模式
'abslvl'或'lvlabs' 使用系數絕對值的scale-by-scale着色模式
'absglb'或'glbabs' 使用系數絕對值并考慮所有尺度的着色模式
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); % 将尺度轉換為頻率 頻率在0-500Hz取1024
coefs = cwt(Au(totalscal,1),scals,wavename); % 求連續小波系數
t=0:1/Fs:(totalscal-1)/Fs;
figure
imagesc(t,f,abs(coefs));
set(gca,'YDir','normal')
colorbar;
xlabel('時間 t/s');
ylabel('頻率 f/Hz');
title('小波時頻圖'); 複制