短时傅里叶变换代码

### 短时傅里叶变换 (STFT) 的实现代码 以下是基于 MATLAB 和 Python 的两种短时傅里叶变换 (STFT) 实现方法： #### 基于 MATLAB 的 STFT 实现 MATLAB 提供了内置函数 `spectrogram` 来计算信号的短时傅里叶变换。如果需要手动实现，可以按照以下方式编写代码。 ```matlab function [S, f, t] = stft(x, fs, windowSize, overlap, fftSize) % 输入参数说明： % x: 输入信号 % fs: 采样率 % windowSize: 时间窗口大小 % overlap: 窗口重叠量 % fftSize: FFT 大小 win = hamming(windowSize); % 使用汉明窗 stepSize = windowSize - overlap; % 计算步长 numFrames = floor((length(x) - overlap) / stepSize); S = []; % 初始化频谱矩阵 for i = 1:numFrames startIdx = (i-1)*stepSize + 1; endIdx = startIdx + windowSize - 1; frame = x(startIdx:endIdx).*win'; % 应用窗函数 specFrame = fft(frame, fftSize); % 对帧进行FFT S = [S, specFrame]; % 将当前帧的频谱加入总频谱矩阵 end f = (0:fftSize/2)*(fs/fftSize); % 频率轴 t = ((windowSize/2):stepSize:(length(x)-windowSize/2))/fs; % 时间轴 ``` 此代码实现了自定义的短时傅里叶变换功能[^3]。 --- #### 基于 Python 的 STFT 实现 Python 中可以通过 NumPy 和 SciPy 库来实现短时傅里叶变换。下面是一个简单的示例代码： ```python import numpy as np from scipy.signal import get_window import matplotlib.pyplot as plt def stft(signal, sample_rate, window_size=1024, hop_length=512, fft_size=1024): """ 参数说明： signal: 输入音频信号 sample_rate: 采样率 window_size: 窗口大小 hop_length: 步幅长度 fft_size: FFT 变换大小 返回值： spectrogram: 谱图 freqs: 频率向量 times: 时间向量 """ win = get_window('hamming', window_size) # 获取汉明窗 num_frames = int(np.ceil(len(signal)/hop_length)) # 总帧数 pad_signal_len = num_frames * hop_length + window_size // 2 padded_signal = np.pad(signal, (int(window_size//2), int(pad_signal_len - len(signal))), mode='constant') indices = np.arange(0, window_size).reshape(-1, 1) + np.arange(0, num_frames*hop_length, hop_length).reshape(1,-1) frames = padded_signal[indices.astype(int)] # 切分帧 framed_signal = frames.T * win # 加窗 spectrum = np.fft.rfft(framed_signal, n=fft_size, axis=-1) # 进行FFT magnitude_spectrum = np.abs(spectrum) # 幅度谱 time_axis = np.linspace(0, len(signal)/sample_rate, num_frames) # 时间轴 frequency_axis = np.linspace(0, sample_rate/2, fft_size//2+1) # 频率轴 return magnitude_spectrum, frequency_axis, time_axis # 测试代码 if __name__ == "__main__": sr = 44100 # 采样率 duration = 1.0 # 持续时间 t = np.linspace(0., duration, int(sr * duration)) audio_signal = 0.5*np.sin(2 * np.pi * 440 * t) # 合成正弦波信号 mag_spec, freqs, times = stft(audio_signal, sr) plt.figure(figsize=(8, 6)) plt.pcolormesh(times, freqs, 20 * np.log10(mag_spec + 1e-7), shading='gouraud') # 绘制对数幅度谱 plt.ylabel('Frequency [Hz]') plt.xlabel('Time [sec]') plt.title('Short-Time Fourier Transform Spectrogram') plt.colorbar(label='Magnitude [dB]') plt.show() ``` 这段代码展示了如何通过 Python 手动实现短时傅里叶变换并绘制对应的频谱图[^1]。 --- ### 相关概念补充 短时傅里叶变换的核心在于将长时间序列分解为多个较短的时间片段，并分别对其进行傅里叶变换。这种方法适用于非平稳信号的分析，因为它能够在时间和频率两个维度上提供局部化信息。 ---