Swift实现4阶Butterworth带通滤波器遇滤波效果不佳问题求助

我需要实现一个截止频率0.05Hz至0.15Hz的4阶Butterworth带通滤波器，但当前代码滤波效果极差，仍有大量非目标频率信号通过，尝试级联两个二阶滤波器也未解决问题。以下是我的实现代码：

// 基于中心频率和品质因数计算带通滤波器系数
func makeBandpassFilterWithFcAndQ(filter: inout [Double], Fs: Double, Fc: Double, Q: Double) {
    let K = tan(.pi * Fc / Fs)
    let norm = 1.0 / (1.0 + K / Q + K * K)
    filter[0] = K / Q * norm
    filter[1] = 0.0
    filter[2] = -filter[0]
    filter[3] = 2.0 * (K * K - 1.0) * norm
    filter[4] = (1.0 - K / Q + K * K) * norm
}

// 基于频率范围计算带通滤波器系数
func makeBandpassFilterWithFreqRange(filter: inout [Double], Fs: Double, Fbtm: Double, Ftop: Double) {
    //let Fc = sqrt(Fbtm * Ftop)
    let Fc = 0.10
    let Q = Fc / (Ftop - Fbtm)
    //let Q = 0.99
    makeBandpassFilterWithFcAndQ(filter: &filter, Fs: Fs, Fc: Fc, Q: Q)
}

// 对信号应用带通滤波
func applyBandpassFilter(signal: [Double], fs: Double, filterCoefficients: [Double]) -> [Double] {
    let numRawDoubles = vDSP_Length(signal.count)
    var rawDoubles = signal
    var filteredDoubles = [Double](repeating: 0.0, count: signal.count)

    let rawStride = vDSP_Stride(1)
    let filteredStride = vDSP_Stride(1)

    // 使用vDSP_deq22应用滤波
    vDSP_deq22D(&rawDoubles, rawStride, filterCoefficients, &filteredDoubles, filteredStride, numRawDoubles - 2)

    return filteredDoubles
}

// 应用4阶Butterworth带通滤波
func apply4thOrderButterworthBandpassFilter(signal: [Double], fs: Double, Fbtm: Double, Ftop: Double) -> [Double] {
    // 第一组系数
    var filterCoefficients1 = [Double](repeating: 0.0, count: 5)
    makeBandpassFilterWithFreqRange(filter: &filterCoefficients1, Fs: fs, Fbtm: Fbtm, Ftop: Ftop)
    
    // 第一次滤波
    let filteredSignal1 = applyBandpassFilter(signal: signal, fs: fs, filterCoefficients: filterCoefficients1)
    
    // 第二组系数（和第一组相同）
    var filterCoefficients2 = [Double](repeating: 0.0, count: 5)
    makeBandpassFilterWithFreqRange(filter: &filterCoefficients2, Fs: fs, Fbtm: Fbtm, Ftop: Ftop)
    
    // 第二次滤波
    let filteredSignal2 = applyBandpassFilter(signal: filteredSignal1, fs: fs, filterCoefficients: filterCoefficients2)
    
    return filteredSignal2
}

问题根源分析 #

1. 中心频率计算错误
   你注释掉了正确的几何中心频率计算let Fc = sqrt(Fbtm * Ftop)，硬编码为0.10Hz。对于0.05-0.15Hz的通带，正确的几何中心频率应为sqrt(0.05*0.15) ≈ 0.0866Hz，硬编码会导致通带偏移，无法准确覆盖目标频率范围。

2. 4阶滤波器级联方式错误
   4阶Butterworth带通滤波器是由两个不同的二阶Butterworth带通节级联而成，而非两个完全相同的二阶节。每个二阶节的Q值不同，直接复用相同系数无法达到4阶Butterworth的滚降特性。

3. vDSP_deq22D未处理滤波状态
   vDSP_deq22D实现的是IIR滤波器的直接II型结构，滤波过程依赖前一次的输入和输出状态。当前代码每次调用都从零初始化状态，级联时会导致滤波连续性断裂，不仅影响开头的信号质量，还会削弱滤波效果。

修正方案 #

1. 恢复正确的中心频率计算 #

在makeBandpassFilterWithFreqRange中恢复几何中心频率的计算：

func makeBandpassFilterWithFreqRange(filter: inout [Double], Fs: Double, Fbtm: Double, Ftop: Double) {
    let Fc = sqrt(Fbtm * Ftop) // 恢复正确的几何中心频率
    let Q = Fc / (Ftop - Fbtm)
    makeBandpassFilterWithFcAndQ(filter: &filter, Fs: Fs, Fc: Fc, Q: Q)
}

2. 生成4阶Butterworth的两个二阶节系数 #

4阶Butterworth带通的两个二阶节对应不同的Q值：1/(2*sin(π/8)) ≈ 1.3066和1/(2*sin(3π/8)) ≈ 0.5412，基于中心频率和带宽分别计算两组系数。

3. 修正滤波函数，保留滤波器状态 #

修改滤波函数，添加状态参数以保留滤波过程中的延迟线值，确保级联时的连续性：

func applyBandpassFilter(signal: [Double], filterCoefficients: [Double], state: inout [Double]) -> [Double] {
    let numSamples = vDSP_Length(signal.count)
    var input = signal
    var output = [Double](repeating: 0.0, count: signal.count)
    
    vDSP_deq22D(&input, 1, filterCoefficients, &output, 1, numSamples, &state, 1)
    
    return output
}

完整修正后的4阶滤波函数 #

func apply4thOrderButterworthBandpassFilter(signal: [Double], fs: Double, Fbtm: Double, Ftop: Double) -> [Double] {
    let Fc = sqrt(Fbtm * Ftop)
    let bandwidth = Ftop - Fbtm
    
    // 4阶Butterworth带通的两个二阶节Q值
    let Q1 = 1.0 / (2 * sin(.pi / 8))
    let Q2 = 1.0 / (2 * sin(3 * .pi / 8))
    
    // 生成第一组二阶节系数
    var coeffs1 = [Double](repeating: 0.0, count: 5)
    makeBandpassFilterWithFcAndQ(filter: &coeffs1, Fs: fs, Fc: Fc, Q: Q1)
    
    // 生成第二组二阶节系数
    var coeffs2 = [Double](repeating: 0.0, count: 5)
    makeBandpassFilterWithFcAndQ(filter: &coeffs2, Fs: fs, Fc: Fc, Q: Q2)
    
    // 初始化两个滤波器的状态（每个二阶节需要2个状态值）
    var state1 = [Double](repeating: 0.0, count: 2)
    var state2 = [Double](repeating: 0.0, count: 2)
    
    // 级联滤波
    let filtered1 = applyBandpassFilter(signal: signal, filterCoefficients: coeffs1, state: &state1)
    let filtered2 = applyBandpassFilter(signal: filtered1, filterCoefficients: coeffs2, state: &state2)
    
    return filtered2
}

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