Python技术问询：如何保留基频倍数频率（FFT或filtfilt实现）

你用signal.iircomb时只看到单个频率点的响应，是因为代码里分别针对F0和2F0创建了独立的单峰/单陷波滤波器，而非真正的梳状滤波器。iircomb本身可以生成在基频所有倍数处都有一致响应的梳状特性，只要正确指定基频参数。

一、正确使用IIR梳状滤波器 #

直接基于基频生成梳状滤波器，就能在所有基频倍数处得到峰值响应，保留目标频率、滤除其他分量：

from scipy import signal
from matplotlib import pyplot as plt

BaseFreq = 90.0
fs = 36000
Q0 = 50.0

# 生成针对基频的峰值梳状滤波器（保留基频及其所有倍数）
b, a = signal.iircomb(BaseFreq, Q0, ftype='peak', fs=fs)
freq, h = signal.freqz(b, a, fs=fs)
response = abs(h)

# 绘制响应
plt.figure(figsize=(12,6))
plt.plot(freq, response)
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度响应')
plt.title('IIR峰值梳状滤波器响应（保留基频及其倍数）')
plt.grid(True)
plt.show()

说明：Q值控制峰值尖锐度，Q越大，峰值越窄，对非谐波频率的衰减效果越强。

二、频域FFT实现方案 #

通过FFT修改频谱，直接保留基频倍数分量，其余置零后逆变换：

import numpy as np

# 假设MySignal是你的时域信号
MySignal = np.random.randn(4096)  # 示例信号，替换为你的实际数据
n = len(MySignal)
freq_res = fs / n

# 做FFT并获取频率轴
fft_signal = np.fft.fft(MySignal)
freqs = np.fft.fftfreq(n, 1/fs)

# 保留基频及其倍数的频谱分量
filtered_fft = np.zeros_like(fft_signal)
max_harmonic = int(fs/(2*BaseFreq))  # 最高可保留的谐波频率（奈奎斯特频率内）
for k in range(1, max_harmonic + 1):
    target_freq = BaseFreq * k
    # 找到最接近目标频率的频谱索引
    idx = np.argmin(np.abs(freqs - target_freq))
    filtered_fft[idx] = fft_signal[idx]
    filtered_fft[-idx] = fft_signal[-idx]  # 处理实信号的负频率对称分量

# 逆FFT得到滤波后的时域信号
filtered_signal = np.fft.ifft(filtered_fft).real

注意：若信号长度不是基频周期的整数倍，会出现频谱泄漏，建议先加汉宁窗再做FFT，或调整信号长度为基频周期的整数倍。

三、Kaiser窗FIR梳状滤波器实现 #

用Kaiser窗截断理想梳状脉冲响应，得到有限长FIR滤波器，兼顾过渡带和阻带衰减：

import numpy as np
from scipy import signal
from matplotlib import pyplot as plt

BaseFreq = 90.0
fs = 36000
MySignal = np.random.randn(4096)  # 示例信号

N = round(fs / BaseFreq)  # 基频对应的采样周期数
taps = 10 * N  # 滤波器长度，越长越接近理想梳状特性

# 生成理想梳状脉冲响应（周期冲激串）
ideal_impulse = np.zeros(taps)
ideal_impulse[::N] = 1

# 用Kaiser窗截断，beta控制阻带衰减（beta=14对应约90dB衰减）
beta = 14
kaiser_window = np.kaiser(taps, beta)
fir_coeffs = ideal_impulse * kaiser_window

# 归一化滤波器系数
fir_coeffs /= np.sum(fir_coeffs)

# 查看频率响应
freq, h = signal.freqz(fir_coeffs, 1, fs=fs)
plt.figure(figsize=(12,6))
plt.plot(freq, abs(h))
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度响应')
plt.title('Kaiser窗FIR梳状滤波器响应')
plt.grid(True)
plt.show()

# 对信号滤波
filtered_signal = signal.lfilter(fir_coeffs, 1, MySignal)

说明：beta参数可按需调整，beta越大阻带衰减越强；滤波器长度越长，谐波处的响应越接近1，过渡带越窄。

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