数字滤波器 C 语言实现大全
·
数字滤波器 C 语言实现库,包含 FIR滤波器、IIR滤波器、自适应滤波器、多速率滤波器 等常用算法,适用于嵌入式系统、信号处理、音频处理等场景。
一、项目结构
DigitalFilters/
├── include/
│ ├── fir_filter.h ★ FIR滤波器
│ ├── iir_filter.h ★ IIR滤波器
│ ├── adaptive_filter.h ★ 自适应滤波器
│ ├── multirate_filter.h ★ 多速率滤波器
│ ├── filter_design.h ★ 滤波器设计
│ ├── window_functions.h ★ 窗函数
│ └── math_utils.h ★ 数学工具
├── src/
│ ├── fir_filter.c
│ ├── iir_filter.c
│ ├── adaptive_filter.c
│ ├── multirate_filter.c
│ ├── filter_design.c
│ ├── window_functions.c
│ └── math_utils.c
├── examples/
│ ├── fir_lowpass_example.c
│ ├── iir_bandpass_example.c
│ ├── lms_adaptive_example.c
│ └── decimation_example.c
├── tests/
│ ├── test_fir.c
│ ├── test_iir.c
│ └── test_performance.c
└── Makefile
二、核心代码实现
1. 数学工具头文件 (include/math_utils.h)
#ifndef MATH_UTILS_H
#define MATH_UTILS_H
#include <stdint.h>
#include <stddef.h>
#include <math.h>
// 浮点类型定义
typedef float float32_t;
typedef double float64_t;
// 复数结构
typedef struct {
float32_t real;
float32_t imag;
} complex_float_t;
// 数学常数
#define PI_F 3.14159265358979323846f
#define TWO_PI_F 6.28318530717958647692f
#define SQRT2_F 1.41421356237309504880f
// 数学函数宏
#define ABS(x) ((x) < 0 ? -(x) : (x))
#define MIN(a, b) ((a) < (b) ? (a) : (b))
#define MAX(a, b) ((a) > (b) ? (a) : (b))
#define CLAMP(x, min, max) (MIN(MAX((x), (min)), (max)))
// 正弦余弦近似(嵌入式优化)
float32_t fast_sin(float32_t x);
float32_t fast_cos(float32_t x);
float32_t fast_sqrt(float32_t x);
// 复数运算
complex_float_t complex_add(complex_float_t a, complex_float_t b);
complex_float_t complex_mul(complex_float_t a, complex_float_t b);
float32_t complex_mag(complex_float_t a);
#endif // MATH_UTILS_H
2. FIR 滤波器头文件 (include/fir_filter.h)
#ifndef FIR_FILTER_H
#define FIR_FILTER_H
#include <stdint.h>
#include <stddef.h>
// FIR滤波器结构
typedef struct {
float32_t *coeffs; // 滤波器系数
float32_t *buffer; // 延迟线缓冲区
uint16_t num_taps; // 抽头数
uint16_t buffer_index; // 缓冲区索引
uint8_t initialized; // 初始化标志
} fir_filter_t;
// FIR滤波器配置
typedef struct {
float32_t *coeffs; // 系数数组
uint16_t num_taps; // 抽头数
} fir_config_t;
// 滤波器类型
typedef enum {
FIR_LOWPASS = 0,
FIR_HIGHPASS,
FIR_BANDPASS,
FIR_BANDSTOP,
FIR_ALLPASS
} fir_filter_type_t;
// 窗函数类型
typedef enum {
WINDOW_RECTANGULAR = 0,
WINDOW_HANN,
WINDOW_HAMMING,
WINDOW_BLACKMAN,
WINDOW_KAISER
} window_type_t;
// 函数声明
int fir_init(fir_filter_t *filter, const fir_config_t *config);
float32_t fir_process(fir_filter_t *filter, float32_t input);
void fir_reset(fir_filter_t *filter);
void fir_set_coeffs(fir_filter_t *filter, const float32_t *coeffs);
void fir_get_coeffs(const fir_filter_t *filter, float32_t *coeffs);
void fir_deinit(fir_filter_t *filter);
// 批量处理函数
void fir_process_block(fir_filter_t *filter,
const float32_t *input,
float32_t *output,
uint16_t block_size);
// 滤波器设计函数
int fir_design_lowpass(fir_filter_t *filter,
float32_t cutoff_freq,
float32_t sample_rate,
uint16_t num_taps,
window_type_t window);
int fir_design_highpass(fir_filter_t *filter,
float32_t cutoff_freq,
float32_t sample_rate,
uint16_t num_taps,
window_type_t window);
int fir_design_bandpass(fir_filter_t *filter,
float32_t low_cutoff,
float32_t high_cutoff,
float32_t sample_rate,
uint16_t num_taps,
window_type_t window);
// 对称FIR优化版本(偶数抽头)
float32_t fir_process_symmetrical(fir_filter_t *filter, float32_t input);
// 实时处理宏(高性能)
#define FIR_PROCESS_MACRO(filter, input) \
({ \
float32_t output = 0.0f; \
filter.buffer[filter.buffer_index] = input; \
uint16_t idx = filter.buffer_index; \
for (uint16_t i = 0; i < filter.num_taps; i++) { \
output += filter.coeffs[i] * filter.buffer[idx]; \
if (idx == 0) idx = filter.num_taps - 1; \
else idx--; \
} \
filter.buffer_index = (filter.buffer_index + 1) % filter.num_taps; \
output; \
})
#endif // FIR_FILTER_H
3. FIR 滤波器实现 (src/fir_filter.c)
#include "fir_filter.h"
#include "math_utils.h"
#include "window_functions.h"
#include <stdlib.h>
#include <string.h>
// 初始化FIR滤波器
int fir_init(fir_filter_t *filter, const fir_config_t *config) {
if (!filter || !config || config->num_taps == 0) {
return -1;
}
// 分配内存
filter->coeffs = (float32_t *)malloc(config->num_taps * sizeof(float32_t));
filter->buffer = (float32_t *)malloc(config->num_taps * sizeof(float32_t));
if (!filter->coeffs || !filter->buffer) {
fir_deinit(filter);
return -1;
}
// 复制系数
memcpy(filter->coeffs, config->coeffs, config->num_taps * sizeof(float32_t));
// 清零缓冲区
memset(filter->buffer, 0, config->num_taps * sizeof(float32_t));
filter->num_taps = config->num_taps;
filter->buffer_index = 0;
filter->initialized = 1;
return 0;
}
// 处理单个样本
float32_t fir_process(fir_filter_t *filter, float32_t input) {
if (!filter || !filter->initialized) {
return 0.0f;
}
float32_t output = 0.0f;
uint16_t i;
// 将新样本存入缓冲区
filter->buffer[filter->buffer_index] = input;
// 卷积计算
uint16_t buffer_idx = filter->buffer_index;
for (i = 0; i < filter->num_taps; i++) {
output += filter->coeffs[i] * filter->buffer[buffer_idx];
// 循环缓冲区索引
if (buffer_idx == 0) {
buffer_idx = filter->num_taps - 1;
} else {
buffer_idx--;
}
}
// 更新缓冲区索引
filter->buffer_index = (filter->buffer_index + 1) % filter->num_taps;
return output;
}
// 处理数据块
void fir_process_block(fir_filter_t *filter,
const float32_t *input,
float32_t *output,
uint16_t block_size) {
if (!filter || !input || !output || !filter->initialized) {
return;
}
for (uint16_t i = 0; i < block_size; i++) {
output[i] = fir_process(filter, input[i]);
}
}
// 对称FIR优化处理(适用于线性相位滤波器)
float32_t fir_process_symmetrical(fir_filter_t *filter, float32_t input) {
if (!filter || !filter->initialized) {
return 0.0f;
}
float32_t output = 0.0f;
uint16_t half_taps = filter->num_taps / 2;
// 存储新样本
filter->buffer[filter->buffer_index] = input;
// 对称卷积
uint16_t idx1 = filter->buffer_index;
uint16_t idx2 = (filter->buffer_index + half_taps) % filter->num_taps;
for (uint16_t i = 0; i < half_taps; i++) {
output += filter->coeffs[i] * (filter->buffer[idx1] + filter->buffer[idx2]);
// 更新索引
idx1 = (idx1 + 1) % filter->num_taps;
if (idx2 == 0) {
idx2 = filter->num_taps - 1;
} else {
idx2--;
}
}
// 处理中心样本(奇数抽头)
if (filter->num_taps % 2 == 1) {
output += filter->coeffs[half_taps] * filter->buffer[idx1];
}
filter->buffer_index = (filter->buffer_index + 1) % filter->num_taps;
return output;
}
// 重置滤波器状态
void fir_reset(fir_filter_t *filter) {
if (!filter || !filter->initialized) {
return;
}
memset(filter->buffer, 0, filter->num_taps * sizeof(float32_t));
filter->buffer_index = 0;
}
// 设置滤波器系数
void fir_set_coeffs(fir_filter_t *filter, const float32_t *coeffs) {
if (!filter || !coeffs || !filter->initialized) {
return;
}
memcpy(filter->coeffs, coeffs, filter->num_taps * sizeof(float32_t));
}
// 获取滤波器系数
void fir_get_coeffs(const fir_filter_t *filter, float32_t *coeffs) {
if (!filter || !coeffs || !filter->initialized) {
return;
}
memcpy(coeffs, filter->coeffs, filter->num_taps * sizeof(float32_t));
}
// 设计低通FIR滤波器
int fir_design_lowpass(fir_filter_t *filter,
float32_t cutoff_freq,
float32_t sample_rate,
uint16_t num_taps,
window_type_t window) {
if (!filter || num_taps == 0 || cutoff_freq <= 0 || sample_rate <= 0) {
return -1;
}
// 归一化截止频率
float32_t norm_cutoff = cutoff_freq / sample_rate;
// 生成理想低通滤波器系数
float32_t *coeffs = (float32_t *)malloc(num_taps * sizeof(float32_t));
if (!coeffs) {
return -1;
}
int center = num_taps / 2;
float32_t sum = 0.0f;
for (int i = 0; i < num_taps; i++) {
if (i == center) {
coeffs[i] = 2.0f * norm_cutoff;
} else {
float32_t angle = PI_F * (i - center) * 2.0f * norm_cutoff;
coeffs[i] = sinf(angle) / (PI_F * (i - center));
}
sum += coeffs[i];
}
// 归一化
for (int i = 0; i < num_taps; i++) {
coeffs[i] /= sum;
}
// 应用窗函数
apply_window(coeffs, num_taps, window);
// 初始化滤波器
fir_config_t config = {
.coeffs = coeffs,
.num_taps = num_taps
};
int result = fir_init(filter, &config);
free(coeffs);
return result;
}
// 设计高通FIR滤波器
int fir_design_highpass(fir_filter_t *filter,
float32_t cutoff_freq,
float32_t sample_rate,
uint16_t num_taps,
window_type_t window) {
if (!filter || num_taps == 0 || cutoff_freq <= 0 || sample_rate <= 0) {
return -1;
}
// 先设计低通滤波器
if (fir_design_lowpass(filter, cutoff_freq, sample_rate, num_taps, window) != 0) {
return -1;
}
// 转换为高通:系数取反(除了中心系数)
int center = num_taps / 2;
for (int i = 0; i < num_taps; i++) {
if (i != center) {
filter->coeffs[i] = -filter->coeffs[i];
}
}
return 0;
}
// 设计带通FIR滤波器
int fir_design_bandpass(fir_filter_t *filter,
float32_t low_cutoff,
float32_t high_cutoff,
float32_t sample_rate,
uint16_t num_taps,
window_type_t window) {
if (!filter || num_taps == 0 || low_cutoff <= 0 || high_cutoff <= 0 ||
low_cutoff >= high_cutoff || sample_rate <= 0) {
return -1;
}
// 设计两个低通滤波器
fir_filter_t low_filter1, low_filter2;
if (fir_design_lowpass(&low_filter1, high_cutoff, sample_rate, num_taps, window) != 0 ||
fir_design_lowpass(&low_filter2, low_cutoff, sample_rate, num_taps, window) != 0) {
return -1;
}
// 带通 = 高通(low_cutoff) - 低通(high_cutoff)
for (int i = 0; i < num_taps; i++) {
filter->coeffs[i] = low_filter2.coeffs[i] - low_filter1.coeffs[i];
}
// 释放临时滤波器
fir_deinit(&low_filter1);
fir_deinit(&low_filter2);
return 0;
}
// 释放滤波器资源
void fir_deinit(fir_filter_t *filter) {
if (!filter) {
return;
}
if (filter->coeffs) {
free(filter->coeffs);
filter->coeffs = NULL;
}
if (filter->buffer) {
free(filter->buffer);
filter->buffer = NULL;
}
filter->initialized = 0;
filter->num_taps = 0;
filter->buffer_index = 0;
}
4. IIR 滤波器头文件 (include/iir_filter.h)
#ifndef IIR_FILTER_H
#define IIR_FILTER_H
#include <stdint.h>
#include <stddef.h>
// IIR滤波器结构
typedef struct {
float32_t *b_coeffs; // 分子系数 (feedforward)
float32_t *a_coeffs; // 分母系数 (feedback)
float32_t *x_buffer; // 输入延迟线
float32_t *y_buffer; // 输出延迟线
uint8_t num_b; // 分子阶数
uint8_t num_a; // 分母阶数
uint8_t x_index; // 输入缓冲区索引
uint8_t y_index; // 输出缓冲区索引
uint8_t initialized; // 初始化标志
} iir_filter_t;
// IIR滤波器配置
typedef struct {
float32_t *b_coeffs; // 分子系数
float32_t *a_coeffs; // 分母系数
uint8_t num_b; // 分子阶数
uint8_t num_a; // 分母阶数
} iir_config_t;
// IIR滤波器类型
typedef enum {
IIR_LOWPASS = 0,
IIR_HIGHPASS,
IIR_BANDPASS,
IIR_BANDSTOP,
IIR_NOTCH
} iir_filter_type_t;
// 双二阶滤波器结构(二阶节)
typedef struct {
float32_t b0, b1, b2; // 分子系数
float32_t a1, a2; // 分母系数(a0=1)
float32_t x1, x2; // 输入延迟
float32_t y1, y2; // 输出延迟
} biquad_t;
// 函数声明
int iir_init(iir_filter_t *filter, const iir_config_t *config);
float32_t iir_process(iir_filter_t *filter, float32_t input);
void iir_reset(iir_filter_t *filter);
void iir_deinit(iir_filter_t *filter);
// 批量处理
void iir_process_block(iir_filter_t *filter,
const float32_t *input,
float32_t *output,
uint16_t block_size);
// 双二阶滤波器(更高效)
int biquad_init(biquad_t *filter,
float32_t freq,
float32_t q,
float32_t gain,
iir_filter_type_t type,
float32_t sample_rate);
float32_t biquad_process(biquad_t *filter, float32_t input);
void biquad_reset(biquad_t *filter);
// 级联双二阶滤波器
typedef struct {
biquad_t *stages;
uint8_t num_stages;
} cascade_biquad_t;
int cascade_biquad_init(cascade_biquad_t *cascade,
const float32_t *freqs,
const float32_t *qs,
const float32_t *gains,
const iir_filter_type_t *types,
uint8_t num_stages,
float32_t sample_rate);
float32_t cascade_biquad_process(cascade_biquad_t *cascade, float32_t input);
void cascade_biquad_reset(cascade_biquad_t *cascade);
void cascade_biquad_deinit(cascade_biquad_t *cascade);
// 滤波器设计
int iir_design_lowpass(iir_filter_t *filter,
float32_t cutoff_freq,
float32_t sample_rate,
uint8_t order);
int iir_design_highpass(iir_filter_t *filter,
float32_t cutoff_freq,
float32_t sample_rate,
uint8_t order);
// 直接形式II(规范型)实现
float32_t iir_process_df2(iir_filter_t *filter, float32_t input);
#endif // IIR_FILTER_H
5. IIR 滤波器实现 (src/iir_filter.c)
#include "iir_filter.h"
#include "math_utils.h"
#include <stdlib.h>
#include <string.h>
// 初始化IIR滤波器
int iir_init(iir_filter_t *filter, const iir_config_t *config) {
if (!filter || !config || config->num_b == 0 || config->num_a == 0) {
return -1;
}
// 分配内存
filter->b_coeffs = (float32_t *)malloc(config->num_b * sizeof(float32_t));
filter->a_coeffs = (float32_t *)malloc(config->num_a * sizeof(float32_t));
filter->x_buffer = (float32_t *)malloc(config->num_b * sizeof(float32_t));
filter->y_buffer = (float32_t *)malloc((config->num_a - 1) * sizeof(float32_t));
if (!filter->b_coeffs || !filter->a_coeffs || !filter->x_buffer || !filter->y_buffer) {
iir_deinit(filter);
return -1;
}
// 复制系数
memcpy(filter->b_coeffs, config->b_coeffs, config->num_b * sizeof(float32_t));
memcpy(filter->a_coeffs, config->a_coeffs, config->num_a * sizeof(float32_t));
// 清零缓冲区
memset(filter->x_buffer, 0, config->num_b * sizeof(float32_t));
memset(filter->y_buffer, 0, (config->num_a - 1) * sizeof(float32_t));
filter->num_b = config->num_b;
filter->num_a = config->num_a;
filter->x_index = 0;
filter->y_index = 0;
filter->initialized = 1;
return 0;
}
// 处理单个样本(直接形式I)
float32_t iir_process(iir_filter_t *filter, float32_t input) {
if (!filter || !filter->initialized) {
return 0.0f;
}
float32_t output = 0.0f;
uint8_t i;
// 存储输入样本
filter->x_buffer[filter->x_index] = input;
// 计算输出:y[n] = sum(b_k * x[n-k]) - sum(a_k * y[n-k])
uint8_t x_idx = filter->x_index;
uint8_t y_idx = filter->y_index;
// 分子部分(前馈)
for (i = 0; i < filter->num_b; i++) {
output += filter->b_coeffs[i] * filter->x_buffer[x_idx];
if (x_idx == 0) {
x_idx = filter->num_b - 1;
} else {
x_idx--;
}
}
// 分母部分(反馈)
for (i = 0; i < filter->num_a - 1; i++) {
output -= filter->a_coeffs[i + 1] * filter->y_buffer[y_idx];
if (y_idx == 0) {
y_idx = filter->num_a - 2;
} else {
y_idx--;
}
}
// 除以a0(通常a0=1,这里省略除法以提高效率)
output /= filter->a_coeffs[0];
// 存储输出样本
filter->y_buffer[filter->y_index] = output;
// 更新索引
filter->x_index = (filter->x_index + 1) % filter->num_b;
filter->y_index = (filter->y_index + 1) % (filter->num_a - 1);
return output;
}
// 直接形式II(规范型)实现
float32_t iir_process_df2(iir_filter_t *filter, float32_t input) {
if (!filter || !filter->initialized) {
return 0.0f;
}
float32_t w = input;
float32_t output = 0.0f;
uint8_t i;
// 计算反馈部分
for (i = 1; i < filter->num_a; i++) {
w -= filter->a_coeffs[i] * filter->y_buffer[(filter->y_index + i - 1) % (filter->num_a - 1)];
}
// 计算前馈部分
for (i = 0; i < filter->num_b; i++) {
output += filter->b_coeffs[i] * filter->y_buffer[(filter->y_index + i) % (filter->num_b - 1)];
}
// 更新延迟线
filter->y_buffer[filter->y_index] = w;
filter->y_index = (filter->y_index + 1) % (filter->num_a - 1);
return output;
}
// 双二阶滤波器处理
float32_t biquad_process(biquad_t *filter, float32_t input) {
if (!filter) {
return 0.0f;
}
// 直接形式II转置
float32_t output = filter->b0 * input + filter->x1;
filter->x1 = filter->b1 * input - filter->a1 * output + filter->x2;
filter->x2 = filter->b2 * input - filter->a2 * output;
return output;
}
// 设计双二阶低通滤波器
int biquad_init(biquad_t *filter,
float32_t freq,
float32_t q,
float32_t gain,
iir_filter_type_t type,
float32_t sample_rate) {
if (!filter || freq <= 0 || q <= 0 || sample_rate <= 0) {
return -1;
}
float32_t omega = TWO_PI_F * freq / sample_rate;
float32_t sn = sinf(omega);
float32_t cs = cosf(omega);
float32_t alpha = sn / (2.0f * q);
float32_t A = powf(10.0f, gain / 40.0f);
float32_t beta = sqrtf(A + A);
switch (type) {
case IIR_LOWPASS:
filter->b0 = (1.0f - cs) / 2.0f;
filter->b1 = 1.0f - cs;
filter->b2 = (1.0f - cs) / 2.0f;
filter->a1 = -2.0f * cs;
filter->a2 = 1.0f - alpha;
break;
case IIR_HIGHPASS:
filter->b0 = (1.0f + cs) / 2.0f;
filter->b1 = -(1.0f + cs);
filter->b2 = (1.0f + cs) / 2.0f;
filter->a1 = -2.0f * cs;
filter->a2 = 1.0f - alpha;
break;
case IIR_BANDPASS:
filter->b0 = alpha;
filter->b1 = 0.0f;
filter->b2 = -alpha;
filter->a1 = -2.0f * cs;
filter->a2 = 1.0f - alpha;
break;
case IIR_NOTCH:
filter->b0 = 1.0f;
filter->b1 = -2.0f * cs;
filter->b2 = 1.0f;
filter->a1 = -2.0f * cs;
filter->a2 = 1.0f - alpha;
break;
default:
return -1;
}
// 归一化
float32_t a0 = 1.0f + alpha;
filter->b0 /= a0;
filter->b1 /= a0;
filter->b2 /= a0;
filter->a1 /= a0;
filter->a2 /= a0;
// 清零延迟
filter->x1 = filter->x2 = filter->y1 = filter->y2 = 0.0f;
return 0;
}
// 级联双二阶滤波器处理
float32_t cascade_biquad_process(cascade_biquad_t *cascade, float32_t input) {
if (!cascade || !cascade->stages) {
return 0.0f;
}
float32_t output = input;
for (uint8_t i = 0; i < cascade->num_stages; i++) {
output = biquad_process(&cascade->stages[i], output);
}
return output;
}
// 重置IIR滤波器
void iir_reset(iir_filter_t *filter) {
if (!filter || !filter->initialized) {
return;
}
memset(filter->x_buffer, 0, filter->num_b * sizeof(float32_t));
memset(filter->y_buffer, 0, (filter->num_a - 1) * sizeof(float32_t));
filter->x_index = 0;
filter->y_index = 0;
}
// 释放IIR滤波器资源
void iir_deinit(iir_filter_t *filter) {
if (!filter) {
return;
}
if (filter->b_coeffs) {
free(filter->b_coeffs);
filter->b_coeffs = NULL;
}
if (filter->a_coeffs) {
free(filter->a_coeffs);
filter->a_coeffs = NULL;
}
if (filter->x_buffer) {
free(filter->x_buffer);
filter->x_buffer = NULL;
}
if (filter->y_buffer) {
free(filter->y_buffer);
filter->y_buffer = NULL;
}
filter->initialized = 0;
}
6. 自适应滤波器 (include/adaptive_filter.h)
#ifndef ADAPTIVE_FILTER_H
#define ADAPTIVE_FILTER_H
#include <stdint.h>
#include <stddef.h>
// 自适应滤波器类型
typedef enum {
ADAPTIVE_LMS = 0, // 最小均方
ADAPTIVE_NLMS, // 归一化LMS
ADAPTIVE_RLS, // 递归最小二乘
ADAPTIVE_APA // 仿射投影算法
} adaptive_algorithm_t;
// LMS自适应滤波器结构
typedef struct {
float32_t *weights; // 权系数
float32_t *buffer; // 输入缓冲区
float32_t step_size; // 步长
float32_t leakage; // 泄漏因子
uint16_t order; // 滤波器阶数
uint16_t buffer_index; // 缓冲区索引
uint8_t initialized; // 初始化标志
} lms_filter_t;
// NLMS自适应滤波器结构
typedef struct {
float32_t *weights; // 权系数
float32_t *buffer; // 输入缓冲区
float32_t step_size; // 步长
float32_t eps; // 正则化参数
uint16_t order; // 滤波器阶数
uint16_t buffer_index; // 缓冲区索引
uint8_t initialized; // 初始化标志
} nlms_filter_t;
// RLS自适应滤波器结构
typedef struct {
float32_t *weights; // 权系数
float32_t *p_matrix; // 逆相关矩阵
float32_t *buffer; // 输入缓冲区
float32_t lambda; // 遗忘因子
float32_t delta; // 初始化参数
uint16_t order; // 滤波器阶数
uint16_t buffer_index; // 缓冲区索引
uint8_t initialized; // 初始化标志
} rls_filter_t;
// 函数声明
int lms_init(lms_filter_t *filter, uint16_t order, float32_t step_size);
float32_t lms_process(lms_filter_t *filter, float32_t input, float32_t desired);
void lms_reset(lms_filter_t *filter);
void lms_deinit(lms_filter_t *filter);
int nlms_init(nlms_filter_t *filter, uint16_t order, float32_t step_size, float32_t eps);
float32_t nlms_process(nlms_filter_t *filter, float32_t input, float32_t desired);
void nlms_reset(nlms_filter_t *filter);
void nlms_deinit(nlms_filter_t *filter);
int rls_init(rls_filter_t *filter, uint16_t order, float32_t lambda, float32_t delta);
float32_t rls_process(rls_filter_t *filter, float32_t input, float32_t desired);
void rls_reset(rls_filter_t *filter);
void rls_deinit(rls_filter_t *filter);
// 噪声消除器
typedef struct {
lms_filter_t adaptive_filter;
float32_t *noise_reference;
uint16_t buffer_size;
} noise_canceller_t;
int noise_canceller_init(noise_canceller_t *nc, uint16_t order, float32_t step_size, uint16_t buffer_size);
float32_t noise_canceller_process(noise_canceller_t *nc, float32_t signal_plus_noise, float32_t noise_reference);
void noise_canceller_deinit(noise_canceller_t *nc);
#endif // ADAPTIVE_FILTER_H
7. 自适应滤波器实现 (src/adaptive_filter.c)
#include "adaptive_filter.h"
#include "math_utils.h"
#include <stdlib.h>
#include <string.h>
// 初始化LMS滤波器
int lms_init(lms_filter_t *filter, uint16_t order, float32_t step_size) {
if (!filter || order == 0 || step_size <= 0) {
return -1;
}
filter->weights = (float32_t *)malloc(order * sizeof(float32_t));
filter->buffer = (float32_t *)malloc(order * sizeof(float32_t));
if (!filter->weights || !filter->buffer) {
lms_deinit(filter);
return -1;
}
memset(filter->weights, 0, order * sizeof(float32_t));
memset(filter->buffer, 0, order * sizeof(float32_t));
filter->step_size = step_size;
filter->leakage = 1.0f;
filter->order = order;
filter->buffer_index = 0;
filter->initialized = 1;
return 0;
}
// LMS处理
float32_t lms_process(lms_filter_t *filter, float32_t input, float32_t desired) {
if (!filter || !filter->initialized) {
return 0.0f;
}
// 存储输入样本
filter->buffer[filter->buffer_index] = input;
// 计算输出
float32_t output = 0.0f;
uint16_t idx = filter->buffer_index;
for (uint16_t i = 0; i < filter->order; i++) {
output += filter->weights[i] * filter->buffer[idx];
if (idx == 0) {
idx = filter->order - 1;
} else {
idx--;
}
}
// 计算误差
float32_t error = desired - output;
// 更新权系数(带泄漏)
idx = filter->buffer_index;
for (uint16_t i = 0; i < filter->order; i++) {
filter->weights[i] = filter->leakage * filter->weights[i] +
filter->step_size * error * filter->buffer[idx];
if (idx == 0) {
idx = filter->order - 1;
} else {
idx--;
}
}
// 更新缓冲区索引
filter->buffer_index = (filter->buffer_index + 1) % filter->order;
return output;
}
// 初始化NLMS滤波器
int nlms_init(nlms_filter_t *filter, uint16_t order, float32_t step_size, float32_t eps) {
if (!filter || order == 0 || step_size <= 0) {
return -1;
}
filter->weights = (float32_t *)malloc(order * sizeof(float32_t));
filter->buffer = (float32_t *)malloc(order * sizeof(float32_t));
if (!filter->weights || !filter->buffer) {
nlms_deinit(filter);
return -1;
}
memset(filter->weights, 0, order * sizeof(float32_t));
memset(filter->buffer, 0, order * sizeof(float32_t));
filter->step_size = step_size;
filter->eps = eps;
filter->order = order;
filter->buffer_index = 0;
filter->initialized = 1;
return 0;
}
// NLMS处理
float32_t nlms_process(nlms_filter_t *filter, float32_t input, float32_t desired) {
if (!filter || !filter->initialized) {
return 0.0f;
}
// 存储输入样本
filter->buffer[filter->buffer_index] = input;
// 计算输出
float32_t output = 0.0f;
uint16_t idx = filter->buffer_index;
for (uint16_t i = 0; i < filter->order; i++) {
output += filter->weights[i] * filter->buffer[idx];
if (idx == 0) {
idx = filter->order - 1;
} else {
idx--;
}
}
// 计算误差
float32_t error = desired - output;
// 计算输入功率
float32_t power = filter->eps;
idx = filter->buffer_index;
for (uint16_t i = 0; i < filter->order; i++) {
power += filter->buffer[idx] * filter->buffer[idx];
if (idx == 0) {
idx = filter->order - 1;
} else {
idx--;
}
}
// 更新权系数
float32_t mu_norm = filter->step_size / power;
idx = filter->buffer_index;
for (uint16_t i = 0; i < filter->order; i++) {
filter->weights[i] += mu_norm * error * filter->buffer[idx];
if (idx == 0) {
idx = filter->order - 1;
} else {
idx--;
}
}
// 更新缓冲区索引
filter->buffer_index = (filter->buffer_index + 1) % filter->order;
return output;
}
// 噪声消除器处理
float32_t noise_canceller_process(noise_canceller_t *nc, float32_t signal_plus_noise, float32_t noise_reference) {
if (!nc || !nc->adaptive_filter.initialized) {
return 0.0f;
}
// 使用自适应滤波器估计噪声
float32_t estimated_noise = lms_process(&nc->adaptive_filter, noise_reference, signal_plus_noise);
// 从混合信号中减去估计的噪声
float32_t clean_signal = signal_plus_noise - estimated_noise;
return clean_signal;
}
// 释放LMS滤波器资源
void lms_deinit(lms_filter_t *filter) {
if (!filter) {
return;
}
if (filter->weights) {
free(filter->weights);
filter->weights = NULL;
}
if (filter->buffer) {
free(filter->buffer);
filter->buffer = NULL;
}
filter->initialized = 0;
}
8. 多速率滤波器 (include/multirate_filter.h)
#ifndef MULTIRATE_FILTER_H
#define MULTIRATE_FILTER_H
#include <stdint.h>
#include <stddef.h>
#include "fir_filter.h"
// 抽取器结构
typedef struct {
fir_filter_t antialiasing_filter; // 抗混叠滤波器
uint16_t decimation_factor; // 抽取因子
uint16_t phase; // 相位
uint8_t initialized; // 初始化标志
} decimator_t;
// 插值器结构
typedef struct {
fir_filter_t interpolation_filter; // 插值滤波器
uint16_t interpolation_factor; // 插值因子
uint16_t phase; // 相位
uint8_t initialized; // 初始化标志
} interpolator_t;
// 多相滤波器结构
typedef struct {
fir_filter_t **polyphase_filters; // 多相滤波器组
uint16_t num_phases; // 相位数
uint16_t current_phase; // 当前相位
uint8_t initialized; // 初始化标志
} polyphase_filter_t;
// 分数倍采样率转换器
typedef struct {
interpolator_t interpolator; // 插值器
decimator_t decimator; // 抽取器
float32_t rate_ratio; // 采样率比
uint8_t initialized; // 初始化标志
} resampler_t;
// 函数声明
int decimator_init(decimator_t *decimator,
uint16_t decimation_factor,
uint16_t filter_order,
float32_t cutoff_freq,
float32_t sample_rate);
float32_t decimator_process(decimator_t *decimator, float32_t input);
void decimator_process_block(decimator_t *decimator,
const float32_t *input,
float32_t *output,
uint16_t input_size,
uint16_t *output_size);
void decimator_deinit(decimator_t *decimator);
int interpolator_init(interpolator_t *interpolator,
uint16_t interpolation_factor,
uint16_t filter_order,
float32_t cutoff_freq,
float32_t sample_rate);
float32_t interpolator_process(interpolator_t *interpolator, float32_t input);
void interpolator_process_block(interpolator_t *interpolator,
const float32_t *input,
float32_t *output,
uint16_t input_size,
uint16_t *output_size);
void interpolator_deinit(interpolator_t *interpolator);
int resampler_init(resampler_t *resampler,
float32_t input_rate,
float32_t output_rate,
uint16_t filter_order);
float32_t resampler_process(resampler_t *resampler, float32_t input);
void resampler_deinit(resampler_t *resampler);
// 高效的多相实现
int polyphase_filter_init(polyphase_filter_t *filter,
const float32_t *coeffs,
uint16_t num_coeffs,
uint16_t num_phases);
float32_t polyphase_filter_process(polyphase_filter_t *filter, float32_t input);
void polyphase_filter_deinit(polyphase_filter_t *filter);
#endif // MULTIRATE_FILTER_H
9. 窗函数 (include/window_functions.h)
#ifndef WINDOW_FUNCTIONS_H
#define WINDOW_FUNCTIONS_H
#include <stdint.h>
#include <stddef.h>
// 窗函数类型
typedef enum {
WINDOW_RECTANGULAR = 0,
WINDOW_HANN,
WINDOW_HAMMING,
WINDOW_BLACKMAN,
WINDOW_KAISER,
WINDOW_GAUSSIAN,
WINDOW_TRIANGULAR
} window_type_t;
// 函数声明
void apply_window(float32_t *data, uint16_t length, window_type_t window);
void rectangular_window(float32_t *data, uint16_t length);
void hann_window(float32_t *data, uint16_t length);
void hamming_window(float32_t *data, uint16_t length);
void blackman_window(float32_t *data, uint16_t length);
void kaiser_window(float32_t *data, uint16_t length, float32_t beta);
void gaussian_window(float32_t *data, uint16_t length, float32_t sigma);
// 计算窗函数的参数
float32_t calculate_kaiser_beta(float32_t ripple_db);
uint16_t calculate_filter_order(float32_t transition_width, float32_t sampling_freq, window_type_t window);
#endif // WINDOW_FUNCTIONS_H
10. 窗函数实现 (src/window_functions.c)
#include "window_functions.h"
#include "math_utils.h"
#include <math.h>
// 应用窗函数
void apply_window(float32_t *data, uint16_t length, window_type_t window) {
if (!data || length == 0) {
return;
}
switch (window) {
case WINDOW_RECTANGULAR:
rectangular_window(data, length);
break;
case WINDOW_HANN:
hann_window(data, length);
break;
case WINDOW_HAMMING:
hamming_window(data, length);
break;
case WINDOW_BLACKMAN:
blackman_window(data, length);
break;
case WINDOW_KAISER:
kaiser_window(data, length, 8.6f); // 默认beta值
break;
case WINDOW_GAUSSIAN:
gaussian_window(data, length, 0.5f); // 默认sigma
break;
default:
rectangular_window(data, length);
break;
}
}
// 矩形窗
void rectangular_window(float32_t *data, uint16_t length) {
for (uint16_t i = 0; i < length; i++) {
data[i] = 1.0f;
}
}
// 汉宁窗
void hann_window(float32_t *data, uint16_t length) {
for (uint16_t i = 0; i < length; i++) {
float32_t factor = 2.0f * PI_F * i / (length - 1);
data[i] = 0.5f * (1.0f - cosf(factor));
}
}
// 汉明窗
void hamming_window(float32_t *data, uint16_t length) {
for (uint16_t i = 0; i < length; i++) {
float32_t factor = 2.0f * PI_F * i / (length - 1);
data[i] = 0.54f - 0.46f * cosf(factor);
}
}
// 布莱克曼窗
void blackman_window(float32_t *data, uint16_t length) {
for (uint16_t i = 0; i < length; i++) {
float32_t factor = 2.0f * PI_F * i / (length - 1);
data[i] = 0.42f - 0.5f * cosf(factor) + 0.08f * cosf(2.0f * factor);
}
}
// 凯泽窗
void kaiser_window(float32_t *data, uint16_t length, float32_t beta) {
float32_t i0_beta = bessel_i0(beta);
for (uint16_t i = 0; i < length; i++) {
float32_t arg = beta * sqrtf(1.0f - powf(2.0f * i / (length - 1) - 1.0f, 2.0f));
data[i] = bessel_i0(arg) / i0_beta;
}
}
// 高斯窗
void gaussian_window(float32_t *data, uint16_t length, float32_t sigma) {
float32_t scale = -0.5f / (sigma * sigma);
for (uint16_t i = 0; i < length; i++) {
float32_t x = (i - length / 2.0f) / (length / 2.0f);
data[i] = expf(scale * x * x);
}
}
// 贝塞尔函数I0(用于凯泽窗)
float32_t bessel_i0(float32_t x) {
float32_t sum = 1.0f;
float32_t term = 1.0f;
float32_t half_x = x / 2.0f;
for (uint16_t k = 1; k <= 20; k++) {
term *= (half_x * half_x) / (k * k);
sum += term;
}
return sum;
}
// 计算凯泽窗的beta参数
float32_t calculate_kaiser_beta(float32_t ripple_db) {
if (ripple_db > 50) {
return 0.1102f * (ripple_db - 8.7f);
} else if (ripple_db >= 21) {
return 0.5842f * powf(ripple_db - 21.0f, 0.4f) + 0.07886f * (ripple_db - 21.0f);
} else {
return 0.0f;
}
}
// 计算滤波器阶数
uint16_t calculate_filter_order(float32_t transition_width, float32_t sampling_freq, window_type_t window) {
float32_t normalized_transition = transition_width / sampling_freq;
switch (window) {
case WINDOW_RECTANGULAR:
return (uint16_t)(1.8f / normalized_transition);
case WINDOW_HANN:
return (uint16_t)(6.2f / normalized_transition);
case WINDOW_HAMMING:
return (uint16_t)(6.6f / normalized_transition);
case WINDOW_BLACKMAN:
return (uint16_t)(11.0f / normalized_transition);
case WINDOW_KAISER:
return (uint16_t)(10.0f / normalized_transition);
default:
return (uint16_t)(6.2f / normalized_transition);
}
}
11. 示例程序 (examples/fir_lowpass_example.c)
#include <stdio.h>
#include <stdlib.h>
#include "fir_filter.h"
#include "window_functions.h"
#define SAMPLE_RATE 48000.0f
#define CUTOFF_FREQ 5000.0f
#define NUM_TAPS 101
int main() {
printf("FIR低通滤波器示例\n");
printf("采样率: %.0f Hz\n", SAMPLE_RATE);
printf("截止频率: %.0f Hz\n", CUTOFF_FREQ);
printf("抽头数: %d\n\n", NUM_TAPS);
// 创建滤波器
fir_filter_t filter;
fir_config_t config;
// 设计滤波器
if (fir_design_lowpass(&filter, CUTOFF_FREQ, SAMPLE_RATE, NUM_TAPS, WINDOW_HAMMING) != 0) {
printf("滤波器设计失败!\n");
return -1;
}
printf("滤波器设计成功!\n");
printf("滤波器系数:\n");
// 打印前10个系数
for (int i = 0; i < 10; i++) {
printf("h[%d] = %.6f\n", i, filter.coeffs[i]);
}
printf("...\n");
// 测试信号:包含1kHz和10kHz的正弦波
printf("\n测试信号处理:\n");
printf("输入: 1kHz + 10kHz 正弦波\n");
printf("期望: 只保留1kHz成分\n\n");
float32_t input_signal[1000];
float32_t output_signal[1000];
// 生成测试信号
for (int i = 0; i < 1000; i++) {
float32_t t = (float32_t)i / SAMPLE_RATE;
input_signal[i] = sinf(2.0f * PI_F * 1000.0f * t) + // 1kHz
sinf(2.0f * PI_F * 10000.0f * t); // 10kHz
}
// 处理信号
fir_process_block(&filter, input_signal, output_signal, 1000);
// 显示前10个输出样本
printf("前10个输出样本:\n");
for (int i = 0; i < 10; i++) {
printf("out[%d] = %.6f\n", i, output_signal[i]);
}
// 计算频率响应
printf("\n频率响应(幅度):\n");
for (float32_t freq = 0; freq <= SAMPLE_RATE/2; freq += 1000) {
float32_t response = 0.0f;
for (int i = 0; i < filter.num_taps; i++) {
float32_t phase = 2.0f * PI_F * freq * i / SAMPLE_RATE;
response += filter.coeffs[i] * cosf(phase);
}
printf("%.0f Hz: %.4f\n", freq, response);
}
// 清理资源
fir_deinit(&filter);
return 0;
}
12. Makefile
CC = gcc
CFLAGS = -Wall -O2 -std=c99 -fPIC
AR = ar
ARFLAGS = rcs
# 目录
INCLUDE_DIR = include
SRC_DIR = src
EXAMPLES_DIR = examples
TESTS_DIR = tests
BUILD_DIR = build
LIB_DIR = lib
# 目标
LIBRARY = $(LIB_DIR)/libdigitalfilters.a
SHARED_LIB = $(LIB_DIR)/libdigitalfilters.so
# 源文件
SRCS = $(wildcard $(SRC_DIR)/*.c)
OBJS = $(SRCS:$(SRC_DIR)/%.c=$(BUILD_DIR)/%.o)
# 示例程序
EXAMPLES = fir_lowpass_example iir_bandpass_example lms_adaptive_example decimation_example
.PHONY: all clean directories library examples test
all: directories library examples
directories:
mkdir -p $(BUILD_DIR) $(LIB_DIR)
library: $(LIBRARY) $(SHARED_LIB)
$(LIBRARY): $(OBJS)
$(AR) $(ARFLAGS) $@ $^
$(SHARED_LIB): $(OBJS)
$(CC) -shared -o $@ $^
$(BUILD_DIR)/%.o: $(SRC_DIR)/%.c
$(CC) $(CFLAGS) -I$(INCLUDE_DIR) -c $< -o $@
examples: $(EXAMPLES)
fir_lowpass_example: $(BUILD_DIR)/fir_lowpass_example.o $(LIBRARY)
$(CC) -o $@ $^ -lm
iir_bandpass_example: $(BUILD_DIR)/iir_bandpass_example.o $(LIBRARY)
$(CC) -o $@ $^ -lm
lms_adaptive_example: $(BUILD_DIR)/lms_adaptive_example.o $(LIBRARY)
$(CC) -o $@ $^ -lm
decimation_example: $(BUILD_DIR)/decimation_example.o $(LIBRARY)
$(CC) -o $@ $^ -lm
$(BUILD_DIR)/%.o: $(EXAMPLES_DIR)/%.c
$(CC) $(CFLAGS) -I$(INCLUDE_DIR) -c $< -o $@
test: $(TESTS_DIR)/test_fir.c $(LIBRARY)
$(CC) $(CFLAGS) -I$(INCLUDE_DIR) -o test_fir $(TESTS_DIR)/test_fir.c $(LIBRARY) -lm
./test_fir
clean:
rm -rf $(BUILD_DIR) $(LIB_DIR) $(EXAMPLES) test_fir
install: library
cp $(LIBRARY) /usr/local/lib/
cp $(SHARED_LIB) /usr/local/lib/
cp -r $(INCLUDE_DIR)/*.h /usr/local/include/
uninstall:
rm -f /usr/local/lib/libdigitalfilters.*
rm -f /usr/local/include/fir_filter.h
rm -f /usr/local/include/iir_filter.h
rm -f /usr/local/include/adaptive_filter.h
rm -f /usr/local/include/multirate_filter.h
rm -f /usr/local/include/window_functions.h
rm -f /usr/local/include/math_utils.h
help:
@echo "可用目标:"
@echo " all - 编译所有内容"
@echo " library - 编译静态和动态库"
@echo " examples - 编译示例程序"
@echo " test - 运行测试"
@echo " clean - 清理编译文件"
@echo " install - 安装库和头文件"
@echo " uninstall - 卸载库和头文件"
参考代码 数字滤波器的C语言实现 www.youwenfan.com/contentcsv/112860.html
三、使用说明
1. 基本用法
#include "fir_filter.h"
int main() {
// 1. 创建滤波器
fir_filter_t filter;
fir_config_t config;
// 2. 设计低通滤波器
fir_design_lowpass(&filter, 1000.0f, 8000.0f, 51, WINDOW_HAMMING);
// 3. 处理信号
float32_t input = 1.0f;
float32_t output = fir_process(&filter, input);
// 4. 清理
fir_deinit(&filter);
return 0;
}
2. 滤波器类型对比
| 滤波器类型 | 优点 | 缺点 | 适用场景 |
|---|---|---|---|
| FIR | 线性相位、稳定、设计简单 | 阶数高、延迟大 | 音频处理、图像处理 |
| IIR | 阶数低、效率高 | 非线性相位、可能不稳定 | 实时处理、嵌入式 |
| 自适应 | 自动跟踪变化 | 收敛慢、计算复杂 | 噪声消除、回声消除 |
| 多速率 | 改变采样率 | 设计复杂 | 采样率转换、多速率系统 |
3. 性能优化建议
- 使用对称FIR:减少50%乘法运算
- 直接形式II:减少内存使用
- 定点运算:使用Q格式提高嵌入式性能
- SIMD指令:使用NEON/SSE加速
- 查表法:预先计算三角函数等
4. 嵌入式移植
// 嵌入式版本(无动态内存分配)
static float32_t coeffs[51];
static float32_t buffer[51];
static fir_filter_t filter = {
.coeffs = coeffs,
.buffer = buffer,
.num_taps = 51,
.buffer_index = 0,
.initialized = 1
};
// 初始化
void init_filter(void) {
// 手动设置系数
coeffs[0] = 0.001f;
// ... 设置其他系数
}
// 处理
float32_t process_sample(float32_t input) {
return fir_process(&filter, input);
}
更多推荐



所有评论(0)