STFT短时傅里叶变换代码

### 短时傅里叶变换 (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]。 --- ### 相关概念补充 短时傅里叶变换的核心在于将长时间序列分解为多个较短的时间片段，并分别对其进行傅里叶变换。这种方法适用于非平稳信号的分析，因为它能够在时间和频率两个维度上提供局部化信息。 ---

### 实现短时傅里叶变换 (STFT) 的 MATLAB 代码示例 MATLAB 提供了内置函数 `spectrogram` 来计算并可视化信号的频谱图，这实际上是 STFT 的一种应用形式[^1]。然而，为了更深入理解 STFT 的工作原理，下面展示了一个手动实现 STFT 的 MATLAB 代码： ```matlab function [S, f, t] = stft(x, fs, windowSize, overlap, fftSize) % 输入参数说明： % x : 输入信号向量 % fs : 采样率(Hz) % windowSize: 窗口大小(样本数) % overlap : 每次移动窗口重叠的比例（0到windowSize之间的整数值） % fftSize : FFT长度 win = hamming(windowSize); % 使用汉明窗作为加权窗口 stepSize = windowSize - overlap; % 初始化输出矩阵 S 和时间轴 t numFrames = floor((length(x)-overlap)/(windowSize-overlap)); S = zeros(fftSize,numFrames); t = ((1:numFrames)*stepSize)/fs; for i=1:numFrames frameStart = (i-1)*(windowSize-overlap)+1; frameEnd = min(frameStart+windowSize-1,length(x)); if frameEnd >= length(x), break; end % 对当前帧施加窗口函数，并执行FFT currentFrame = x(frameStart:frameEnd).*win(1:frameEnd-frameStart+1)'; X = fft(currentFrame, fftSize); % 存储幅度谱 S(:,i)=abs(X); end % 计算频率轴 f=(0:(fftSize/2))*fs/fftSize; end ``` 此代码定义了一个名为 `stft` 的函数来计算输入音频片段 `x` 在不同时间段内的频域表示。通过调整 `windowSize`, `overlap` 及其他参数可以改变分辨率及时延特性。