%短時傅里葉變換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（2562^nextpow2（nsc））;%N點采樣長度
%也可手動設置
nsc=100;%海明窗的長度即每個窗口的長度
nov=0;%重疊率
nff=max（2562^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（rowncoln）;%初始化最終結果
%對每個片段進行fft變換
index=1;%當前片段第一個信號位置在原始信號中的索引
for?i=1:coln
????%提取當前片段信號值并用海明窗進行加權
????temp_S=S（index:index+nsc-1）.*h‘;
????%進行N點FFT變換
????temp_X=fft（temp_Snff）;
????%取一半
????STFT_X（:i）=temp_X（1:rown）‘;
????%將索引后移
????index=index+（nsc-nov）;
end

%%?matlab自帶函數
spectrogram（Shamming（nsc?‘periodic‘）novnffFs）;
title（‘spectro