跳转至

数字滤波器

学习目标

完成本模块后,你将能够: - 理解数字滤波器的基本原理和分类 - 掌握FIR滤波器的设计和实现 - 掌握IIR滤波器的设计和实现 - 选择合适的滤波器类型和参数 - 在嵌入式系统中高效实现数字滤波器 - 在医疗器械中应用数字滤波技术

前置知识

  • 信号处理基础
  • 数字信号处理理论
  • C语言编程
  • 定点数运算
  • ADC/DAC原理

内容

概念介绍

数字滤波器是医疗器械信号处理的核心技术,用于去除噪声、提取特征信号、抑制干扰。与模拟滤波器相比,数字滤波器具有精度高、稳定性好、易于调整等优点。

数字滤波器的应用: - 心电信号处理:去除工频干扰、基线漂移 - 血氧信号处理:去除运动伪影、环境光干扰 - 血压信号处理:平滑脉搏波形 - 脑电信号处理:提取特定频段的脑电波 - 超声信号处理:去除噪声、增强对比度

数字滤波器分类

按脉冲响应分类

  1. FIR滤波器(Finite Impulse Response)
  2. 有限脉冲响应
  3. 线性相位特性
  4. 稳定性好
  5. 计算量较大

  6. IIR滤波器(Infinite Impulse Response)

  7. 无限脉冲响应
  8. 非线性相位
  9. 计算效率高
  10. 可能不稳定

按频率特性分类

  • 低通滤波器(LPF):通过低频,抑制高频
  • 高通滤波器(HPF):通过高频,抑制低频
  • 带通滤波器(BPF):通过特定频段
  • 带阻滤波器(BSF):抑制特定频段

FIR滤波器

FIR滤波器原理

FIR滤波器的差分方程:

y[n] = b₀×x[n] + b₁×x[n-1] + b₂×x[n-2] + ... + bₙ×x[n-N]

其中: - y[n]:输出信号 - x[n]:输入信号 - b₀, b₁, ..., bₙ:滤波器系数 - N:滤波器阶数

FIR滤波器实现

// FIR滤波器结构
typedef struct {
    float* coeffs;        // 滤波器系数
    float* buffer;        // 延迟线缓冲区
    uint16_t length;      // 滤波器长度(阶数+1)
    uint16_t index;       // 当前索引
} FIR_Filter_t;

// 初始化FIR滤波器
void FIR_Init(FIR_Filter_t* fir, float* coeffs, uint16_t length) {
    fir->coeffs = coeffs;
    fir->length = length;
    fir->index = 0;

    // 分配并清零缓冲区
    fir->buffer = (float*)malloc(length * sizeof(float));
    memset(fir->buffer, 0, length * sizeof(float));
}

// FIR滤波器处理(浮点实现)
float FIR_Process(FIR_Filter_t* fir, float input) {
    float output = 0.0f;

    // 将新样本存入缓冲区
    fir->buffer[fir->index] = input;

    // 计算输出(卷积)
    uint16_t buf_idx = fir->index;
    for (uint16_t i = 0; i < fir->length; i++) {
        output += fir->coeffs[i] * fir->buffer[buf_idx];

        // 循环索引
        if (buf_idx == 0) {
            buf_idx = fir->length - 1;
        } else {
            buf_idx--;
        }
    }

    // 更新索引
    fir->index++;
    if (fir->index >= fir->length) {
        fir->index = 0;
    }

    return output;
}

// 释放FIR滤波器
void FIR_Free(FIR_Filter_t* fir) {
    if (fir->buffer != NULL) {
        free(fir->buffer);
        fir->buffer = NULL;
    }
}

代码说明: - 使用循环缓冲区存储历史样本 - 通过卷积计算输出 - 浮点实现,精度高但计算量大

FIR滤波器定点实现

// 定点FIR滤波器(Q15格式)
typedef struct {
    int16_t* coeffs;      // 滤波器系数(Q15)
    int16_t* buffer;      // 延迟线缓冲区(Q15)
    uint16_t length;      // 滤波器长度
    uint16_t index;       // 当前索引
} FIR_Filter_Q15_t;

// 初始化定点FIR滤波器
void FIR_Q15_Init(FIR_Filter_Q15_t* fir, int16_t* coeffs, uint16_t length) {
    fir->coeffs = coeffs;
    fir->length = length;
    fir->index = 0;

    fir->buffer = (int16_t*)malloc(length * sizeof(int16_t));
    memset(fir->buffer, 0, length * sizeof(int16_t));
}

// 定点FIR滤波器处理
int16_t FIR_Q15_Process(FIR_Filter_Q15_t* fir, int16_t input) {
    int32_t acc = 0;  // 32位累加器

    // 存入新样本
    fir->buffer[fir->index] = input;

    // 计算输出
    uint16_t buf_idx = fir->index;
    for (uint16_t i = 0; i < fir->length; i++) {
        // Q15 × Q15 = Q30,需要右移15位得到Q15
        acc += ((int32_t)fir->coeffs[i] * (int32_t)fir->buffer[buf_idx]);

        if (buf_idx == 0) {
            buf_idx = fir->length - 1;
        } else {
            buf_idx--;
        }
    }

    // 右移15位,从Q30转换为Q15
    int16_t output = (int16_t)(acc >> 15);

    // 更新索引
    fir->index++;
    if (fir->index >= fir->length) {
        fir->index = 0;
    }

    return output;
}

代码说明: - 使用Q15定点格式(16位有符号整数) - Q15表示范围:-1.0 到 +0.999969482421875 - 乘法结果为Q30,需要右移15位 - 定点运算速度快,适合嵌入式系统

低通FIR滤波器设计示例

// 50Hz低通滤波器(采样率1000Hz,31阶)
// 使用窗函数法设计
const float lpf_50hz_coeffs[31] = {
    -0.0018f, -0.0027f, -0.0025f, -0.0006f,  0.0032f,
     0.0084f,  0.0138f,  0.0176f,  0.0181f,  0.0143f,
     0.0062f, -0.0051f, -0.0171f, -0.0272f, -0.0327f,
    -0.0314f, -0.0220f, -0.0050f,  0.0171f,  0.0416f,
     0.0649f,  0.0832f,  0.0935f,  0.0935f,  0.0832f,
     0.0649f,  0.0416f,  0.0171f, -0.0050f, -0.0220f,
    -0.0314f
};

// 创建50Hz低通滤波器
FIR_Filter_t lpf_50hz;

void Init_LPF_50Hz(void) {
    FIR_Init(&lpf_50hz, (float*)lpf_50hz_coeffs, 31);
}

// 使用示例
float filtered_signal = FIR_Process(&lpf_50hz, raw_signal);

代码说明: - 截止频率50Hz,采样率1000Hz - 31阶滤波器,过渡带较窄 - 系数通过MATLAB或Python设计工具生成

IIR滤波器

IIR滤波器原理

IIR滤波器的差分方程:

y[n] = b₀×x[n] + b₁×x[n-1] + ... + bₘ×x[n-M]
       - a₁×y[n-1] - a₂×y[n-2] - ... - aₙ×y[n-N]

其中: - y[n]:输出信号 - x[n]:输入信号 - b₀, b₁, ..., bₘ:前向系数 - a₁, a₂, ..., aₙ:反馈系数

二阶IIR滤波器(Biquad)

二阶IIR滤波器是最常用的IIR滤波器结构,多个二阶节级联可以实现高阶滤波器。

// 二阶IIR滤波器(Biquad)结构
typedef struct {
    float b0, b1, b2;     // 前向系数
    float a1, a2;         // 反馈系数(a0=1)
    float x1, x2;         // 输入历史
    float y1, y2;         // 输出历史
} Biquad_Filter_t;

// 初始化Biquad滤波器
void Biquad_Init(Biquad_Filter_t* bq, float b0, float b1, float b2, 
                 float a1, float a2) {
    bq->b0 = b0;
    bq->b1 = b1;
    bq->b2 = b2;
    bq->a1 = a1;
    bq->a2 = a2;

    // 清零历史
    bq->x1 = 0.0f;
    bq->x2 = 0.0f;
    bq->y1 = 0.0f;
    bq->y2 = 0.0f;
}

// Biquad滤波器处理(Direct Form I)
float Biquad_Process(Biquad_Filter_t* bq, float input) {
    // 计算输出
    float output = bq->b0 * input + bq->b1 * bq->x1 + bq->b2 * bq->x2
                   - bq->a1 * bq->y1 - bq->a2 * bq->y2;

    // 更新历史
    bq->x2 = bq->x1;
    bq->x1 = input;
    bq->y2 = bq->y1;
    bq->y1 = output;

    return output;
}

// Biquad滤波器处理(Direct Form II,更节省内存)
float Biquad_Process_DF2(Biquad_Filter_t* bq, float input) {
    // Direct Form II只需要2个状态变量
    float w = input - bq->a1 * bq->x1 - bq->a2 * bq->x2;
    float output = bq->b0 * w + bq->b1 * bq->x1 + bq->b2 * bq->x2;

    // 更新状态
    bq->x2 = bq->x1;
    bq->x1 = w;

    return output;
}

代码说明: - Direct Form I:直观,但需要4个状态变量 - Direct Form II:节省内存,只需2个状态变量 - 二阶节稳定性好,易于实现

巴特沃斯低通滤波器

// 巴特沃斯低通滤波器(2阶,截止频率50Hz,采样率1000Hz)
void Init_Butterworth_LPF_50Hz(Biquad_Filter_t* bq) {
    // 系数通过设计工具计算
    float b0 = 0.0201f;
    float b1 = 0.0402f;
    float b2 = 0.0201f;
    float a1 = -1.5610f;
    float a2 = 0.6414f;

    Biquad_Init(bq, b0, b1, b2, a1, a2);
}

// 使用示例
Biquad_Filter_t butterworth_lpf;
Init_Butterworth_LPF_50Hz(&butterworth_lpf);

float filtered = Biquad_Process(&butterworth_lpf, raw_signal);

代码说明: - 巴特沃斯滤波器:通带平坦,过渡带平滑 - 2阶滤波器计算量小,适合实时处理 - 系数通过MATLAB或Python计算

陷波滤波器(Notch Filter)

陷波滤波器用于抑制特定频率的干扰,如50Hz/60Hz工频干扰。

// 50Hz陷波滤波器(采样率1000Hz,Q=30)
void Init_Notch_50Hz(Biquad_Filter_t* bq) {
    float fs = 1000.0f;   // 采样率
    float f0 = 50.0f;     // 陷波频率
    float Q = 30.0f;      // 品质因数(Q越大,陷波越窄)

    // 计算中间变量
    float w0 = 2.0f * M_PI * f0 / fs;
    float alpha = sinf(w0) / (2.0f * Q);

    // 计算系数
    float b0 = 1.0f;
    float b1 = -2.0f * cosf(w0);
    float b2 = 1.0f;
    float a0 = 1.0f + alpha;
    float a1 = -2.0f * cosf(w0);
    float a2 = 1.0f - alpha;

    // 归一化
    b0 /= a0;
    b1 /= a0;
    b2 /= a0;
    a1 /= a0;
    a2 /= a0;

    Biquad_Init(bq, b0, b1, b2, a1, a2);
}

// 使用示例:去除心电信号中的50Hz工频干扰
Biquad_Filter_t notch_50hz;
Init_Notch_50Hz(&notch_50hz);

float ecg_filtered = Biquad_Process(&notch_50hz, ecg_raw);

代码说明: - 陷波滤波器在特定频率处衰减最大 - Q值越大,陷波越窄,对其他频率影响越小 - 常用于去除工频干扰

移动平均滤波器

移动平均滤波器是最简单的低通滤波器,计算效率高。

// 移动平均滤波器
typedef struct {
    int16_t* buffer;
    uint16_t length;
    uint16_t index;
    int32_t sum;
} MovingAverage_t;

// 初始化移动平均滤波器
void MovingAverage_Init(MovingAverage_t* ma, uint16_t length) {
    ma->length = length;
    ma->index = 0;
    ma->sum = 0;

    ma->buffer = (int16_t*)malloc(length * sizeof(int16_t));
    memset(ma->buffer, 0, length * sizeof(int16_t));
}

// 移动平均滤波器处理
int16_t MovingAverage_Process(MovingAverage_t* ma, int16_t input) {
    // 减去最旧的样本
    ma->sum -= ma->buffer[ma->index];

    // 加上新样本
    ma->buffer[ma->index] = input;
    ma->sum += input;

    // 更新索引
    ma->index++;
    if (ma->index >= ma->length) {
        ma->index = 0;
    }

    // 返回平均值
    return (int16_t)(ma->sum / ma->length);
}

// 优化版本:长度为2的幂次,使用位移代替除法
int16_t MovingAverage_Process_Fast(MovingAverage_t* ma, int16_t input) {
    // 假设length = 16 (2^4)
    ma->sum -= ma->buffer[ma->index];
    ma->buffer[ma->index] = input;
    ma->sum += input;

    ma->index = (ma->index + 1) & 0x0F;  // 等价于 % 16

    return (int16_t)(ma->sum >> 4);  // 等价于 / 16
}

代码说明: - 移动平均滤波器计算简单,适合实时处理 - 使用累加和,避免重复计算 - 长度为2的幂次时,可以使用位运算优化

医疗器械滤波应用

心电信号滤波

// 心电信号滤波器组
typedef struct {
    Biquad_Filter_t notch_50hz;    // 50Hz陷波(工频干扰)
    Biquad_Filter_t hpf_05hz;      // 0.5Hz高通(基线漂移)
    Biquad_Filter_t lpf_40hz;      // 40Hz低通(高频噪声)
} ECG_FilterBank_t;

ECG_FilterBank_t ecg_filters;

// 初始化心电滤波器组
void ECG_FilterBank_Init(ECG_FilterBank_t* fb) {
    // 50Hz陷波滤波器
    Init_Notch_50Hz(&fb->notch_50hz);

    // 0.5Hz高通滤波器(去除基线漂移)
    float b0_hp = 0.9899f;
    float b1_hp = -1.9799f;
    float b2_hp = 0.9899f;
    float a1_hp = -1.9798f;
    float a2_hp = 0.9798f;
    Biquad_Init(&fb->hpf_05hz, b0_hp, b1_hp, b2_hp, a1_hp, a2_hp);

    // 40Hz低通滤波器(去除高频噪声)
    float b0_lp = 0.0675f;
    float b1_lp = 0.1349f;
    float b2_lp = 0.0675f;
    float a1_lp = -1.1430f;
    float a2_lp = 0.4128f;
    Biquad_Init(&fb->lpf_40hz, b0_lp, b1_lp, b2_lp, a1_lp, a2_lp);
}

// 心电信号滤波处理
float ECG_Filter_Process(ECG_FilterBank_t* fb, float ecg_raw) {
    float signal = ecg_raw;

    // 级联滤波
    signal = Biquad_Process(&fb->hpf_05hz, signal);   // 去基线漂移
    signal = Biquad_Process(&fb->notch_50hz, signal); // 去工频干扰
    signal = Biquad_Process(&fb->lpf_40hz, signal);   // 去高频噪声

    return signal;
}

// 使用示例
void ECG_Processing_Task(void) {
    // 从ADC读取原始心电信号
    float ecg_raw = ADC_ReadECG();

    // 滤波处理
    float ecg_filtered = ECG_Filter_Process(&ecg_filters, ecg_raw);

    // 后续处理(QRS检测等)
    ECG_QRS_Detection(ecg_filtered);
}

代码说明: - 级联多个滤波器,分别处理不同类型的干扰 - 高通滤波器去除基线漂移 - 陷波滤波器去除工频干扰 - 低通滤波器去除高频噪声

血氧信号滤波

// 血氧信号滤波器
typedef struct {
    Biquad_Filter_t lpf_5hz;       // 5Hz低通(平滑脉搏波)
    MovingAverage_t ma_dc;         // DC分量提取
} SpO2_Filter_t;

SpO2_Filter_t spo2_filter;

// 初始化血氧滤波器
void SpO2_Filter_Init(SpO2_Filter_t* sf) {
    // 5Hz低通滤波器
    float b0 = 0.0013f;
    float b1 = 0.0025f;
    float b2 = 0.0013f;
    float a1 = -1.8226f;
    float a2 = 0.8278f;
    Biquad_Init(&sf->lpf_5hz, b0, b1, b2, a1, a2);

    // 移动平均滤波器(提取DC分量)
    MovingAverage_Init(&sf->ma_dc, 64);
}

// 血氧信号处理
void SpO2_Process(SpO2_Filter_t* sf, float red_raw, float ir_raw) {
    // 低通滤波
    float red_ac = Biquad_Process(&sf->lpf_5hz, red_raw);
    float ir_ac = Biquad_Process(&sf->lpf_5hz, ir_raw);

    // 提取DC分量
    float red_dc = MovingAverage_Process(&sf->ma_dc, (int16_t)red_raw);
    float ir_dc = MovingAverage_Process(&sf->ma_dc, (int16_t)ir_raw);

    // 计算AC/DC比值
    float ratio = (red_ac / red_dc) / (ir_ac / ir_dc);

    // 根据比值计算SpO2
    float spo2 = 110.0f - 25.0f * ratio;

    // 限制范围
    if (spo2 > 100.0f) spo2 = 100.0f;
    if (spo2 < 0.0f) spo2 = 0.0f;
}

代码说明: - 低通滤波器提取AC分量(脉搏波动) - 移动平均滤波器提取DC分量(平均光强) - AC/DC比值用于计算血氧饱和度

滤波器设计工具

使用Python设计滤波器

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

# 设计FIR低通滤波器
def design_fir_lpf(cutoff_hz, fs_hz, numtaps=31):
    """
    设计FIR低通滤波器
    cutoff_hz: 截止频率(Hz)
    fs_hz: 采样率(Hz)
    numtaps: 滤波器阶数+1
    """
    # 归一化截止频率
    nyquist = fs_hz / 2.0
    cutoff_norm = cutoff_hz / nyquist

    # 使用窗函数法设计
    coeffs = signal.firwin(numtaps, cutoff_norm, window='hamming')

    return coeffs

# 设计IIR巴特沃斯低通滤波器
def design_butterworth_lpf(cutoff_hz, fs_hz, order=2):
    """
    设计IIR巴特沃斯低通滤波器
    cutoff_hz: 截止频率(Hz)
    fs_hz: 采样率(Hz)
    order: 滤波器阶数
    """
    # 归一化截止频率
    nyquist = fs_hz / 2.0
    cutoff_norm = cutoff_hz / nyquist

    # 设计滤波器
    b, a = signal.butter(order, cutoff_norm, btype='low')

    # 转换为二阶节(SOS)形式
    sos = signal.butter(order, cutoff_norm, btype='low', output='sos')

    return b, a, sos

# 设计陷波滤波器
def design_notch_filter(f0_hz, fs_hz, Q=30):
    """
    设计陷波滤波器
    f0_hz: 陷波频率(Hz)
    fs_hz: 采样率(Hz)
    Q: 品质因数
    """
    # 归一化频率
    w0 = f0_hz / (fs_hz / 2.0)

    # 设计陷波滤波器
    b, a = signal.iirnotch(w0, Q)

    return b, a

# 生成C代码
def generate_c_code(coeffs, name):
    """生成C语言数组定义"""
    print(f"const float {name}[{len(coeffs)}] = {{")
    for i in range(0, len(coeffs), 5):
        line = "    "
        for j in range(5):
            if i + j < len(coeffs):
                line += f"{coeffs[i+j]:8.4f}f, "
        print(line)
    print("};")

# 示例:设计50Hz低通滤波器
fs = 1000  # 采样率1000Hz
cutoff = 50  # 截止频率50Hz

# FIR滤波器
fir_coeffs = design_fir_lpf(cutoff, fs, numtaps=31)
print("FIR Low-Pass Filter Coefficients:")
generate_c_code(fir_coeffs, "lpf_50hz_coeffs")

# IIR滤波器
b, a, sos = design_butterworth_lpf(cutoff, fs, order=2)
print("\nIIR Butterworth Filter Coefficients:")
print(f"b0={b[0]:.4f}, b1={b[1]:.4f}, b2={b[2]:.4f}")
print(f"a0={a[0]:.4f}, a1={a[1]:.4f}, a2={a[2]:.4f}")

代码说明: - 使用scipy.signal库设计滤波器 - 自动生成C语言代码 - 支持FIR和IIR滤波器设计

性能优化

使用CMSIS-DSP库

#include "arm_math.h"

// 使用CMSIS-DSP的FIR滤波器
arm_fir_instance_f32 fir_instance;
float fir_state[31 + 1 - 1];  // 状态缓冲区
float fir_coeffs[31];          // 滤波器系数

// 初始化
void Init_CMSIS_FIR(void) {
    arm_fir_init_f32(&fir_instance, 31, fir_coeffs, fir_state, 1);
}

// 处理
void Process_CMSIS_FIR(float* input, float* output, uint32_t blockSize) {
    arm_fir_f32(&fir_instance, input, output, blockSize);
}

// 使用CMSIS-DSP的Biquad滤波器
arm_biquad_casd_df1_inst_f32 biquad_instance;
float biquad_state[4];  // 2阶滤波器需要4个状态
float biquad_coeffs[5]; // b0, b1, b2, a1, a2

// 初始化
void Init_CMSIS_Biquad(void) {
    arm_biquad_cascade_df1_init_f32(&biquad_instance, 1, 
                                     biquad_coeffs, biquad_state);
}

// 处理
void Process_CMSIS_Biquad(float* input, float* output, uint32_t blockSize) {
    arm_biquad_cascade_df1_f32(&biquad_instance, input, output, blockSize);
}

代码说明: - CMSIS-DSP库提供优化的DSP函数 - 利用ARM Cortex-M的DSP指令 - 性能比普通C代码快2-10倍

最佳实践

数字滤波器设计最佳实践

  • 选择合适的滤波器类型:FIR线性相位,IIR计算效率高
  • 合理设计参数:截止频率、阶数、采样率等
  • 使用设计工具:MATLAB、Python等工具辅助设计
  • 定点实现:嵌入式系统优先使用定点运算
  • 使用优化库:CMSIS-DSP等优化库提高性能
  • 级联设计:复杂滤波器可以级联多个简单滤波器
  • 测试验证:使用实际信号测试滤波效果
  • 考虑延迟:FIR滤波器有固定延迟,需要考虑实时性
  • 稳定性检查:IIR滤波器需要检查极点位置,确保稳定
  • 量化误差:定点实现需要考虑量化误差和溢出

常见陷阱

注意事项

  • 系数精度不足:定点实现时系数量化误差大
  • 溢出问题:累加器溢出导致结果错误
  • 不稳定的IIR:极点在单位圆外,滤波器不稳定
  • 过高的阶数:计算量大,实时性差
  • 采样率不足:违反奈奎斯特定理,产生混叠
  • 初始化错误:历史状态未清零,产生瞬态响应
  • 相位失真:IIR滤波器非线性相位,可能影响波形
  • 边界效应:滤波器启动时的瞬态响应
  • 数值精度:浮点运算精度问题
  • 内存泄漏:动态分配的缓冲区未释放

实践练习

  1. FIR滤波器设计:设计并实现一个50Hz低通FIR滤波器
  2. IIR滤波器设计:设计并实现一个50Hz陷波IIR滤波器
  3. 心电信号滤波:实现完整的心电信号滤波器组
  4. 性能优化:使用CMSIS-DSP库优化滤波器性能
  5. 定点实现:将浮点滤波器转换为定点实现

自测问题

FIR和IIR滤波器有什么区别?如何选择?

FIR和IIR滤波器各有优缺点,需要根据应用场景选择。

答案

主要区别

特性 FIR滤波器 IIR滤波器
脉冲响应 有限 无限
相位特性 线性相位 非线性相位
稳定性 总是稳定 可能不稳定
计算量 较大 较小
延迟 固定延迟 延迟较小
设计难度 简单 较复杂

选择建议

  1. 选择FIR的场景
  2. 需要线性相位(如音频处理)
  3. 稳定性要求高
  4. 计算资源充足
  5. 可以接受较大延迟

  6. 选择IIR的场景

  7. 计算资源有限
  8. 需要陡峭的过渡带
  9. 实时性要求高
  10. 不关心相位失真

  11. 医疗器械应用

  12. 心电信号:IIR(实时性要求高)
  13. 音频信号:FIR(需要线性相位)
  14. 血氧信号:IIR(计算简单)
  15. 脑电信号:根据具体需求选择
如何设计一个50Hz陷波滤波器去除工频干扰?

陷波滤波器用于抑制特定频率的干扰。

答案

设计步骤

  1. 确定参数
  2. 陷波频率:f0 = 50Hz(或60Hz)
  3. 采样率:fs = 1000Hz
  4. 品质因数:Q = 30(Q越大,陷波越窄)

  5. 计算系数

    w0 = 2π × f0 / fs
    α = sin(w0) / (2Q)
    
    b0 = 1
    b1 = -2 × cos(w0)
    b2 = 1
    a0 = 1 + α
    a1 = -2 × cos(w0)
    a2 = 1 - α
    
    归一化:除以a0
    

  6. 实现代码

    void Init_Notch_50Hz(Biquad_Filter_t* bq) {
        float fs = 1000.0f;
        float f0 = 50.0f;
        float Q = 30.0f;
    
        float w0 = 2.0f * M_PI * f0 / fs;
        float alpha = sinf(w0) / (2.0f * Q);
    
        float b0 = 1.0f;
        float b1 = -2.0f * cosf(w0);
        float b2 = 1.0f;
        float a0 = 1.0f + alpha;
        float a1 = -2.0f * cosf(w0);
        float a2 = 1.0f - alpha;
    
        b0 /= a0; b1 /= a0; b2 /= a0;
        a1 /= a0; a2 /= a0;
    
        Biquad_Init(bq, b0, b1, b2, a1, a2);
    }
    

  7. 验证效果

  8. 在50Hz处衰减应>40dB
  9. 其他频率影响应<1dB
  10. 使用频谱分析验证

注意事项: - Q值不宜过大(>50),可能导致不稳定 - 需要精确的采样率 - 可以级联多个陷波滤波器(50Hz和100Hz)

在嵌入式系统中如何优化滤波器性能?

嵌入式系统资源有限,需要优化滤波器实现。

答案

优化方法

  1. 使用定点运算
  2. 浮点运算慢,定点运算快
  3. 使用Q15或Q31格式
  4. 注意溢出和精度
// Q15定点乘法
int16_t q15_mult(int16_t a, int16_t b) {
    return (int16_t)(((int32_t)a * b) >> 15);
}
  1. 使用CMSIS-DSP库
  2. ARM官方优化的DSP库
  3. 利用SIMD指令
  4. 性能提升2-10倍
#include "arm_math.h"
arm_fir_f32(&fir_instance, input, output, blockSize);
  1. 降低滤波器阶数
  2. 使用IIR代替高阶FIR
  3. 级联多个低阶滤波器
  4. 权衡性能和精度

  5. 批量处理

  6. 一次处理多个样本
  7. 减少函数调用开销
  8. 提高缓存命中率

  9. 使用DMA

  10. DMA传输ADC数据
  11. CPU可以处理其他任务
  12. 降低CPU负载

  13. 优化数据结构

  14. 使用循环缓冲区
  15. 避免数据搬移
  16. 内存对齐

  17. 编译器优化

  18. 使用-O2或-O3优化级别
  19. 使用inline函数
  20. 循环展开

性能对比: - 浮点FIR(31阶):~1000 cycles - 定点FIR(31阶):~500 cycles - CMSIS FIR(31阶):~200 cycles - IIR(2阶):~50 cycles

如何处理滤波器的初始瞬态响应?

滤波器启动时会有瞬态响应,可能影响结果。

答案

瞬态响应问题

  1. 产生原因
  2. 滤波器历史状态为零
  3. 输入信号突变
  4. 初始条件不匹配

  5. 影响

  6. 输出信号初期失真
  7. 可能触发误报警
  8. 影响测量精度

处理方法

  1. 预填充历史

    void FIR_Prefill(FIR_Filter_t* fir, float initial_value) {
        for (int i = 0; i < fir->length; i++) {
            fir->buffer[i] = initial_value;
        }
    }
    

  2. 丢弃初始样本

    // 丢弃前N个样本
    for (int i = 0; i < fir->length; i++) {
        FIR_Process(fir, input[i]);  // 不使用输出
    }
    
    // 之后的样本才有效
    for (int i = fir->length; i < total_samples; i++) {
        output[i] = FIR_Process(fir, input[i]);
    }
    

  3. 使用稳态初始化

    // 使用输入信号的平均值初始化
    float avg = CalculateAverage(input, 100);
    FIR_Prefill(fir, avg);
    

  4. 渐进启动

    // 逐渐增加滤波器权重
    float weight = 0.0f;
    for (int i = 0; i < fir->length; i++) {
        weight = (float)i / fir->length;
        output[i] = weight * FIR_Process(fir, input[i]) + 
                   (1.0f - weight) * input[i];
    }
    

医疗器械应用: - 心电监护:丢弃前2秒数据 - 血氧测量:预填充平均值 - 血压测量:使用渐进启动

如何验证滤波器的正确性和性能?

滤波器设计完成后需要充分验证。

答案

验证方法

  1. 频率响应测试
  2. 输入不同频率的正弦波
  3. 测量输出幅度和相位
  4. 绘制频率响应曲线
// 测试频率响应
for (float freq = 1.0f; freq <= 500.0f; freq += 1.0f) {
    float amplitude = TestFrequencyResponse(filter, freq, fs);
    printf("%.1f Hz: %.2f dB\n", freq, 20*log10(amplitude));
}
  1. 脉冲响应测试
  2. 输入单位脉冲
  3. 观察输出波形
  4. 验证FIR长度和IIR衰减

  5. 阶跃响应测试

  6. 输入阶跃信号
  7. 测量上升时间和超调
  8. 验证瞬态特性

  9. 实际信号测试

  10. 使用真实的医疗信号
  11. 对比滤波前后的效果
  12. 验证噪声抑制能力

  13. 性能测试

  14. 测量执行时间
  15. 测量内存占用
  16. 验证实时性
uint32_t start = DWT->CYCCNT;
float output = FIR_Process(&fir, input);
uint32_t cycles = DWT->CYCCNT - start;
  1. 稳定性测试
  2. 长时间运行
  3. 极端输入测试
  4. 验证无溢出和发散

  5. 精度测试

  6. 对比浮点和定点实现
  7. 测量量化误差
  8. 验证满足精度要求

验证工具: - MATLAB/Python:离线分析 - 示波器:实时波形观察 - 频谱分析仪:频域分析 - 逻辑分析仪:时序分析

相关资源

参考文献

  1. "Digital Signal Processing" - John G. Proakis, Dimitris G. Manolakis
  2. "The Scientist and Engineer's Guide to Digital Signal Processing" - Steven W. Smith
  3. CMSIS-DSP Software Library - ARM
  4. IEC 60601-2-27:2011 - Particular requirements for electrocardiographic monitoring equipment
  5. "Embedded Signal Processing with the Micro Signal Architecture" - Woon-Seng Gan, Sen M. Kuo

💬 讨论区

欢迎在这里分享您的想法、提出问题或参与讨论。需要 GitHub 账号登录。