这个问题很典型，我来拆解下理论可行性、实操难点，再给你解决代码里除法的困惑：

理论层面：有前提的可行性 #

这个方法**仅在系统是线性时不变（LTI）**的前提下成立：

• LTI系统的核心特性就是输出频谱等于输入频谱乘以系统频率响应H(w)，即 (A'(w) = H(w) \cdot A(w))，因此理论上可以通过 (H(w) = A'(w)/A(w)) 反推系统特性，再用 (B(w) = B'(w)/H(w)) 还原输入B，最后逆傅里叶变换得到时域信号。
• 如果系统不是LTI（比如存在非线性、时变特性，或者有记忆性的非线性变换），那么H(w)不是固定的频率响应，这个方法完全失效。

实操层面：存在诸多限制 #

理论可行不代表实操顺畅，实际落地会遇到几个关键问题：

• 频谱零点/极小值问题：输入信号A的频谱A(w)可能在某些频率点为0（或因噪声接近0），直接做除法会导致数值爆炸，还原信号完全失真。
• 噪声放大：实际采集的信号A、A'、B'都会带噪声，当A(w)幅值较小时，噪声会被除法操作放大，严重干扰B的还原精度。
• 频谱泄漏误差：用DFT/FFT计算频谱时，窗函数会导致频谱泄漏，使得估计出的H(w)和真实系统响应有偏差，误差积累后还原的B会偏离真实值。
• 非理想LTI系统：现实中几乎没有严格的LTI系统，轻微的时变、非线性都会让H(w)的估计出现偏差，进一步影响还原效果。

代码中除法操作的解决思路 #

不管用Python还是Matlab，直接逐点做频谱除法都会踩除以0的坑，这里给你实用的处理方案：

Python 实现示例（基于numpy/scipy） #

import numpy as np
from scipy.signal import butter, lfilter
import matplotlib.pyplot as plt

# 1. 生成模拟信号与系统输出
fs = 1000  # 采样率
t = np.linspace(0, 1, fs, endpoint=False)
A = np.sin(2*np.pi*50*t) + np.sin(2*np.pi*120*t)  # 已知输入A

# 模拟未知LTI系统（4阶低通滤波）
b, a = butter(4, 150, 'low', fs=fs)
A_prime = lfilter(b, a, A)  # A的输出A'
B_true = np.sin(2*np.pi*80*t) + np.sin(2*np.pi*180*t)  # 真实输入B
B_prime = lfilter(b, a, B_true)  # B的输出B'

# 2. 计算频谱
N = len(A)
A_fft = np.fft.fft(A)
A_prime_fft = np.fft.fft(A_prime)
B_prime_fft = np.fft.fft(B_prime)

# 3. 安全处理除法：添加正则项避免除以0
epsilon = 1e-8  # 极小的正数，根据实际噪声调整
H_fft = A_prime_fft / (A_fft + epsilon)  # 估计系统频率响应
B_fft = B_prime_fft / (H_fft + epsilon)  # 还原B的频谱

# 4. 逆傅里叶变换得到时域信号
B_recon = np.fft.ifft(B_fft).real  # 取实部，消除数值误差带来的虚部

# 5. 可视化对比
plt.figure(figsize=(12, 6))
plt.subplot(211)
plt.plot(t, B_true, label='真实输入B', linewidth=2)
plt.plot(t, B_recon, label='还原的B', alpha=0.7, linestyle='--')
plt.legend()
plt.title('时域信号对比')
plt.subplot(212)
freq = np.fft.fftfreq(N, 1/fs)
plt.plot(freq[:N//2], np.abs(B_fft)[:N//2], label='还原的B频谱')
plt.xlim(0, 200)
plt.legend()
plt.title('频谱对比')
plt.tight_layout()
plt.show()

Matlab 实现示例 #

fs = 1000;
t = linspace(0, 1, fs, 'Endpoint', false);
A = sin(2*pi*50*t) + sin(2*pi*120*t); % 已知输入A

% 模拟未知LTI系统（零相位低通滤波，避免相位失真）
[b, a] = butter(4, 150/(fs/2), 'low');
A_prime = filtfilt(b, a, A); % A的输出A'
B_true = sin(2*pi*80*t) + sin(2*pi*180*t); % 真实输入B
B_prime = filtfilt(b, a, B_true); % B的输出B'

% 计算频谱
N = length(A);
A_fft = fft(A);
A_prime_fft = fft(A_prime);
B_prime_fft = fft(B_prime);

% 安全除法：添加正则项
epsilon = 1e-8;
H_fft = A_prime_fft ./ (A_fft + epsilon);
B_fft = B_prime_fft ./ (H_fft + epsilon);

% 逆FFT还原时域信号
B_recon = real(ifft(B_fft));

% 可视化对比
figure;
subplot(2,1,1);
plot(t, B_true, 'b-', t, B_recon, 'r--');
legend('真实输入B', '还原的B');
title('时域信号对比');
subplot(2,1,2);
freq = (0:N-1)*(fs/N);
plot(freq(1:N/2), abs(B_fft(1:N/2)));
xlim([0 200]);
title('频谱对比');

核心处理技巧 #

• 正则项epsilon：给分母加一个极小的正数，避免除以0或极小值导致的数值溢出，epsilon的大小可以根据信号噪声水平调整（噪声大可以适当调大，比如1e-6）。
• 实部提取：因为原信号是实信号，逆FFT后会有微小的虚部（数值计算误差导致），必须用real()提取实部得到正确的时域信号。
• 频谱范围限制：如果知道系统的有效工作频率范围，可以只在该范围内做除法，其他频率点直接置0，减少噪声和无效频谱的干扰。
• 相位校准：如果需要还原信号的准确相位，要确保H(w)的相位估计准确（比如用零相位滤波模拟系统，避免相位偏移）。

内容的提问来源于stack exchange，提问作者P.Alipoor