在語音與音樂處理過程中，常用到短時傅里葉變換（Short Time Fourier Transformation, STFT）。在一些學習路徑中，STFT也是學習小波之前的預備知識。本文簡單實現了 Matlab 中 Spectrogram 函數，在沒有小波知識支撐下，討論了參數的選擇，以及分辨率相關的問題。
參考博客：
來源：CSDN 作者：風翼冰舟
來源：知乎 作者：咚懂咚懂咚
來源：CSDN 作者：水可馬二
來源：CSDN 作者：沈子恆
來源：CSDN 作者：gent__chen
短時傅里葉變換簡介
自己寫一下對STFT的認識，不一定對。
短時傅里葉變換，其實還是傅里葉變換，只不過把一段長信號按信號長度（nsc）、重疊點數（nov）重新采樣。原始信號的每一個數據點，都有一個電壓值對應，即只有一個直流自由度。重新采樣之后，數據點個數變少，每一個數據點由（nsc）個點組成，即得到了其他頻率的自由度。但是這個點的頻譜是有上下限的——上限受采樣頻率界定，下限受數據時長界定，和 FFT 一樣。
這么處理數據是有緣由的。語音與音樂信號中，信號頻率常常變化，而且頻率成分豐富，導致簡單的傅里葉變換不能很好的描述信號。然而我們人耳處理聲音的能力很強，可以對很短的一段聲音進行精確的頻率分析。所以在語音識別以及語音信號預處理過程中，STFT 是 FFT 的仿生改良版（個人理解）。當然，在其他聲學相關方向中，STFT也是蠻有用的。沒查過相關文獻不敢亂說，模態分析、瞬態特性等應該能用上（個人猜測）。
上圖是我用本文的例程實現的語音“XX大學”的 STFT 結果。語音是自己拿錄音筆錄的，時間域上沒有做截斷，但由於自己聲線太低，頻率域截到 2000 Hz。可以看到，我在開始之后等了一段時間才開始說話，后面有幾段“嘎嘣脆”的噪聲。還可以看到，我的元音發音還是很清晰的，但聲線真心低。理論上能從這張圖還原出原始信號，不在話下。
Spectrogram 函數函數實現
簡單編了一個 STFT 函數如下：
function [ S , W , T ] = mf_spectrogram...
( signal , nsc , nov , fs )
%MF_SPECTROGRAM Short-time Fourier transform realization
% Input
% signal - Signal vector
% nsc - Abb. of Number SeCtion
% nov - Abb. of Number OverLap
% fs - Abb. of Frequency of Sample
% Output
% S - A matrix that each colum is a FFT for time of nsc
% W - A vector labeling frequency
% T - A vector labeling time
% Signal Preprocessing
h = hamming(nsc, 'periodic'); % Hamming weight function
L = length(signal); % Length of Signal
nst = nsc-nov; % Number of STep per colum
ncl = fix( (L-nsc)/nst ) + 1; % Number of CoLum
nff = 2^nextpow2(nsc); % Number of Fast Fourier Transformation
Ps = zeros(nsc,ncl);
for n = 1:ncl % Ps means Processed Signal
Ps(:,n) = signal( (n-1)*nst+1 : (n-1)*nst+nsc ).*h';
end % Ps is a matrix
% Short-time Fourier transform
STFT0 = fft(Ps,nff);
% Turn into DFT in dB
STFT1 = abs(STFT0/nff);
STFT2 = STFT1(1:nff/2+1,:); % Symmetric results of FFT
STFT2(2:end-1,:) = 2*STFT2(2:end-1,:); % Add the value of the other half
STFT3 = 20*log10(STFT2); % Turn sound pressure into dB level
% Axis Generating
fax = fs*(0:(nff/2))/nff; % Frequency axis setting
tax = ( .5*nsc : nst : nst*(ncl-1)+.5*nsc ) / fs; % Time axis generating
[ffax,ttax] = meshgrid(tax,fax); % Generating grid of figure
% Output
W = ffax;
T = ttax;
S = STFT3;
end
為了節省代碼量，我搞了一個繪圖函數：
function [ ] = my_pcolor( f , t , s )
%MY_PCOLOR 繪圖函數
% Input
% f - 頻率軸矩陣
% t - 時間軸矩陣
% s - 幅度軸矩陣
gca = pcolor(f,t,s); % 繪制偽彩色圖
set(gca, 'LineStyle','none'); % 取消網格，避免一片黑
handl = colorbar; % 彩圖坐標尺
set(handl, 'FontName', 'Times New Roman', 'FontSize', 14)
ylabel(handl, 'Magnitude, dB') % 坐標尺名稱
end
下面是實現的例程：
clc
clear
close all
% 基本參數
fa = [ 50 800 ]; % 掃描頻率上下限
fs = 6400; % 采樣頻率
% 分辨率相關設定參數
te = 1; % [s] 掃頻時間長度
fre = 8; % [s] 頻率分辨率
tre = 0.002; % [Hz] 時間分辨率
% Chirp 信號生成
t = 0:1/fs:te; % [s] 時間序列
sc = chirp(t,fa(1),te,fa(2)); % 信號生成
% 分辨率相關輸入參數
nsc = floor(fs/fre);
% nov = floor(nsc-(fs*tre));
nov = floor(nsc*0.9);
nff = max(256,2^nextpow2(nsc));
% 計算與繪制結果
subplot(1,3,1) % 繪制自編函數結果
[S,W,T] = mf_spectrogram(sc,nsc,nov,fs);
my_pcolor(W,T,S)
caxis([-130.86,-13.667]);
title('自編函數');
xlabel('時間 second');
ylabel('頻率 Hz');
subplot(1,3,2) % 繪制 Matlab 函數結果
s = spectrogram(sc,hamming(nsc),nov,nff,fs,'yaxis');
% Turn into DFT in dB
s1 = abs(s/nff);
s2 = s1(1:nff/2+1,:); % Symmetric results of FFT
s2(2:end-1,:) = 2*s2(2:end-1,:); % Add the value of the other half
s3 = 20*log10(s2); % Turn sound pressure into dB level
my_pcolor(W,T,s3 )
caxis([-130.86,-13.667]);
title('Matlab 自帶函數');
xlabel('時間 second');
ylabel('頻率 Hz');
subplot(1,3,3) % 兩者誤差
my_pcolor(W,T,20*log10(abs(10.^(s3/20)-10.^(S/20))))
caxis([-180,-13.667]);
title('error');
ylabel('頻率 Hz');
xlabel('時間 second');
suptitle('Spectrogram 實現與比較');
跑的結果我就不貼了，Demo確定是沒問題的。網上也有很多相關話題的，例程都比較簡單，但我非常善於把問題復雜化：te（掃描長度）、fre（頻率下限）、tre（時域分辨率）、nsc（單段數據長度）、nov（重疊數據點數）、nff（FFT點數） 等參數都是為了看設定什么樣的參數能夠使得 STFT 結果最優。我只是科普性地了解過小波，公式推導還沒展開，所以在此只能“實驗性”地討論“尺度”相關的話題。
觀察頻率上限受到采樣頻率限制，頻率下限受到nsc限制。nsc越大，單段數據時間跨度越長，在該段時間內頻率如果變化快，會導致頻域分辨率降低（te = 1 ; fre = 2 ; nov = 99%）。另一方面，為了保持 FFT 頻率分辨率，設置 nff 不低於256，這使得當nsc較小時，單段信號中，低頻成分可能也就幾個周期就結束了，等同於時域加了一個矩形窗，最終造成了頻域產生旁瓣。一句話，單段信號過短，低頻效果不好；單段信號過長，捕捉不到頻率快速變化的信號。
數據長度、FFT 點數對結果的影響
實驗性地比較了一下不同數據長度、FFT點數對結果的影響：
% This script is demonstrating zero effects
% Basic parameter
fs = 100; % 采樣頻率
Ndata = [ 30 60 500 ]; % 數據長度
Nff = [ 32 64 512 ]; % FFT的數據長度
% Signal Generating
t1 = ( 0:Ndata(1)-1 )/fs;
t2 = ( 0:Ndata(2)-1 )/fs;
t3 = ( 0:Ndata(3)-1 )/fs;
x1 = 0.5*sin(2*pi*15*t1)+2*sin(2*pi*40*t1); % 時間域信號
x2 = 0.5*sin(2*pi*15*t2)+2*sin(2*pi*40*t2); % 時間域信號
x3 = 0.5*sin(2*pi*15*t3)+2*sin(2*pi*40*t3); % 時間域信號
% FFT and Plot
for n = 1:3 % Ndata sweep
xname = ['x',num2str(n)];
x = eval(xname);
for m = 1:3 % Nff sweep
name = ['y',num2str(n),num2str(m)];
eval([name '=fft(x,Nff(m));'])
y = eval(name);
Y = abs(y);
f = (0:Nff(m)-1)*fs/Nff(m); % 真實頻率
subplot(3,3,(n-1)*3+m)
plot(f(1:Nff(m)/2),[Y(1),Y(2:Nff(m)/2)*2]/min(Ndata(n),Nff(m)));
xlabel('頻率/Hz');
ylabel('振幅');
ylim([0,2]);
title(['Ndata=',num2str(Ndata(n)),' Nfft=',num2str(Nff(m))]);
grid on;
end
end
為了節省代碼，同時運行結束后保留計算數據，使用了 eval 函數對變量和字符串進行相互轉化。繪制了9個圖如下：
橫向是提高 FFT 點數的比較，縱向是提高數據長度的比較。可以看到，提高 FFT 點數會使得頻域分辨率提高；增長數據長度，能夠減少旁瓣的生成。
END