短时傅里叶变换STFT（非使用fft函数）

离散的短时傅里叶变换

Matlab代码

function [STFT, fre,t] = stft(x, win, hop, fs,freqlow,freqhigh,alpha)
% 计算短时傅里叶变换
% Input:
%   x - 一维信号
%   win - 窗函数
%   hop - hop - 移动长度
%   fs - 采样率
%   freqlow - 最小频率
%   freqhigh - 最大频率
%   alpha - 频率分辨率
%
% Output:
%   STFT - STFT-矩阵 [T, F]
%   fre - 频率向量
%   t - 分割后时间向量
% Jinshun Shen,Xidian University, March 2022
%jsshen@stu.xidian.edu.cn
% This program can not be used for commercialization without the authorization of its author. 

% 把 x 变为列向量
x = x(:);
win=win(:);
xlen = length(x);
wlen = length(win);

% 窗口数目 L
L =fix(xlen/hop);
fre=freqlow:alpha:freqhigh;
N=length(fre);
STFT = zeros(N, L);
    
% STFT    
flag1=0;%记录第一阶段加的窗个数
flag2=0;
for i = 0:L-1
    % 加窗
    if i*hop <= 0.5*wlen%刚开始，窗必须从最高的地方（窗的中心）开始，
        %否则0.5*wlen长度的信号的频率会消失，因为gauss窗的前0.5*wlen权重小
        K=round(0.5*wlen+i*hop);%窗的宽度,最大为wlen
        
        xw = x(1:K).*win(wlen-K+1:end);
        t=(1: K)/fs;
        t=t(:);
        for k=1:N
            f=fre(k);
             STFT(k,i+1)=sum(xw.*exp(-1i*2*pi*f*t));
        end
        flag1=flag1+1;
    elseif i*hop > 0.5*wlen && i*hop<= xlen-0.5*wlen%窗函数不需要截断
        xw = x(1+(i-flag1)*hop : wlen+(i-flag1)*hop).*win;
        t=(1+(i-flag1)*hop : wlen+(i-flag1)*hop)/fs;
        t=t(:);
        for k=1:N
            f=fre(k);
             STFT(k,i+1)=sum(xw.*exp(-1i*2*pi*f*t));
        end
        flag2=flag2+1;
    elseif i*hop > xlen-0.5*wlen
        K=wlen-(i-flag1-flag2)*hop;%窗的宽度,最大为wlen
        xw = x(xlen-K+1:end).*win(1:K);
        t=(xlen-K+1:xlen)/fs;
        t=t(:);
        for k=1:N
            f=fre(k);
             STFT(k,i+1)=sum(xw.*exp(-1i*2*pi*f*t));
        end
    end
end
% 取每个窗口中点的时间点
    t = (0:hop:xlen-1)/fs;  
end

实验结果