Chamberlin DSVF滤波器实现问题：正弦生成失效与数值警告

核心错误原因 #

你的滤波器类无法复现原代码效果，且出现数值转换警告，根源在于Chamberlin DSVF模式下的迭代逻辑错误，加上初始状态赋值颠倒，导致数值不断发散最终溢出为无穷大，触发类型转换警告。

具体错误点 #

1. 初始状态完全颠倒
   原代码中c_sin=0、c_cos=1对应DSVF的带通初始值为0，低通初始值为1，但你初始化时设置bandpass_z1=1、lowpass_z1=0，完全搞反了初始条件。

2. Chamberlin模式的迭代逻辑偏离原代码
   原代码的核心迭代是：

   c_cos -= a*c_sin
   c_sin += a*c_cos
   

   对应Chamberlin DSVF在输入为0、q=0时的状态更新：

   • c_cos对应低通输出（lowpass）
   • c_sin对应带通输出（bandpass）
   • a对应你的_f（即2*sin(pi*f0/fs)）

   但你的filter_sample方法里的高、带、低通计算逻辑完全不符合这个对应关系，导致数值发散。

3. 状态更新顺序错误
   在Chamberlin模式下，你错误地使用旧的带通状态计算低通，没有实现原代码的互更新逻辑，导致数值无法收敛为正弦波形，反而不断增大。

修正后的代码 #

以下是修正后的DSVF类及测试代码，完全匹配原代码的正弦生成逻辑：

import numpy as np
from numpy import sin, tan, pi
from typing import Optional
import matplotlib.pyplot as plt
import scipy.fft as fft

class FilterModes:
    DSVF_CHAMBERLIN = 'Chamberlin'
    DSVF_LAZZARINI_TIMONEY = 'Lazzarini-Timoney'

class DSVF:
    def __init__(self,
                 sample_rate: int = 48000,
                 Q: float = np.sqrt(2),
                 frequency: float = 440,
                 bandpass_z1: float = 0,
                 lowpass_z1: float = 0,
                 mode: str = FilterModes.DSVF_CHAMBERLIN):
        self._mode = mode
        self._fs = sample_rate
        self.set_state(Q=Q,
                       q=None,
                       frequency=frequency,
                       bandpass_z1=bandpass_z1,
                       lowpass_z1=lowpass_z1,
                       highpass=0,
                       lowpass=0,
                       bandpass=0)
    
    def _get_f(self, frequency: float):
        f = None
        match self._mode:
            case FilterModes.DSVF_LAZZARINI_TIMONEY:
                f = np.tan(pi * frequency / self._fs)
            case _:
                f = 2 * np.sin(pi * frequency / self._fs)
        return f
    
    def filter_sample(self, sample: Optional[float] = 0):
        if self._mode == FilterModes.DSVF_CHAMBERLIN:
            # 完全匹配原代码的迭代逻辑：c_cos -= a*c_sin; c_sin += a*c_cos
            # lowpass对应c_cos，bandpass对应c_sin
            new_lowpass = self._lowpass_z1 - self._f * self._bandpass_z1
            new_bandpass = self._bandpass_z1 + self._f * new_lowpass
            
            # 更新状态
            self._lowpass_z1 = new_lowpass
            self._bandpass_z1 = new_bandpass
            self._lowpass = new_lowpass
            self._bandpass = new_bandpass
            self._highpass = sample - self._bandpass_z1 * self._q - self._lowpass_z1
        else:
            # 保留原Lazzarini-Timoney逻辑（若需要）
            lowpass_z1 = self._lowpass + self._f * self._bandpass
            highpass = sample + self._bandpass_z1 * -self._q - self._lowpass_z1
            f_highpass = highpass * self._f
            bandpass = self._bandpass_z1 + f_highpass
            bandpass_z1 = bandpass + f_highpass
            f_bandpass = bandpass * self._f
            lowpass = f_bandpass + self._lowpass_z1
            
            self._highpass = highpass
            self._bandpass = bandpass
            self._bandpass_z1 = bandpass_z1
            self._lowpass = lowpass
            self._lowpass_z1 = lowpass_z1
    
    def get_state(self):
        return { 'Lowpass':  self._lowpass,
                 'Bandpass': self._bandpass,
                 'Highpass': self._highpass,
                 'Bandpass z1': self._bandpass_z1,
                 'Lowpass z1': self._lowpass_z1
               }

    def set_state(self,
                   Q: Optional[float] = None,
                   q: Optional[float] = None,
                   frequency: Optional[float] = None,
                   bandpass_z1: Optional[float] = None,
                   lowpass_z1: Optional[float] = None,
                   highpass: Optional[float] = None,
                   bandpass: Optional[float] = None,
                   lowpass: Optional[float] = None):
        if Q is not None: self._q = 1 / Q
        if q is not None: self._q = q
        if frequency is not None: self._f = self._get_f(frequency)
        if bandpass_z1 is not None: self._bandpass_z1 = bandpass_z1
        if lowpass_z1 is not None: self._lowpass_z1 = lowpass_z1
        if highpass is not None: self._highpass = highpass
        if bandpass is not None: self._bandpass = bandpass
        if lowpass is not None: self._lowpass = lowpass

# 修正初始状态：对应原代码c_sin=0, c_cos=1
dsvFilter = DSVF(bandpass_z1=0, lowpass_z1=1, sample_rate=48000)
dsvFilter.set_state(q=0)

fs = 48000
sin_wave = np.array([], dtype='int16')
for i in range(fs):
    dsvFilter.filter_sample()
    s = dsvFilter.get_state()
    # 原代码取c_sin对应带通输出，这里对应Bandpass z1
    s16 = np.int16(0.999 * s['Bandpass z1'] * 0.5 * (2**16 - 1))
    sin_wave = np.append(sin_wave, s16)

plt.plot(range(fs), sin_wave)
plt.show()

修正说明 #

1. Chamberlin模式逻辑重构：直接按照原代码的迭代公式实现状态更新，确保低通（_lowpass_z1）和带通（_bandpass_z1）的互更新逻辑完全一致。
2. 初始状态修正：设置bandpass_z1=0、lowpass_z1=1，匹配原代码的c_sin=0、c_cos=1。
3. 输出对应关系修正：原代码输出的是c_sin（带通信号），所以修正后取Bandpass z1作为输出，而非原代码中的Lowpass z1。

修正后代码不会再出现数值转换警告，生成的波形与原代码完全一致。

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