如何在 R 中构造与 scipy.signal.firwin 一致的多带阻 FIR 核？

先核对边界计算，再检查频率与增益是否配对。给出的表达式不会在 100 Hz 处生成重复边界。即使两个相邻点频率相同，只要它们的 gain 都是 0，当前代码也不会进入除法分支。从右向左叠加子核可以用于多阻带，但按不同长度生成子核再嵌入，通常不能保证得到与一次调用 scipy.signal.firwin 相同的系数。

先核对边界：代码注释中的数值有误

按给出的表达式计算，实际结果如下，单位均为 Hz：

100 Hz 陷波的阻带边界是 99.75 和 100.25，没有重复；150 Hz 陷波的边界也没有倒置。这三个陷波相隔较远，相邻阻带没有相接。

可以用 R 重算，不要把注释里的数值当成实际输入：

sfreq <- 1024
freqs <- c(50, 100, 150)
tb_2 <- 0.5
nw <- freqs / 200

lows <- freqs - nw / 2 - tb_2
highs <- freqs + nw / 2 + tb_2
f_p1 <- lows + tb_2
f_p2 <- highs - tb_2

print(data.frame(freqs, nw, lows, f_p1, f_p2, highs))

如果运行时确实出现了两个 100.0，需要追踪实际数组的生成过程，检查是否存在舍入、手工录入或其他未展示的计算。

重复频率要结合 gain 判断

原函数只在下面的条件成立时计算 transition：

if (this_gain != prev_gain) {

两个相邻的 (100, 0) 不满足该条件，因此不会直接导致 3.3 / 0。实际除零意味着参与计算的两个点频率相同、增益不同，或者运行代码与展示代码存在差异。

处理方式取决于重复点的含义：

• 相同频率、相同增益：可以合并为一个点，同时保留频率与增益的配对。
• 相同频率、不同增益：表示零宽度跳变，当前依赖正过渡带宽的长度公式无法处理。应检查边界构造，或重新定义过渡区。
• 两个阻带真正重叠或相接：如果设计意图是在整个并集上阻带增益为 0，可以合并阻带，并重新生成外侧过渡边界。
• 两个阻带之间仍要求通带：不能直接合并，否则会删除该通带。需要保留通带宽度，并检查所选滤波核长度能否满足响应要求。

将频率和增益组成配对记录后一起排序，不要只对 freq 排序，再按原顺序填写 gain：

edges <- data.frame(
  freq = c(0, lows, f_p1, f_p2, highs, sfreq / 2),
  gain = c(1, rep(1, length(lows)),
           rep(0, length(f_p1)), rep(0, length(f_p2)),
           rep(1, length(highs)), 1)
)
edges <- edges[order(edges$freq), , drop = FALSE]

这个构造适用于当前互不重叠的三个陷波。阻带重叠时，只排序无法解决通带与阻带要求之间的冲突，需要先整理频带定义。

给现有扫描算法增加边界检查

在计算长度之前拒绝无效的过渡宽度，可以直接定位问题，避免错误传到索引计算：

if (this_gain != prev_gain) {
  delta <- prev_freq - this_freq
  if (!is.finite(delta) || delta <= 0) {
    stop("Gain changes require a positive transition width")
  }

  transition <- delta / 2
  this_N <- round(3.3 / transition)
  this_N <- this_N + 1 - this_N %% 2

  if (!is.finite(this_N) || this_N > N) {
    stop("Master kernel is too short for this transition")
  }
  this_N <- as.integer(this_N)

  # Continue with cutoff calculation and centered embedding.
}

频率单位也要统一。如果传入的 freq 已除以 sfreq / 2，那么 delta / 2 等于以采样率归一化的过渡宽度；如果传入的是 Hz，就不能直接套用同一个 3.3 / transition 表达式。这个经验长度公式本身也不保证指定的通带纹波或阻带衰减。

当前居中嵌入方式要求 N 和 this_N 都为奇数，且 this_N <= N。不要添加一个很小的数来绕过零宽度，这会产生极大的子核长度，也会掩盖输入问题。

为什么扫描叠加不一定等于 firwin

扫描方向本身没有问题。理想多带滤波器可以表示为若干低通响应的加减，因此可以从右向左处理增益变化。

差别在于子核的构造方式。firwin 在同一个长度为 N 的采样网格上累加各通带的 sinc 项，再统一乘窗并进行整体缩放，不会为每条过渡边界分别计算一个 this_N。这一点可在 SciPy 的 firwin 实现中核对。

不同长度的 Hamming 子核即使居中嵌入，窗形状也不同。如果 .firwin_lowpass() 还对各个子核分别归一化，结果也可能与最后只归一化一次不同。单阻带时一致，并不能说明多阻带时也一定一致。

按统一长度构造多带阻核

如果目标是对齐一次 firwin 调用，可以直接指定严格递增的截止频率。firwin 的 cutoff 表示半幅值截止位置，并非完整的通带、阻带边界数组，且不能包含 0、Nyquist 或重复值。当前带阻设计在 Nyquist 处保留通带，因此需要奇数长度。参数约束见 SciPy 文档。

若沿用原扫描算法的约定，以过渡区中点作为截止位置，六个截止频率为：

cutoff_hz <- as.vector(rbind(
  (lows + f_p1) / 2,
  (f_p2 + highs) / 2
))
# 49.625 50.375 99.5 100.5 149.375 150.625

下面的 R 函数适用于 DC 和 Nyquist 均为通带、内部有若干阻带的情况。它使用统一长度、对称 Hamming window 和整体 DC 归一化：

firwin_multistop <- function(N, cutoff_hz, sfreq, scale = TRUE) {
  stopifnot(
    length(N) == 1L, is.finite(N),
    N >= 3, N == floor(N), N %% 2 == 1,
    length(sfreq) == 1L, is.finite(sfreq), sfreq > 0,
    length(cutoff_hz) >= 2L,
    length(cutoff_hz) %% 2L == 0L,
    all(is.finite(cutoff_hz)),
    all(cutoff_hz > 0), all(cutoff_hz < sfreq / 2),
    all(diff(cutoff_hz) > 0)
  )

  sinc <- function(x) {
    y <- rep(1, length(x))
    nz <- x != 0
    y[nz] <- sin(pi * x[nz]) / (pi * x[nz])
    y
  }

  m <- (seq_len(N) - 1) - (N - 1) / 2
  bounds <- c(0, cutoff_hz / (sfreq / 2), 1)
  h <- numeric(N)

  for (j in seq.int(1L, length(bounds) - 1L, by = 2L)) {
    left <- bounds[j]
    right <- bounds[j + 1L]
    h <- h + right * sinc(right * m) - left * sinc(left * m)
  }

  w <- 0.54 - 0.46 * cos(2 * pi * (seq_len(N) - 1) / (N - 1))
  h <- h * w
  if (scale) h <- h / sum(h)
  h
}

对应的 Python 参考调用是：

from scipy import signal

# N must be the same odd number used in R.
cutoff_hz = [49.625, 50.375, 99.5, 100.5, 149.375, 150.625]
h_ref = signal.firwin(
    numtaps=N,
    cutoff=cutoff_hz,
    fs=1024,
    window="hamming",
    pass_zero="bandstop",
    scale=True,
)

上述实现采用统一网格的 sinc 累加方式。这里没有给出实际跨语言运行结果，不能承诺系数逐位相同。比较时应统一 N、截止频率、对称窗和 scale，再检查系数最大绝对差及频率响应。即使系数对齐，也仍需确认所选长度能否达到所需的陷波深度和过渡带要求。

备注：内容仅供参考。