在R中构建匹配scipy.signal.firwin的多带FIR阻带核及问题排查

我正在用R从零实现基于firwin/Hamming窗的FIR滤波，单阻带场景下的实现和scipy.signal.firwin完全匹配，但将多个阻带合并为单个核时遇到了问题。

多陷波实现步骤 #

以50、100、150Hz陷波为例，采样率sfreq=1024，过渡带宽trans_bandwidth=1Hz：

步骤1：计算每个陷波的边缘 #

参考Python代码：

tb_2  = 0.5  # 半过渡带宽
freqs = [50, 100, 150]
nw    = [f/200 for f in freqs]   # 陷波宽度 = [0.25, 0.5, 0.75]
lows  = [f - n/2 - tb_2 for f,n in zip(freqs,nw)]  # [49.375, 99.5,  149.625]
highs = [f + n/2 + tb_2 for f,n in zip(freqs,nw)]  # [50.625, 100.5, 150.375]
f_p1 = [l + tb_2 for l in lows]   # 阻带内边缘: [49.875, 100.0, 150.125]
f_p2 = [h - tb_2 for h in highs]  # 阻带内边缘: [50.125, 100.0, 149.875]

步骤2：构建合并后的排序频率/增益数组 #

freq = sorted(f_p1 + lows + highs + f_p2)
gain = # lows/highs对应通带边缘设为1，f_p1/f_p2对应阻带边缘设为0
# 前置0，后置奈奎斯特频率（512Hz）

步骤3：用firwin方法构建单个核 #

通过从右到左扫描，将sinc×Hamming子核嵌入长度为N的主核。

遇到的问题 #

注意100Hz陷波的f_p1[1] = 100.0且f_p2[1] = 100.0，频率数组中出现重复值，导致零宽过渡带，进而在核构建函数中触发除以零错误（3.3 / transition，其中transition = 0）。

核构建函数代码：

.firwin_kernel <- function(N, freq, gain) {
  h      <- numeric(N)
  center <- N %/% 2L + 1L
  if (gain[length(gain)] == 1) h[center] <- 1.0
  prev_freq <- freq[length(freq)]
  prev_gain <- gain[length(gain)]
  for (i in seq(length(freq) - 1L, 1L)) {
    this_freq <- freq[i]
    this_gain <- gain[i]
    if (this_gain != prev_gain) {
      transition <- (prev_freq - this_freq) / 2.0
      this_N     <- as.integer(round(3.3 / transition))  # transition=0时崩溃
      this_N     <- this_N + 1L - this_N %% 2L
      cutoff     <- (prev_freq + this_freq) / 2.0
      this_h     <- .firwin_lowpass(this_N, cutoff)
      offset     <- (N - this_N) %/% 2L
      idx        <- seq(offset + 1L, offset + this_N)
      if (this_gain == 0L) h[idx] <- h[idx] - this_h
      else                 h[idx] <- h[idx] + this_h
    }
    prev_gain <- this_gain
    prev_freq <- this_freq
  }
  h
}

疑问与解答 #

疑问1：相邻阻带产生重复频率值时该如何处理？ #

1. 根源排查：100Hz陷波的阻带内边缘重合，本质是该陷波的阻带宽度为0，这在实际中是不合理的——FIR滤波器无法实现无限窄的阻带。建议先调整陷波宽度nw或过渡带宽，确保阻带宽度大于0。
2. 频率点去重与合并：如果确实需要处理相邻/重叠频带，构建freq数组时先对所有频率点去重再排序，同时重新计算对应gain值，确保相邻频率点之间存在非零过渡带。
3. 代码保护机制：在核构建函数中加入判断，当transition == 0时跳过该过渡带的处理（零宽度过渡带无需构建子核），避免除以零错误：

   if (transition <= 0) {
     prev_gain <- this_gain
     prev_freq <- this_freq
     next
   }
   

疑问2：从右到左扫描嵌入子核的方法在多带场景下是否仍有效？ #

这个方法（firwin的叠加法，从全通滤波器开始，通过加减低通子核构建多带滤波器）在多带场景下是有效的，但需要满足：

• 频率点严格递增，且相邻频率点之间的过渡带宽度非零；
• 当频带接近但未重叠时，只要过渡带宽度正常，就不会失效；
• 若出现重复频率点（零宽过渡带），则会触发错误，需要提前处理重复点或在代码中加入保护逻辑；
• 额外建议：加入子核长度this_N不超过主核N的判断，避免索引越界问题。

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