手动操作实现：信号的分帧、加窗、傅里叶变换和逆变换，重建音频信号

niuguangshuo

已于 2024-11-18 17:47:57 修改

阅读量2.3k

点赞数 23

分类专栏：音频算法文章标签：音视频 android

于 2024-11-12 20:05:01 首次发布

本文链接：https://blog.youkuaiyun.com/qq_34941290/article/details/143722579

版权

音频算法专栏收录该内容

13 篇文章

订阅专栏

引言

在音频信号处理中，信号的分析和重建是核心任务之一。通过将信号分帧、加窗、进行傅里叶变换（STFT）以及逆变换，我们能够提取信号的频率特征并重建出原始信号。本文将逐步介绍这一过程，并展示每个步骤的可视化结果。

1. 原始信号与分帧

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

为什么要分帧？

分帧的主要原因是音频信号通常是非平稳的，即信号的特性会随着时间变化。通过将信号切分为小的时间段（帧），我们能够在每一帧内近似地认为信号是静态的，从而使得频率分析更为准确。这样，我们可以捕捉到信号在不同时间段的频率特征，便于后续的处理和分析。

帧之间的重叠

帧之间的重叠是为了减少在分帧处理时可能丢失的边界信息。当信号被切分为帧时，帧的边缘可能会导致信息的丢失，特别是在信号变化较快的情况下。通过让帧之间有一定的重叠，我们可以确保在分析过程中更好地保留信号的连续性和完整性，从而提高频率分析的精度。

在这里插入图片描述

2. 分帧信号加窗

在音频信号处理中，分帧是将连续的信号切分为若干个小的时间段（帧），以便于分析。然而，直接对每一帧进行傅里叶变换可能会导致频谱泄漏现象。这是因为在信号的每一帧的边缘，信号值可能会发生突变，导致频谱中出现不必要的伪影。

加窗的作用

加窗的主要作用是平滑每一帧的信号，减少在傅里叶变换过程中出现的频谱泄漏。具体来说，加窗有以下几个重要作用：

减少频谱泄漏：通过对每一帧施加窗函数，可以降低帧边缘的突变，使得信号在帧的开始和结束处的值更接近于零。这种平滑处理有助于减少在频域中出现的伪影，从而提升频谱的准确性。
提高频率分辨率：窗函数能够更好地保留信号的频率成分，从而提高频率分析的精度。经过加窗的信号在频域中的表现会更加清晰，能够有效区分不同频率成分。
控制能量分布：不同的窗函数（如汉明窗、汉宁窗、矩形窗等）具有不同的特性，它们可以影响信号在频域中的能量分布。选择合适的窗函数可以优化信号处理的效果，确保重要频率成分的能量得到适当的保留。
改善重建效果：在进行逆变换时，加窗有助于更好地恢复信号的幅度，确保重建信号尽可能接近于原始信号。通过加窗，我们可以在重建过程中有效地修正信号的幅度，从而提高重建的准确性。

重叠与加窗的联系

重叠和加窗是音频信号处理中的两个重要概念，它们相辅相成。在进行分帧时，重叠可以确保在帧之间保留更多的信号信息，而加窗则可以平滑每一帧的边缘，减少频谱泄漏。通过重叠和加窗的结合，我们能够更有效地分析信号的频率特征，并在后续的重建过程中获得更高的准确性。

常见的窗函数

在这里插入图片描述

1. 汉明窗（Hamming Window）

$\cos\left(\frac{2\pi n}{N-1}\right), \quad n = 0, 1, \ldots, N-1$

其中：

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

2. 汉宁窗（Hanning Window）

$\frac{1}{2} \left(1 - \cos\left(\frac{2\pi n}{N-1}\right)\right), \quad n = 0, 1, \ldots, N-1$

其中：

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

3. 矩形窗（Rectangular Window）

公式：
矩形窗的公式为：

$\quad n = 0, 1, \ldots, N-1$

其中：

( N ) 是窗的长度。
( n ) 是当前样本的索引。
汉明窗 和 汉宁窗 都是加窗处理的常用窗函数，能够有效减少频谱泄漏，且在频域分析中表现良好。
矩形窗 不进行加窗处理，直接使用帧信号，但容易导致频谱泄漏，因此在实际应用中应谨慎使用。

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

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

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

$\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）是将频域信息转换回时域的过程。通过对每一帧的频谱进行逆变换，我们可以重建信号。在此过程中，需要注意恢复信号的幅度。
在这里插入图片描述

在重建信号时，我们对频域信号进行逆变换，得到时域信号。具体公式为：

$\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. 恢复原始信号幅度

在进行逆变换后，重建的信号幅度需要通过窗函数进行恢复。由于窗函数在信号处理过程中影响了信号的幅度，因此在重建信号时，必须对恢复的幅度进行适当的归一化处理，以确保信号的真实幅度。
在这里插入图片描述

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

在重建完整信号时，需要考虑帧之间的重叠。对于重叠的样本，必须根据重叠次数对幅度进行归一化，以确保每个样本的幅度恢复正确。通过这种方式，我们可以得到更接近于原始信号的重建结果。
在这里插入图片描述

7. 重建信号和原始信号差异对比

最后，我们可以通过对比重建信号与原始信号之间的差异来评估重建效果。通过绘制两者之间的差异，我们能够直观地观察到重建过程的准确性和有效性。
在这里插入图片描述

完整代码

import numpy as np
import matplotlib.pyplot as plt
import librosa

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

# 使用 librosa.util.frame 切分信号
frame_length = 320  # 帧长度
hop_length = 160  # 帧之间的跳跃长度
frames = librosa.util.frame(y, frame_length=frame_length, hop_length=hop_length)

# 使用窗函数（汉明窗）
window = np.hanning(frame_length)

# 可视化
num_frames = frames.shape[1]  # 帧的数量

# 原始信号
plt.figure(figsize=(12, 2 * (num_frames + 1)))
plt.subplot(num_frames + 1, 1, 1)
t = np.linspace(0, duration, num_samples)
plt.plot(t, y, label='Original Signal', color='blue')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.xlim(0, 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.xlim(0, 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, 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'Windowed Frame {i + 1}', linestyle='--')
    plt.ylabel('Amplitude')
    plt.xlim(0, duration)
    plt.ylim(-3, 3)  # 设置 y 轴范围

# 绘制加窗后每个帧的频谱，重建信号
plt.figure(figsize=(12, 2 * (num_frames + 1)))
reconstructed_frames = []
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)

    # 重建信号
    reconstructed_frame = np.fft.ifft(spectrum).real
    reconstructed_frames.append(reconstructed_frame)

    # 绘制频谱
    plt.plot(freqs[:len(freqs) // 2], np.abs(spectrum)[:len(spectrum) // 2], alpha=0.7,
             label=f'Spectrum of Frame {i + 1}')
    plt.ylabel('Magnitude')
    plt.xlim(0, sr / 2)  # 只显示到 Nyquist 频率
plt.xlabel('Frequency (Hz)')

# 绘制每一帧加窗的重建信号
plt.figure(figsize=(12, 2 * (num_frames + 1)))  # 增加子图的高度
for i in range(num_frames):
    plt.subplot(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)

    # 绘制重建后的信号
    plt.plot(t_frame, reconstructed_frames[i], alpha=0.7, label=f'Reconstructed Frame {i + 1}', linestyle='--')
    plt.ylabel('Amplitude')
    plt.xlim(0, duration)
    plt.ylim(-3, 3)  # 设置 y 轴范围
plt.xlabel('Time (s)')

# 利用窗函数恢复幅度
plt.figure(figsize=(12, 2 * (num_frames + 1)))  # 增加子图的高度
for i in range(num_frames):
    plt.subplot(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)

    # 恢复幅度
    recovered_frame = reconstructed_frames[i] / (window + 1e-10)  # 避免除零错误

    # 绘制恢复后的信号
    plt.plot(t_frame, recovered_frame, alpha=0.7, label=f'Recovered Frame {i + 1}', linestyle='--')
    plt.ylabel('Amplitude')
    plt.xlim(0, duration)
    plt.ylim(-3, 3)  # 设置 y 轴范围
plt.xlabel('Time (s)')

# 重建完整信号
reconstructed_signal = np.zeros(num_samples)
# 计算重叠次数
overlap_count = (frame_length - hop_length) // hop_length + 1
for i in range(num_frames):
    start = i * hop_length
    recovered_frame = reconstructed_frames[i] / (window + 1e-12)
    if i == 0:
        recovered_frame[hop_length:] /= overlap_count
    elif i == num_frames - 1:
        recovered_frame[:frame_length - hop_length] /= overlap_count
        pass
    else:
        recovered_frame = recovered_frame / overlap_count
    # 处理重叠部分
    reconstructed_signal[start:start + frame_length] += recovered_frame

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

plt.subplot(2, 1, 2)
plt.plot(np.linspace(0, duration, num_samples), reconstructed_signal, label='Reconstructed Signal', color='orange')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.xlim(0, duration)

plt.tight_layout()

# 对比重建信号与原始信号的差异
plt.figure(figsize=(12, 6))
plt.plot(np.linspace(0, duration, num_samples), y - reconstructed_signal, 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, duration)
plt.ylim(-3, 3)  # 设置 y 轴范围
plt.legend()
plt.tight_layout()
plt.show()