基于MATLAB的短时傅里叶变换及其逆变换的实现

基于MATLAB的短时傅里叶变换（STFT）及其逆变换（ISTFT）的实现

一、核心代码实现

%% 清空环境
clear; clc; close all;

%% 参数设置
fs = 1000;          % 采样频率
t = 0:1/fs:2;       % 时间向量（2秒）
f1 = 50; f2 = 200;  % 信号频率成分
x = sin(2*pi*f1*t) + 0.5*sin(2*pi*f2*t);  % 原始信号

%% STFT参数
win_len = 256;      % 窗长（影响频率分辨率）
overlap = 128;      % 重叠长度（影响时间分辨率）
nfft = 512;         % FFT点数
window = hamming(win_len);  % 窗函数选择

%% 执行STFT
[S, f, t_stft] = stft_custom(x, fs, window, overlap, nfft);

%% 可视化STFT结果
figure;
imagesc(t_stft, f, 20*log10(abs(S)+eps));
axis xy; shading interp;
xlabel('时间 (s)'); ylabel('频率 (Hz)');
title('STFT时频谱图 (Hamming窗)');
colorbar;

%% 逆变换重建信号
x_recon = istft_custom(S, window, overlap, nfft, fs);

%% 重建效果验证
figure;
subplot(2,1,1);
plot(t, x); hold on;
plot(t, x_recon, '--r');
title('原始信号与重建信号对比');
xlabel('时间 (s)'); ylabel('幅值');
legend('原始信号', '重建信号');

subplot(2,1,2);
plot(t, x - x_recon);
title('重建误差');
xlabel('时间 (s)'); ylabel('误差幅值');

%% 自定义STFT函数
function [S, f, t] = stft_custom(x, fs, window, overlap, nfft)
    % 输入参数：
    % x: 输入信号
    % fs: 采样频率
    % window: 窗函数
    % overlap: 重叠长度
    % nfft: FFT点数

    [win_len, ~] = size(window);
    hop = win_len - overlap;  % 步长

    % 分帧处理
    num_frames = 1 + floor((length(x) - win_len)/hop);
    frames = zeros(num_frames, win_len);

    for i = 1:num_frames
        start_idx = (i-1)*hop + 1;
        frames(i,:) = x(start_idx:start_idx+win_len-1) .* window';
    end

    % 计算STFT
    S = fft(frames, nfft, 2);
    S = S(:,1:win_len/2+1);  % 单边频谱

    % 频率轴
    f = (0:win_len/2) * fs / win_len;
    % 时间轴
    t = (0:num_frames-1) * hop / fs;
end

%% 自定义ISTFT函数（最小二乘优化）
function x_recon = istft_custom(S, window, overlap, nfft, fs)
    % 输入参数：
    % S: STFT频谱
    % window: 窗函数
    % overlap: 重叠长度
    % nfft: FFT点数
    % fs: 采样频率

    [n_freq, n_frames] = size(S);
    win_len = length(window);
    hop = win_len - overlap;

    % 初始化重建信号
    x_recon = zeros((n_frames-1)*hop + win_len, 1);

    % 重叠相加法
    for i = 1:n_frames
        % 逆FFT
        frame = ifft(S(:,i) .* conj(S(:,i)), nfft);
        frame = real(frame(1:win_len));

        % 重叠相加
        start_idx = (i-1)*hop + 1;
        x_recon(start_idx:start_idx+win_len-1) = ...
            x_recon(start_idx:start_idx+win_len-1) + frame;
    end

    % 窗函数能量归一化
    x_recon = x_recon ./ sum(window.^2);
end

二、关键参数说明

三、应用场景示例

1. 语音信号分析

% 加载语音信号
[语音, fs] = audioread('speech.wav');
语音 =语音(:,1); % 单声道

% STFT参数设置
win_len = 512; overlap = 256;

% 计算并绘制语谱图
[S,f,t] = stft_custom(语音, fs, hamming(win_len), overlap, 1024);
imagesc(t, f, 20*log10(abs(S)+eps));
title('语音信号语谱图');

2. 机械故障检测

% 加载振动信号
[vibration, fs] = load_vibration_data();

% STFT分析
[S,f,t] = stft_custom(vibration, fs, hanning(1024), 512, 2048);

% 特征提取（频域能量）
energy = sum(abs(S), 1);
plot(t, energy);
title('振动信号频域能量分布');

四、注意事项

1. 窗函数选择 短时分析需满足常数重叠相加（COLA）条件，推荐使用Hann/Hamming窗 窗长与信号频率成分匹配：高频信号用短窗，低频信号用长窗

2. 频谱泄漏抑制

   • 增加窗长可减少频谱泄漏

   • 使用加窗函数后需进行能量归一化：

     S = S .* sum(window.^2);
     

3. 实时性处理 对实时信号需设置固定帧长和重叠 使用overlap-add方法保证时域连续性

参考代码 短时傅立叶变换及其逆变换源程序 www.youwenfan.com/contentale/78965.html

该方法通过合理选择窗函数和参数设置，实现了高精度的时频分析。