已知频率高斯峰滤波：Butterworth滤波器异常问题排查

我有一组间隔已知的高斯峰数据，叠加在低频噪声之上，每个高斯峰的高度服从Poisson分布，峰高各不相同。尝试用Butterworth高通滤波器去除低频噪声，但结果出现三个异常：

• 信号低于0；
• 峰未处于平坦基线，推测和峰幅差异有关；
• 目标峰左右出现“幽灵峰”。

附上实现代码，请求排查错误：

from scipy.stats import norm
from scipy.stats import skewnorm
from scipy.stats import poisson
from astropy.stats import freedman_bin_width
from scipy.interpolate import interp1d

Delta_X = 5
gaus_width_0 = 0.5
gaus_width_1 = 0.02
mu=2


######
# DEFINE HIGH PASS BUTTERWORTH FILTER 
######

def butter_highpass(cutoff, fs, order=5):
   nyq = 0.5 * fs
   normal_cutoff = cutoff / nyq
   b, a = signal.butter(order, normal_cutoff, btype='high', analog=False)
   return b, a

def butter_highpass_filter(data, cutoff, fs, order=5):
   b, a = butter_highpass(cutoff, fs, order=order)
   y = signal.filtfilt(b, a, data)
   return y


######
# GENERATE NOISE - GAUSSIANS SEPARATED BY Delta_X, WITH SOME SMEARING. EACH PEAK HAS A POISSON PROBABILITY
######

data = []
for i in range(10):
  mean = float(i)*Delta_X
  std = gaus_width_0 + gaus_width_1*i
  gaus = norm.rvs(mean, std, size=int(20000*poisson.pmf(i, mu)))
  data.append(gaus)

######
# GENERATE NOISE - A SKEWED NORMAL DISTRIBUTION
######

noise = 4*Delta_X*skewnorm.rvs(5, size=100000)
data.append(noise)

data = np.concatenate(data)


######
# GENERATE A HISTOGRAM
####

bw, bins = freedman_bin_width(data, return_bins=True)
counts, _ = np.histogram(data, bins=bins)
bin_centres =  (bins[:-1] + bins[1:])/2.
nbins = len(bin_centres)
bin_numbers = np.arange(0, nbins)

######
# INTERPOLATE BETWEEN BIN UNITS AND MEASURED UNITS TO FIND THE FREQUENCY, IN BINS
####


f_inv = interp1d(bin_centres, bin_numbers)
freq_bins = f_inv(Delta_X) - f_inv(0)


######
# APPLY THE BUTTERWORTH FILTER
####

filtered_counts = butter_highpass_filter(counts, freq_bins, nbins)

plt.figure(figsize=(10,10))
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.xlabel("x", fontsize=20)
plt.plot(bin_centres, counts, lw=3, label="Original Histogram")
plt.plot(bin_centres, filtered_counts, lw=3, label="Original Histogram")

plt.grid(which="both")
plt.show()

运行后得到异常结果：

核心错误分析 #

1. 采样频率与截止频率的概念完全混淆 #

你把峰间隔对应的bin数当作截止频率，把总bin数当作采样频率，这不符合数字滤波器的频率定义：

• 数字滤波器的采样频率fs是单位x轴区间内的采样点数量，即1/bw（bw是直方图的bin宽度），而非总bin数nbins。
• 截止频率应设为需要滤除的低频成分的频率：你的峰间隔是5，对应特征频率为1/5=0.2，截止频率要小于这个值，才能保留峰的成分。

2. 未去除直流分量 #

原始直方图带有正的直流基线，高通滤波器会滤除该直流分量，直接导致信号出现负值；同时低频噪声与直流分量混合，会让滤波后的基线无法平坦。

3. 振铃效应（幽灵峰）的来源 #

• 滤波器阶数过高（默认5阶）会让频率响应过渡带过窄，结合错误的截止频率，引发振铃效应。
• 直方图边缘数据不连续，滤波时的边缘效应会放大振铃。

修正后的代码 #

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.stats import skewnorm
from scipy.stats import poisson
from astropy.stats import freedman_bin_width
from scipy.interpolate import interp1d

Delta_X = 5
gaus_width_0 = 0.5
gaus_width_1 = 0.02
mu=2


######
# 定义Butterworth高通滤波器（降低阶数减少振铃）
######

def butter_highpass(cutoff, fs, order=3):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = signal.butter(order, normal_cutoff, btype='high', analog=False)
    return b, a

def butter_highpass_filter(data, cutoff, fs, order=3):
    b, a = butter_highpass(cutoff, fs, order=order)
    # 先去除直流分量，避免滤波后出现负值
    data_dc_removed = data - np.mean(data)
    y = signal.filtfilt(b, a, data_dc_removed)
    return y


######
# 生成带Poisson峰高的高斯峰数据
######

data = []
for i in range(10):
    mean = float(i)*Delta_X
    std = gaus_width_0 + gaus_width_1*i
    gaus = norm.rvs(mean, std, size=int(20000*poisson.pmf(i, mu)))
    data.append(gaus)

######
# 生成低频偏态噪声
######

noise = 4*Delta_X*skewnorm.rvs(5, size=100000)
data.append(noise)
data = np.concatenate(data)


######
# 生成直方图
####

bw, bins = freedman_bin_width(data, return_bins=True)
counts, _ = np.histogram(data, bins=bins)
bin_centres =  (bins[:-1] + bins[1:])/2.

######
# 修正频率参数：采样频率为单位x的bin数，截止频率低于峰特征频率
####

fs = 1 / bw  # 每单位x轴长度对应的采样点数量
cutoff_freq = 0.1  # 峰特征频率为0.2，截止频率设为0.1确保保留峰

######
# 应用滤波器
####

filtered_counts = butter_highpass_filter(counts, cutoff_freq, fs)
# 可选：若需要非负结果，可加上原始数据的最小正值恢复基线
# filtered_counts = filtered_counts + np.min(counts[counts>0])

plt.figure(figsize=(10,10))
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.xlabel("x", fontsize=20)
plt.plot(bin_centres, counts, lw=3, label="原始直方图")
plt.plot(bin_centres, filtered_counts, lw=3, label="滤波后结果")
plt.legend(fontsize=16)
plt.grid(which="both")
plt.show()

关键修正说明 #

• 补全缺失的依赖包导入：numpy、scipy.signal、matplotlib.pyplot，否则代码无法运行。
• 修正频率计算逻辑：采样频率fs设为1/bw，截止频率设为小于峰特征频率的0.1，确保峰成分被保留。
• 添加直流分量去除：滤波前减去数据均值，避免出现负值，同时让基线更平坦。
• 降低滤波器阶数：从5阶改为3阶，减少振铃效应。
• 可选基线恢复：若需要非负信号，可在滤波后加上原始数据的最小正值，贴合直方图的物理意义。

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