短時傅里葉變換

設時域信號為x(l)，分幀加窗處理后得到的第n幀信號為x_n(m)，則x_n(m)滿足下式：

其中N是每一幀信號的長度，n是幀序號，m是一幀中資料的序號，

時域信號x(l)的離散短時傅里葉變換為：

其中k是譜線號，

當N是2的整數倍時，這個離散短時傅里葉變換可以使用FFT來計算，

MATLAB程式

MATLAB程式演示信號分幀、加窗、求離散短時傅里葉變換，并最終使用三維圖展示結果，

其中打開的test.wav檔案是一個8kHz采樣率的音頻檔案，這個檔案可以自由錄制得到，

clc
clear
close all

[x,fs]=audioread('test.wav');% 讀入資料檔案
t=(1:length(x))/fs;

wlen=256;           % 幀長
inc=128;            % 幀移
win=hanning(wlen);  % 窗函式
nfft=wlen;          % 短時DFT的點數

y=STFFT(x,win,nfft,inc);% 求短時傅里葉變換
fn=size(y,2);% 幀數 = y的列數

freq=(0:wlen/2)*fs/wlen;% 計算FFT后的頻率刻度
frameTime=(((1:fn)-1)*inc+wlen/2)/fs;% 分幀后計算每幀對應的時間


figure;
plot(t,x);
xlabel('時間/s');ylabel('幅度');
title('原始時域信號');

figure;
surf(frameTime,freq,abs(y));
colormap(jet);
xlabel('時間/s');ylabel('頻率/Hz');
title('時頻三維圖');

function d=STFFT(x,win,nfft,inc)
xn=enframe(x,win,inc)';
y=fft(xn,nfft);
d=y(1:(1+nfft/2),:);
end
 
function frameout=enframe(x,win,inc)
 
nx=length(x(:));            % 取資料長度
nwin=length(win);           % 取窗長
if (nwin == 1)              % 判斷窗長是否為1，若為1，即表示沒有設窗函式
   len = win;               % 是，幀長=win
else
   len = nwin;              % 否，幀長=窗長
end
if (nargin < 3)             % 如果只有兩個引數，設幀inc=幀長
   inc = len;
end
nf = fix((nx-len+inc)/inc); % 計算幀數
frameout=zeros(nf,len);     % 初始化
indf= inc*(0:(nf-1)).';     % 設定每幀在x中的位移量位置
inds = (1:len);             % 每幀資料對應1:len
frameout(:) = x(indf(:,ones(1,len))+inds(ones(nf,1),:));   % 對資料分幀
if (nwin > 1)               % 若引數中包括窗函式，把每幀乘以窗函式
    w = win(:)';            % 把win轉成行資料
    frameout = frameout .* w(ones(nf,1),:);  % 乘窗函式
end
 
end

運行結果

三維圖可以旋轉觀察~

俯視圖可以更直觀看出頻率成分隨時間的變化~

轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/252204.html

標籤：其他

上一篇：嵌入式的一點回憶

下一篇：樹莓派CM4和CM4IO上手

數字信號處理（3）- 短時頻域分析

短時頻域分析

短時傅里葉變換

MATLAB程式

運行結果