基于NumPy手动实现：信号的分帧、加窗、傅里叶变换 / 逆变换、重建

1. 原始信号与分帧

在音频信号处理之前，我们生成一个随机信号作为示例。原始信号在时域上表现为波动。为了进行频率分析，我们将信号切分为多个小的时间段，称为帧。每一帧将被单独分析，以便在每个小段时间内假设信号是静态的。

为什么要分帧？

音频信号通常是非平稳的，即信号特性随时间变化。通过将信号切分为小的时间段（帧），我们可以在每一帧内近似认为信号是静态的，从而使频率分析更为准确。这种方法帮助我们捕捉信号在不同时间段的频率特征，便于后续处理和分析。

帧之间的重叠

帧之间的重叠可以减少在分帧处理时可能丢失的边界信息。信号的快速变化可能导致信息的丢失，而通过让帧之间有一定的重叠，我们可以更好地保留信号的连续性和完整性，从而提高频率分析的精度。

2. 分帧信号加窗

在音频信号处理中，分帧是将连续信号切分为若干小时间段（帧），以便于分析。然而，直接对每一帧进行傅里叶变换可能导致频谱泄漏，尤其在帧的边缘，信号值的突变会引入伪影。加窗的主要目的是平滑每一帧的信号，从而减少频谱泄漏。

1. 减少频谱泄漏：
   在进行傅里叶变换时，假设信号在分析窗口内是周期性的。如果信号在窗口的边缘处发生突变（例如，信号值突然从一个值跳到另一个值），则傅里叶变换会将这种突变视为频率成分，导致频谱中出现伪影。通过施加窗函数，信号在帧的开始和结束处的值逐渐接近于零，从而减少这种突变，降低频谱泄漏的风险。

2. 提高频率分辨率：
   窗函数的形状直接影响频谱的表现。通过加窗，信号的频率成分能够更好地保留，经过加窗的信号在频域中的表现更加清晰。这使得不同频率成分之间的区分更加有效，尤其是在频率接近的情况下。

3. 控制能量分布：
   不同的窗函数（如汉明窗、汉宁窗、矩形窗等）具有不同的特性，影响信号在频域中的能量分布。例如，汉明窗和汉宁窗能够有效地减少边缘效应，而矩形窗则可能导致较大的频谱泄漏。选择合适的窗函数可以优化信号处理效果，确保重要频率成分的能量得到适当保留。

4. 改善重建效果：
   在进行逆变换时，加窗有助于更好地恢复信号的幅度。通过在每一帧上应用窗函数，重建信号时可以有效地修正信号的幅度，从而提高重建的准确性。这对于确保重建信号尽可能接近原始信号至关重要。

重叠确保在帧之间保留更多信号信息，而加窗则平滑每一帧的边缘，减少频谱泄漏。通过结合重叠和加窗，我们能够更有效地分析信号的频率特征，并在后续重建过程中获得更高的准确性。

常见的窗函数

1. 汉明窗（Hamming Window）

$$w[n]=0.54-0.46\cos \left(\frac{2\pi n}{N-1}\right),n=0,1,\dots ,N-1$$

其中：

• ( N ) 是窗的长度。
• ( n ) 是当前样本的索引。

2. 汉宁窗（Hanning Window）

$$w[n]=\frac{1}{2}\left(1-\cos \left(\frac{2\pi n}{N-1}\right)\right),n=0,1,\dots ,N-1$$

其中：

• ( N ) 是窗的长度。
• ( n ) 是当前样本的索引。

3. 矩形窗（Rectangular Window）

公式：
矩形窗的公式为：

$$w[n]=1,n=0,1,\dots ,N-1$$

其中：

• ( N ) 是窗的长度。

• ( n ) 是当前样本的索引。

• 汉明窗 和 汉宁窗 都是加窗处理的常用窗函数，能够有效减少频谱泄漏，且在频域分析中表现良好。

• 矩形窗 不进行加窗处理，直接使用帧信号，但容易导致频谱泄漏，因此在实际应用中应谨慎使用。

下面是对分帧信号，应用汉明窗的结果。

3. 短时傅里叶变换（STFT）

短时傅里叶变换（STFT）将每一帧信号转换到频域，从而提取频率特征。通过对加窗后的帧进行傅里叶变换，我们能够观察信号在时间和频率上的变化，获得频谱信息。

对于每一帧的信号 ( x[n] )，我们对其应用窗函数 ( w[n] ) 后进行傅里叶变换，得到频域表示。具体公式为：

$$X[k]=\sum _{n=0}^{N-1}x[n]w[n]e^{-j\frac{2\pi }{N}kn}$$

其中：

• ( X[k] ) 是第 ( k ) 个频率分量的复数值。
• ( x[n] ) 是当前帧的时域信号。
• ( w[n] ) 是窗函数（如汉明窗或汉宁窗）。
• ( N ) 是窗的长度。
• ( j ) 是虚数单位。

在代码中，这一过程对应于以下行：

spectrum = np.fft.fft(windowed_frame)

4. 逆短时傅里叶变换（ISTFT）

逆短时傅里叶变换（ISTFT）是将频域信息转换回时域的过程。通过对每一帧的频谱进行逆变换，我们可以重建信号。为了保证加窗后，首尾帧信号的完整性，需要在首、尾各自填充 1 帧静音信号。

在重建信号时，我们对频域信号进行逆变换，得到时域信号。具体公式为：
$$x[n]=\frac{1}{N}\sum _{k=0}^{N-1}X[k]e^{j\frac{2\pi }{N}kn}$$

其中：

• ( x[n] ) 是重建的时域信号。
• ( X[k] ) 是频域信号的复数值。
• ( N ) 是窗的长度。

在代码中，这一过程对应于以下行：

reconstructed_frame = np.fft.ifft(spectrum).real

5. 重建原始信号并正确处理重叠部分

OLA（Overlap-Add，重叠加法）是一种常用的信号处理技术，主要用于将分帧处理的信号重建为完整信号。对于每一帧，将其加到重建信号的相应位置，并考虑重叠部分的处理。

reconstructed_signal = np.zeros(all_samples)
for i in range(all_num_frames):
    start = i * hop_length
    reconstructed_signal[start:start + frame_length] += reconstructed_frames[i] /n

对于汉宁窗（Hanning window），如果帧移（hop length）等于帧长的一半，那么在重建信号时，可以直接叠加每一帧的信号，也就是 n=1。
对于汉宁窗（Hanning window），如果帧移（hop length）等于帧长的1/4，那么在重建信号时，叠加后的信号需要乘以 1/2，，也就是 n=2。

完整代码

import matplotlib.pyplot as plt
import numpy as np

# 生成信号
sr = 16000  # 采样率
signal_duration = 0.04  # 持续时间（0.04秒）
signal_samples = int(sr * signal_duration)  # 样本数
y = np.random.randn(signal_samples)  # 生成均值为0，标准差为1的随机信号

frame_length = 320  # 帧长度
hop_length = 160  # 帧之间的跳跃长度
num_frames = (len(y) - frame_length) // hop_length + 1
frames = np.array([y[i * hop_length:i * hop_length + frame_length] for i in range(num_frames)]).T
# 使用窗函数（汉明窗）
window = np.hanning(frame_length)

# 原始信号
plt.figure(figsize=(12, 2 * (num_frames + 1)))
plt.subplot(num_frames + 1, 1, 1)
t = np.linspace(0, signal_duration, signal_samples)
plt.plot(t, y, label='Original Signal', color='blue')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.xlim(0, signal_duration)
plt.ylim(-3, 3)  # 设置 y 轴范围
# 绘制切分的帧
for i in range(num_frames):
    plt.subplot(num_frames + 1, 1, i + 2)  # 从第二个子图开始
    t_start = i * hop_length / sr
    t_end = t_start + frame_length / sr
    t_frame = np.linspace(t_start, t_end, frame_length)

    # 绘制帧
    plt.plot(t_frame, frames[:, i], alpha=0.7, label=f'Frame {i + 1}', linestyle='--')
    plt.ylabel('Amplitude')
    plt.legend()
    plt.xlim(0, signal_duration)
    plt.ylim(-3, 3)  # 设置 y 轴范围

# 绘制切分的帧加窗后的信号
plt.figure(figsize=(12, 2 * (num_frames + 1)))
plt.subplot(num_frames + 1, 1, 1)
plt.plot(t, y, label='Original Signal', color='blue')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.xlim(0, signal_duration)
plt.ylim(-3, 3)  # 设置 y 轴范围

for i in range(num_frames):
    plt.subplot(num_frames + 1, 1, i + 2)  # 从第二个子图开始
    # 应用窗函数
    windowed_frame = frames[:, i] * window
    t_start = i * hop_length / sr
    t_end = t_start + frame_length / sr
    t_frame = np.linspace(t_start, t_end, frame_length)

    # 绘制加窗后的信号
    plt.plot(t_frame, windowed_frame, alpha=0.7, label=f'Hanning Windowed Frame {i + 1}', linestyle='--')
    plt.ylabel('Amplitude')
    plt.legend()
    plt.xlim(0, signal_duration)
    plt.ylim(-3, 3)  # 设置 y 轴范围

# 绘制加窗后每个帧的频谱
plt.figure(figsize=(12, 2 * (num_frames + 1)))
for i in range(num_frames):
    plt.subplot(num_frames, 1, i + 1)  # 每个帧的频谱
    # 应用窗函数
    windowed_frame = frames[:, i] * window
    # 计算 FFT
    spectrum = np.fft.fft(windowed_frame)
    freqs = np.fft.fftfreq(len(spectrum), 1 / sr)
    # 绘制频谱
    plt.plot(freqs[:len(freqs) // 2], np.abs(spectrum)[:len(spectrum) // 2], alpha=0.7,
             label=f'Spectrum of Hanning Windowed Frame {i + 1}')
    plt.ylabel('Magnitude')
    plt.xlim(0, sr / 2)  # 只显示到 Nyquist 频率
    plt.legend()
plt.xlabel('Frequency (Hz)')

padding_duration = 0.01
padding_samples = int(sr * padding_duration)
y_padding = np.pad(y, (padding_samples, padding_samples), mode='constant', constant_values=0)
all_duration = signal_duration + padding_duration * 2
all_samples = signal_samples + padding_samples * 2
all_num_frames = (len(y_padding) - frame_length) // hop_length + 1
all_frames = np.array([y_padding[i * hop_length:i * hop_length + frame_length] for i in range(all_num_frames)]).T
# 绘制每一帧加窗的重建信号
plt.figure(figsize=(12, 1.5 * (all_num_frames + 1)))  # 增加子图的高度
reconstructed_frames = []
for i in range(all_num_frames):
    plt.subplot(all_num_frames, 1, i + 1)  # 每个帧的重建信号
    t_start = i * hop_length / sr
    t_end = t_start + frame_length / sr
    t_frame = np.linspace(t_start, t_end, frame_length)
    windowed_frame = all_frames[:, i] * window
    spectrum = np.fft.fft(windowed_frame)
    freqs = np.fft.fftfreq(len(spectrum), 1 / sr)
    reconstructed_frame = np.fft.ifft(spectrum).real
    reconstructed_frames.append(reconstructed_frame)
    # 绘制重建后的信号
    plt.plot(t_frame, reconstructed_frame, alpha=0.7,
             label=f'Reconstructed Hanning Windowed Frame {i + 1} After Padding',
             linestyle='--')
    plt.ylabel('Amplitude')
    plt.legend()
    plt.xlim(0, all_duration)
    plt.ylim(-3, 3)  # 设置 y 轴范围
plt.xlabel('Time (s)')
#
# 利用窗函数重建完整信号
reconstructed_signal = np.zeros(all_samples)
for i in range(all_num_frames):
    start = i * hop_length
    reconstructed_signal[start:start + frame_length] += reconstructed_frames[i]

# 10. 可视化重建的完整信号
plt.figure(figsize=(12, 6))
plt.subplot(2, 1, 1)
plt.plot(np.linspace(0, signal_duration, signal_samples), y_padding[padding_samples:-padding_samples],
         label='Original Signal', color='blue', alpha=0.5)
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.xlim(0, signal_duration)

plt.subplot(2, 1, 2)
plt.plot(np.linspace(0, signal_duration, signal_samples), reconstructed_signal[padding_samples:-padding_samples],
         label='Reconstructed Signal', color='orange')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.xlim(0, signal_duration)

plt.tight_layout()

# 对比重建信号与原始信号的差异
plt.figure(figsize=(12, 6))
plt.plot(np.linspace(0, signal_duration, signal_samples),
         y_padding[padding_samples:-padding_samples] - reconstructed_signal[padding_samples:-padding_samples],
         label='Difference', color='red', alpha=0.7)
plt.xlabel('Time (s)')
plt.ylabel('Amplitude Difference')
plt.title('Difference between Original and Reconstructed Signal')
plt.xlim(0, signal_duration)
plt.ylim(-3, 3)  # 设置 y 轴范围
plt.legend()
plt.tight_layout()
plt.show()