前言

上一篇博客講了離散傅里葉變換，里面的實例是對整個信號進行計算，雖然理論上有N點傅里葉變換(本博客就不區分FFT和DFT了，因為它倆就是一個東東，只不過復雜度不同)，但是我個人理解是這個N點是信號前面連續的N個數值，即N點FFT意思就是截取前面N個信號進行FFT，這樣就要求我們的前N個采樣點必須包含當前信號的一個周期，不然提取的余弦波參數與正確的疊加波的參數相差很大。

如果在N點FFT的時候，如果這N個采樣點不包含一個周期呢？或者說我們的信號壓根不是一個周期函數咋辦？或者有一段是噪音數據呢？如果用FFT計算，就會對整體結果影響很大，然后就有人想通過局部來逼近整體，跟微積分的思想很像，將信號分成一小段一小段，然后對每一小段做FFT，就跟分段函數似的，無數個分段函數能逼近任意的曲線((⊙o⊙)…應該沒錯吧)，這樣每一段都不會互相影響到了。

下面的參考博客中有一篇的一句話很不錯：在短時傅里葉變換過程中，窗的長度決定頻譜圖的時間分辨率和頻率分辨率，窗長越長，截取的信號越長，信號越長，傅里葉變換后頻率分辨率越高，時間分辨率越差；相反，窗長越短，截取的信號就越短，頻率分辨率越差，時間分辨率越好，也就是說短時傅里葉變換中，時間分辨率和頻率分辨率之間不能兼得，應該根據具體需求進行取舍。

國際慣例，參考博客：

基于MATLAB短時傅里葉變換和小波變換的時頻分析

小波前奏–短時傅里葉變換

短時傅里葉變換原理解

matlab自帶的短時傅里葉變換函數spectrogram

理論及實現

其實就是多了幾個參數，需要指定的有：

• 每個窗口的長度：nsc
• 每相鄰兩個窗口的重疊率：nov
• 每個窗口的FFT采樣點數：nff

可以計算的有：

• 海明窗：w=hamming(nsc, 'periodic')

• 信號被分成了多少片：$\frac{len(S)?nsc}{nsc?nov}$

• 短時傅里葉變換：
  $$X(k)=\sum _{n=1}^{N}w(n)?x(n)?exp(?j?2\pi ?(k?1)?(n?1)/N),1<=k<=N$$
  其實和FFT的公式一樣，只不過多了個海明窗加權

直接擼代碼：

①先設置參數：

%默認設置： % nsc=floor(L/4.5);%海明窗的長度 % nov=floor(nsc/2);%重疊率 % nff=max(256,2^nextpow2(nsc));%N點采樣長度 %也可手動設置 nsc=100;%海明窗的長度,即每個窗口的長度 nov=30;%重疊率 nff=256;%N點采樣長度

這里面有個默認設置，就是調用matlab自帶的短時傅里葉變換時，如果沒指定相關參數，就會采用默認參數值，這個可以去mathwork官網看。

②計算海明窗以及初始化結果值：

h=hamming(nsc, 'periodic');%計算海明窗的數值，給窗口內的信號加權重 coln = 1+fix((L-nsc)/(nsc-nov));%信號被分成了多少個片段 %如果nfft為偶數，則S的行數為(nfft/2+1)，如果nfft為奇數，則行數為(nfft+1)/2 %因為matlab的FFT結果是對稱的，只需要一半 rown=nff/2+1; STFT_X=zeros(rown,coln);%初始化最終結果

這里的信號被劃分的片段數目可以按照卷積的方法計算

③對每個片段碼公式：

%對每個片段進行fft變換 index=1;%當前片段第一個信號位置在原始信號中的索引 for i=1:coln%提取當前片段信號值,并用海明窗進行加權temp_S=S(index:index+nsc-1).*h';%進行N點FFT變換temp_X=fft(temp_S,nff);%取一半STFT_X(:,i)=temp_X(1:rown)';%將索引后移index=index+(nsc-nov); end

可以發現我這里沒碼公式，因為上一篇博客證明了手擼的DFT與matlab自帶的FFT公式一樣，有高度強迫癥的可以把上一篇博客的DFT寫成一個函數，然后把此處的FFT換成你的函數名即可。注意這里的關鍵操作有兩點：

• 對當前窗口的輸入信號進行海明加權
• 窗口中輸入信號的獲取方法有點類似于卷積，卷積核大小是1*nsc，步長是nsc-nov

④正確性驗證：與matlab自帶的STFT函數spectrogram的結果進行比較：

%% matlab自帶函數 [spec_s,spec_f,spec_t]=spectrogram(S,hamming(nsc, 'periodic'),nov,nff,Fs); %減法，看看差距 plot(abs(spec_s)-abs(STFT_X))

啥也不說了，穩如狗

頻譜解讀

直接計算每個坐標軸的數值就知道其代表的意思了，其實它和一款叫做Praat的軟件所畫的圖很類似，貼一張圖，來源百度

上面是聲線，就是直接audioread聲音得到的數值，下面就是聲譜圖，看著是二維圖形，其實是3D圖形，坐標軸分別代表時間，頻率，以及當前時間在當前頻率上的幅值。

那么在matlab中如何計算這些值？如下：

回顧一下整個快速離散傅里葉變換的流程：

• 利用窗口和重疊率對整個輸入信號進行片段劃分
• 對每個片段的信號做N點傅里葉變換，并利用海明窗加權

接下來解析三個坐標，先規定一下橫坐標表示頻率，縱坐標表示時間，第三個坐標表示幅值：

• 第三個坐標：幅值的計算與FFT的幅值計算一樣，都是計算得到的結果除以采樣點N，再乘以2，只不過這里要利用海明值進行縮放處理

  % 海明窗口的均值 K = sum(hamming(nsc, 'periodic'))/nsc; %獲取幅值(除以采樣長度即可，后面再決定那幾個除以2)，并依據窗口的海明系數縮放 STFT_X = abs(STFT_X)/nsc/K; % correction of the DC & Nyquist component if rem(nff, 2) % 如果是奇數，說明第一個是奈奎斯特點，只需對2:end的數據乘以2STFT_X(2:end, :) = STFT_X(2:end, :).*2; else % 如果是偶數，說明首尾是奈奎斯特點，只需對2:end-1的數據乘以2STFT_X(2:end-1, :) = STFT_X(2:end-1, :).*2; end

• 橫坐標：頻率
  首先要知道當前窗口代表了多少頻率，它是在總頻率Fs的基礎上，對每個窗口進行nff點采樣，需要注意的是進行多少次nff采樣，當前窗口就有多少個頻率值，它們之間是均勻的Fs/nff，這個也通常被稱作頻率到分辨率的比率。這里看論文《Constant-Q transform toolbox for music processing 》中的一句話：

  The discrete short-time Fourier transform gives a constant resolution for each bin or frequency sampled equal to the sampling rate dividedbythewindowsizein samples. This means,for example,if we takea window of 1024 samples with a sampling rate of 32O30 samples/s (reasonable for musical signals),there solution is 31.3Hz

  就是說我們從采樣率為32030Hz的樣本中挑選包含1024個樣本的窗口，那么分辨率就是$\frac{32030}{1024}=31.3Hz$

  所以橫坐標的值就是$i\times \frac{Fs}{nff},(i\in (0,nff))$

  %頻率軸 STFT_f=(0:rown-1)*Fs./nff;

• 縱坐標時間：
  因為采用了窗口，所以縱坐標的時間比實際時間短，每個坐標代表當前窗口中間信號在原始信號中的索引，窗口是nsc，重疊率是nov，那么每次向前挪的步長為nsc-nov，總共挪coln-1次，需要注意的是這個挪是在采樣點上挪，需要將采樣點挪轉換為時間挪，由于整個信號采樣率是Fs，代表每秒的采樣數，那么相鄰的采樣點時間間隔是1/Fs，自然就得到了縱坐標的刻度：

  %時間軸 STFT_t=(nsc/2:(nsc-nov):nsc/2+(coln-1)*(nsc-nov))/Fs;

• 額外信息：賦值轉分貝
  我也不清楚這個意義是啥，但是在matlab中有對應函數，而且軟件praat和默認函數spectrogram的結果圖中都有這個信息，所以還是算一下吧，公式很簡單$y=20?\log 10(x)$，直接碼公式:

  STFT_X = 20*log10(STFT_X + 1e-6);

  這里加個很小的值以后，畫圖好看點

坐標值都得到，直接mesh出來就行，整個代碼如下：

%% 畫頻譜圖 % 海明窗口的均值 K = sum(hamming(nsc, 'periodic'))/nsc; %獲取幅值(除以采樣長度即可，后面再決定哪幾個除以2)，并依據窗口的海明系數縮放 STFT_X = abs(STFT_X)/nsc/K; % correction of the DC & Nyquist component if rem(nff, 2) % 如果是奇數，說明第一個是奈奎斯特點，只需對2:end的數據乘以2STFT_X(2:end, :) = STFT_X(2:end, :).*2; else % 如果是偶數，說明首尾是奈奎斯特點，只需對2:end-1的數據乘以2STFT_X(2:end-1, :) = STFT_X(2:end-1, :).*2; end% convert amplitude spectrum to dB (min = -120 dB) %將幅值轉換成分貝有函數是ydb=mag2db(y)，這里直接算 STFT_X = 20*log10(STFT_X + 1e-6); %時間軸 STFT_t=(nsc/2:(nsc-nov):nsc/2+(coln-1)*(nsc-nov))/Fs; %頻率軸 STFT_f=(0:rown-1)*Fs./nff; % plot the spectrogram figure surf(STFT_f, STFT_t, STFT_X') shading interp axis tight box on view(0, 90) set(gca, 'FontName', 'Times New Roman', 'FontSize', 14) xlabel('Frequency, Hz') ylabel('Time, s') % title('Amplitude spectrogram of the signal') title('手繪圖')handl = colorbar; set(handl, 'FontName', 'Times New Roman', 'FontSize', 14) ylabel(handl, 'Magnitude, dB')

對比一下matlab自帶函數spectrogram和我們手繪圖的正面圖和3D圖：

從圖里面很容易看出來咱們輸入的這個音頻信號由50和100兩個頻率組成，幅值大概在10到20，what？咋少了一半，(⊙o⊙)…，后面再研究為啥少了一半還，按理說乘以2了呀，反正頻率對了

后記

感覺對于FFT的理解告一段落，先把蝶形算法擱著，下一步就是折騰常Q變換(Constant-Q transform)了，目前的用處是一個音樂的一拍可能有很多音組合而成，但是每個音的頻率又不一樣，那么就需要設置不同的窗口進行采樣，相當于進行了多次STFT操作，只不過每次的窗口大小不同罷了，有興趣可以看一波論文：《Calculation of a constant Q spectral transform 》，有張圖介紹了CQT和DFT的區別，具體我還在研究，感覺就是把STFT的for循環里面的N變成動態計算的就行了。

貼一波代碼，直接貼了，懶得貼網盤：

%短時傅里葉變換STFT %依據FFT手動實現STFT clear clc close all Fs = 1000; % Sampling frequency T = 1/Fs; % Sampling period L = 1000; % Length of signal t = (0:L-1)*T; % Time vector S = 20*cos(100*2*pi*t)+40*cos(50*2*pi*t);%0.2-0.7*cos(2*pi*50*t+20/180*pi) + 0.2*cos(2*pi*100*t+70/180*pi) ;%% 所需參數 %主要包含:信號分割長度(默認分割8個窗口)，海明窗口，重疊率，N點采樣 %默認設置： % nsc=floor(L/4.5);%海明窗的長度 % nov=floor(nsc/2);%重疊率 % nff=max(256,2^nextpow2(nsc));%N點采樣長度 %也可手動設置 nsc=100;%海明窗的長度,即每個窗口的長度 nov=0;%重疊率 nff=max(256,2^nextpow2(nsc));%N點采樣長度 %% 手動實現STFT h=hamming(nsc, 'periodic');%計算海明窗的數值，給窗口內的信號加權重 coln = 1+fix((L-nsc)/(nsc-nov));%信號被分成了多少個片段 %如果nfft為偶數，則S的行數為(nfft/2+1)，如果nfft為奇數，則行數為(nfft+1)/2 %因為matlab的FFT結果是對稱的，只需要一半 rown=nff/2+1; STFT_X=zeros(rown,coln);%初始化最終結果 %對每個片段進行fft變換 index=1;%當前片段第一個信號位置在原始信號中的索引 for i=1:coln%提取當前片段信號值,并用海明窗進行加權temp_S=S(index:index+nsc-1).*h';%進行N點FFT變換temp_X=fft(temp_S,nff);%取一半STFT_X(:,i)=temp_X(1:rown)';%將索引后移index=index+(nsc-nov); end%% matlab自帶函數 spectrogram(S,hamming(nsc, 'periodic'),nov,nff,Fs); title('spectrogram函數畫圖') [spec_X,spec_f,spec_t]=spectrogram(S,hamming(nsc, 'periodic'),nov,nff,Fs); %減法，看看差距 % plot(abs(spec_X)-abs(STFT_X))%% 畫頻譜圖 % 海明窗口的均值 K = sum(hamming(nsc, 'periodic'))/nsc; %獲取幅值(除以采樣長度即可，后面再決定哪幾個乘以2)，并依據窗口的海明系數縮放 STFT_X = abs(STFT_X)/nsc/K; % correction of the DC & Nyquist component if rem(nff, 2) % 如果是奇數，說明第一個是奈奎斯特點，只需對2:end的數據乘以2STFT_X(2:end, :) = STFT_X(2:end, :).*2; else % 如果是偶數，說明首尾是奈奎斯特點，只需對2:end-1的數據乘以2STFT_X(2:end-1, :) = STFT_X(2:end-1, :).*2; end% convert amplitude spectrum to dB (min = -120 dB) %將幅值轉換成分貝有函數是ydb=mag2db(y)，這里直接算 STFT_X = 20*log10(STFT_X + 1e-6); %時間軸 STFT_t=(nsc/2:(nsc-nov):nsc/2+(coln-1)*(nsc-nov))/Fs; %頻率軸 STFT_f=(0:rown-1)*Fs./nff; % plot the spectrogram figure surf(STFT_f, STFT_t, STFT_X') shading interp axis tight box on view(0, 90) set(gca, 'FontName', 'Times New Roman', 'FontSize', 14) xlabel('Frequency, Hz') ylabel('Time, s') title('Amplitude spectrogram of the signal')handl = colorbar; set(handl, 'FontName', 'Times New Roman', 'FontSize', 14) ylabel(handl, 'Magnitude, dB')

博客已同步更新到微信公眾號中，有問題可以在微信公眾號中私聊，或者在csdn博客下面評論。

總結

以上是生活随笔為你收集整理的【音频处理】短时傅里叶变换的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。