### 实现短时傅里叶变换（STFT）而不使用 Python 或 MATLAB 中的内置 DFT/FFT 函数 实现短时傅里叶变换（STFT）的核心在于对信号进行分帧、加窗，然后对每一帧计算其离散傅里叶变换（DFT）。在不依赖内置函数的情况下，需要手动实现这些步骤。 #### 1. **自定义窗函数** 首先，定义常见的窗函数，例如汉明窗（Hamming）、汉宁窗（Hanning）和余弦斜坡窗（Cosine Ramp）。这些窗函数可以通过数学公式直接实现： - **汉明窗**： $ w(n) = 0.54 - 0.46 \cdot \cos\left(\frac{2\pi n}{N-1}\right), \quad 0 \leq n \leq N-1 $ - **汉宁窗**： $ w(n) = 0.5 - 0.5 \cdot \cos\left(\frac{2\pi n}{N-1}\right), \quad 0 \leq n \leq N-1 $ - **余弦斜坡窗**： $ w(n) = \cos\left(\frac{\pi n}{N-1} - \frac{\pi}{2}\right), \quad 0 \leq n \leq N-1 $ ```python import numpy as np def window_function(window_type, N): n = np.arange(N) if window_type == 'hamming': return 0.54 - 0.46 * np.cos(2 * np.pi * n / (N - 1)) elif window_type == 'hanning': return 0.5 - 0.5 * np.cos(2 * np.pi * n / (N - 1)) elif window_type == 'cosine_ramp': return np.cos(np.pi * n / (N - 1) - np.pi / 2) else: raise ValueError("Unsupported window type") ``` #### 2. **信号分帧** 接下来，将输入信号划分为多个重叠帧。每帧的长度为 `frame_length`，帧之间的跳跃步长为 `hop_size`。如果信号长度不足以填满最后一帧，则可以进行零填充。 ```python def frame_signal(signal, frame_length, hop_size): num_frames = 1 + (len(signal) - frame_length) // hop_size frames = np.zeros((num_frames, frame_length)) for i in range(num_frames): start = i * hop_size end = start + frame_length frames[i] = signal[start:end] return frames ``` #### 3. **自定义 DFT 实现** 为了计算每一帧的频谱，需要实现离散傅里叶变换（DFT）。DFT 的公式如下： $$ X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-j2\pi kn/N}, \quad k = 0, 1, ..., N-1 $$ 其中 $ x[n] $ 是输入信号，$ X[k] $ 是对应的频率成分。 ```python def dft(x): N = len(x) n = np.arange(N) k = n.reshape((N, 1)) M = np.exp(-2j * np.pi * k * n / N) return np.dot(M, x) ``` #### 4. **组合所有步骤以实现 STFT** 最后，将上述功能组合起来，依次完成信号分帧、加窗和 DFT 计算，从而得到 STFT 结果。 ```python def stft(signal, frame_length, hop_size, window_type='hanning'): # Step 1: 分帧 frames = frame_signal(signal, frame_length, hop_size) # Step 2: 加窗 window = window_function(window_type, frame_length) framed_windows = frames * window # Step 3: 对每一帧计算 DFT stft_result = np.array([dft(frame) for frame in framed_windows]) return stft_result ``` #### 5. **示例用法** 以下是一个简单的示例，展示如何使用上述代码实现 STFT： ```python # 示例信号 fs = 1000 # 采样率 t = np.linspace(0, 1, fs, endpoint=False) signal = np.sin(2 * np.pi * 50 * t) + np.sin(2 * np.pi * 120 * t) # 参数设置 frame_length = 256 hop_size = 128 # 计算 STFT stft_result = stft(signal, frame_length, hop_size, window_type='hanning') # 可视化 import matplotlib.pyplot as plt plt.imshow(np.abs(stft_result[:, :frame_length//2+1].T), origin='lower', aspect='auto', cmap='jet') plt.title('STFT Magnitude') plt.xlabel('Frame') plt.ylabel('Frequency Bin') plt.colorbar() plt.show() ``` ---