基于傅里叶变换求解系统脉冲响应的异常问题及修正需求

问题背景 #

我需要通过傅里叶变换计算系统的脉冲响应，已有存储在.npy文件中的输入、输出矩阵。傅里叶变换实现代码得到了错误的脉冲响应结果，但互相关方法的代码得到了预期结果，现需修正傅里叶变换方法的问题，确保得到正确结果。

错误的傅里叶变换实现代码 #

import numpy as np
import matplotlib.pyplot as plt
import os
import scipy.signal as sg


working_directory = os.path.dirname(os.path.abspath(__file__))
t_path = os.path.join(working_directory, 't.npy')
x_t_path = os.path.join(working_directory, 'x_t.npy')
y_t_path = os.path.join(working_directory, 'y_t.npy')

t = np.load(t_path)
x_t = np.load(x_t_path)
y_t = np.load(y_t_path)
x_t[np.abs(x_t) < 1e-10] = 1e-10


# Compute the Fourier transforms of the input and output
X_f = np.fft.fft(x_t)
Y_f = np.fft.fft(y_t)

# Compute the frequency response
H_f = np.divide(Y_f, X_f, out=np.zeros_like(Y_f), where=X_f!=0)

# Regularize the frequency response to avoid division by zero
H_f[np.abs(H_f) < 1e-10] = 1e-10

# Compute the impulse response
h_t = np.fft.ifft(H_f)

正确的互相关实现代码 #

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


working_directory = os.path.dirname(os.path.abspath(__file__))
t_path = os.path.join(working_directory, 't.npy')
x_t_path = os.path.join(working_directory, 'x_t.npy')
y_t_path = os.path.join(working_directory, 'y_t.npy')

t = np.load(t_path)
x_t = np.load(x_t_path)
y_t = np.load(y_t_path)

# Ensure the input and output signals are aligned and have the same size
assert len(x_t) == len(y_t), "Input and output signals must have the same size"

# Calculate the cross-correlation between the input and output signals
cross_corr = np.correlate(y_t, x_t, mode='full')

# Normalize the cross-correlation to obtain the impulse response
impulse_response = cross_corr / np.sum(x_t**2)

# Calculate the time vector for the impulse response
dt = t[1] - t[0]
t_impulse = np.arange(-len(x_t) + 1, len(x_t)) * dt

相关图表说明 #

• 时域输入输出图：输入输出信号的时域波形
• 频域输入输出图：输入输出信号的频域频谱
• 错误的脉冲响应时域图：傅里叶变换法得到的错误脉冲响应波形
• 预期的脉冲响应图：互相关法得到的正确脉冲响应波形

问题分析与修正方案 #

傅里叶变换法出错的核心原因是混淆了循环卷积与线性卷积的差异，同时缺少归一化和正确的后处理步骤：

1. 循环卷积混叠：np.fft.fft默认计算的是循环卷积，而实际系统的输入输出关系是线性卷积，必须将信号补零到至少len(x)+len(y)-1的长度，避免时域混叠。
2. 未处理实信号特性：脉冲响应是实信号，逆傅里叶变换的虚部为数值计算误差，需提取实部。
3. 缺少归一化：互相关法中做了除以np.sum(x_t**2)的归一化，傅里叶变换法需同步该步骤保证幅值匹配。
4. 时间轴未对齐：未生成与互相关法一致的时间轴，导致结果无法直接对比。

修正后的傅里叶变换代码 #

import numpy as np
import os

working_directory = os.path.dirname(os.path.abspath(__file__))
t_path = os.path.join(working_directory, 't.npy')
x_t_path = os.path.join(working_directory, 'x_t.npy')
y_t_path = os.path.join(working_directory, 'y_t.npy')

t = np.load(t_path)
x_t = np.load(x_t_path)
y_t = np.load(y_t_path)
dt = t[1] - t[0]
N = len(x_t)

# 补零到线性卷积所需长度（2N-1），避免循环卷积混叠
pad_length = 2 * N - 1
x_padded = np.zeros(pad_length)
x_padded[:N] = x_t
y_padded = np.zeros(pad_length)
y_padded[:N] = y_t

# 计算补零后的傅里叶变换
X_f = np.fft.fft(x_padded)
Y_f = np.fft.fft(y_padded)

# 计算频率响应，严格过滤极小值避免除零
H_f = np.divide(Y_f, X_f, out=np.zeros_like(Y_f), where=np.abs(X_f) > 1e-10)

# 逆傅里叶变换后提取实部（实信号的虚部为数值误差）
h_t = np.fft.ifft(H_f).real

# 与互相关法保持一致的归一化步骤
h_t = h_t / np.sum(x_t**2)

# 生成对齐的时间轴
t_impulse = np.arange(-N + 1, N) * dt

验证说明 #

修正后的代码会输出与互相关法一致的脉冲响应：

• 补零操作确保了线性卷积的计算，消除了循环混叠导致的波形畸变；
• 实部提取和归一化保证了结果的幅值和物理意义正确；
• 对齐的时间轴可直接与互相关法的结果对比验证。

内容的提问来源于stack exchange，提问作者anas-abseh