Python实现指数扫频解卷积提取双耳音频HRIR冲激响应求助

核心疑问解答 #

• FFT相除实现逻辑：音频信号的时域卷积等价于频域的乘积，因此解卷积操作在频域即为两个信号FFT结果的复数除法。Python中直接用/运算符即可完成复数相除，需要加一个极小的偏移量（如1e-16）避免原始扫频FFT出现零值时的除零错误，核心代码为hrir_fft = recorded_fft / (sweep_fft + 1e-16)。
• 丢弃后半段逻辑：我们将信号零填充到2倍长度做FFT是为了避免圆周卷积混叠，逆FFT得到的原始结果中，前半段为因果的有效冲激响应，后半段为谐波失真响应、圆周混叠冗余内容，没有实际意义，因此仅需要截取前N点（N为HRIR的预期长度，48kHz采样下通常取128_{512点，对应2}10ms）即可。
• 加窗逻辑：直接截断冲激响应会引入吉布斯效应，导致频域响应出现不必要的旁瓣，因此需要对截取后的HRIR施加汉宁窗、汉明窗等平滑窗函数，让信号首尾平滑过渡到零，降低截断带来的频域失真。

完整实现代码 #

import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt

# 基础参数设置
fs = 48000  # 采样率
sweep_duration = 1  # 扫频信号时长1s
f_start = 20  # 扫频起始频率
f_end = 20000  # 扫频结束频率
hrir_length = 256  # 预期HRIR长度（48kHz下对应约5.3ms，符合人头HRIR的实际长度范围）

# 1. 生成原始指数扫频信号
t = np.linspace(0, sweep_duration, int(fs*sweep_duration), endpoint=False)
original_sweep = sig.chirp(t, f0=f_start, f1=f_end, t1=sweep_duration, method='logarithmic')

# 2. 模拟录制信号：用模拟HRIR卷积原始扫频，加少量环境噪声
simulated_real_hrir = np.random.randn(hrir_length) * np.exp(-np.linspace(0, 5, hrir_length))  # 模拟指数衰减的真实HRIR
recorded_sweep = sig.convolve(original_sweep, simulated_real_hrir)[:len(original_sweep)] + 0.001 * np.random.randn(len(original_sweep))

# 3. 零填充到2倍原长度，避免圆周卷积混叠
n_fft = 2 * len(original_sweep)
sweep_pad = np.pad(original_sweep, (0, n_fft - len(original_sweep)))
recorded_pad = np.pad(recorded_sweep, (0, n_fft - len(recorded_sweep)))

# 4. 频域解卷积：FFT后相除
sweep_fft = np.fft.fft(sweep_pad)
recorded_fft = np.fft.fft(recorded_pad)
hrir_fft = recorded_fft / (sweep_fft + 1e-16)  # 加极小值避免除零错误

# 5. 逆FFT得到时域冲激响应
hrir_raw = np.fft.ifft(hrir_fft).real  # 取实部，去除计算带来的极小虚数噪声

# 6. 截取有效区域+加窗得到最终HRIR
hrir_causal = hrir_raw[:hrir_length]  # 丢弃后半段冗余内容
window = sig.hann(hrir_length)  # 汉宁窗做平滑截断
hrir_final = hrir_causal * window

结果可视化代码 #

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(10, 12))

# 子图1：原始扫频与录制扫频对比
ax1.plot(t, original_sweep, label='原始指数扫频', alpha=0.7)
ax1.plot(t, recorded_sweep, label='模拟录制扫频', alpha=0.7)
ax1.set_title('扫频信号时域对比')
ax1.set_xlabel('时间(s)')
ax1.legend()

# 子图2：逆FFT原始结果与有效区域标注
ax2.plot(hrir_raw[:1000], label='逆FFT输出原始结果')
ax2.axvspan(hrir_length, 1000, color='gray', alpha=0.3, label='丢弃的后半段冗余区域')
ax2.set_title('逆FFT输出与有效区域截取')
ax2.set_xlabel('采样点')
ax2.legend()

# 子图3：解卷积HRIR与真实HRIR对比
ax3.plot(simulated_real_hrir, label='模拟真实HRIR', alpha=0.7)
ax3.plot(hrir_final, label='解卷积得到的HRIR', alpha=0.7)
ax3.set_title('HRIR结果对比')
ax3.set_xlabel('采样点')
ax3.legend()

plt.tight_layout()
plt.show()

可视化结果说明 #

• 第一幅图展示原始干净扫频和经过HRIR卷积、加噪后的录制信号差异
• 第二幅图可以看到逆FFT输出的前256点是幅值较高的有效冲激响应，后续幅度快速衰减到接近零，灰色区域为需要丢弃的冗余部分
• 第三幅图对比解卷积得到的HRIR和模拟的真实HRIR，二者一致性较高，验证了解卷积流程的正确性

内容的提问来源于stack exchange，提问作者Sanket Jain