滤波器数学模型
如图Fig1所示，是IIR二阶数字滤波器的数学计算公式
Fig1.png
转换到离散域，计算公式如下
Fig2.png
其中y(n)为当前采样点滤波后的数值，y(n-1)为上一个采样点滤波后的数值,y(n-2)类推
x(n)为当前采样点滤波前的数值，x(n-1),x(n-2)依次类推
滤波器系数计算
滤波器主要有以下几种：高通/低通/带通。下分别说明此三种滤波器的系统的求取方法。
通常，对一个滤波器的要求，我们主要给出以下技术规格：中心频率frequency，采样频率sampleRate，增益dBgain，品质因数Q。
为计算方便，先定义以下几个值：
Fig3.jpg
高通滤波器系数的计算：
Fig4.jpg
低通滤波器系数的计算：
Fig5.jpg
带通滤波器系统计算：
Fig6.jpg
采样及信号处理
Fig7.png
采样原始信号为一个正弦信号和一个直流分量的相加。直流分量的幅值为1，正弦信号周期为1s，幅值也为1。
初始化采样信号
#define NUM_POINTS 1000 // 采样点1000个
#define PI 3.1415926
double m_signal[NUM_POINTS ]; // 原始采样信号
double m_signal_filter[NUM_POINTS ]; // 滤波后信号
double delta = 0;
for (int i = 0; i < NUM_POINTS; i++)
{
m_signal[i] = 1+sin(2.0*PI*delta); // 直流+正弦信号
delta += 1e-3; // 采样间隔0.001s
}
首先定义数字滤波器的结构体
typedef struct
{
UINT16 type;
float fl;
}FILTERTYPE; // 滤波器类型
typedef struct
{
FILTERTYPE type;
float ts; // 采样间隔时间
float dbgain; // 增益
float q; // 品质因子
float amp; // 幅值
float a[3]; // a0,a1,a2
float b[3]; // b0,b1,b2
float x[3]; // x(n),x(n-1),x(n-2)
float y[2]; // y(n-1),y(n-2)
}FILTERINFO, *LPFILTERINFO;
初始化滤波器，这里我们选择二阶带通滤波器做测试
#define FILTER_BYPASS 0 // 不滤波
#define FILTER_1ST_LP 1 // 一阶低通
#define FILTER_1ST_HP 2 // 一阶高通
#define FILTER_2ND_LP 3 // 2阶低通
#define FILTER_2ND_HP 4 // 2阶高通
#define FILTER_2ND_BP 5 // 2阶带通
#define FILTER_2ND_BS 6 // 2阶带阻
FILTERINFO pre_filter_Vo;
memset(&pre_filter_Vo, 0, sizeof(FILTERINFO));
// 信号的增益不予改变，所以不去计算dbgain
pre_filter_Vo.type.type = FILTER_2ND_BP; // 二阶带通滤波器
pre_filter_Vo.ts = 1.0 / NUM_POINTS;
pre_filter_Vo.type.fl = 1; // 中心频率为1Hz
pre_filter_Vo.q = 1; // 品质因子为1，对中心频率信号滤波效果最佳
Init_Filter(&pre_filter_Vo, 1);
其中，Init_Filter函数定义如下
void CMonitorDlg::Init_Filter(LPFILTERINFO fi, UINT16 zero)
{
double omega = 2 * PI*fi->type.fl*fi->ts;
double sine = sin(omega);
double cosine = cos(omega);
double alpha = sine / (2 * fi->q);
fi->b[0] = fi->b[1] = fi->b[2] = fi->a[0] = fi->a[1] = fi->a[2] = 0;
switch (fi->type.type)
{
case FILTER_1ST_LP:
case FILTER_1ST_HP:
fi->b[0] = fi->type.fl * 2.0 * PI*fi->ts;
fi->a[1] = 1 - fi->b[0];
break;
case FILTER_2ND_LP:
fi->a[0] = 1 + alpha;
fi->a[1] = (-2 * cosine) / fi->a[0];
fi->a[2] = (1 - alpha) / fi->a[0];
fi->b[0] = ((1 - cosine) / 2) / fi->a[0];
fi->b[1] = (1 - cosine) / fi->a[0];
fi->b[2] = ((1 - cosine) / 2) / fi->a[0];
break;
case FILTER_2ND_HP:
fi->a[0] = 1 + alpha;
fi->a[1] = (-2 * cosine) / fi->a[0];
fi->a[2] = (1 - alpha) / fi->a[0];
fi->b[0] = ((1 + cosine) / 2) / fi->a[0];
fi->b[1] = (-(1 + cosine)) / fi->a[0];
fi->b[2] = ((1 + cosine) / 2) / fi->a[0];
break;
case FILTER_2ND_BP:
case FILTER_2ND_BS:
fi->a[0] = 1 + alpha;
fi->a[1] = (-2 * cosine) / fi->a[0];
fi->a[2] = (1 - alpha) / fi->a[0];
fi->b[0] = (sine / 2) / fi->a[0];
fi->b[1] = 0;
fi->b[2] = (-sine / 2) / fi->a[0];
break;
}
if (zero)
fi->x[0] = fi->x[1] = fi->x[2] = fi->y[0] = fi->y[1] = 0;
}
DSP滤波器函数
float Filter(LPFILTERINFO fi, float in)
{
float f = in;
if (fi->type.type == FILTER_1ST_LP || fi->type.type == FILTER_1ST_HP)
{
f = fi->b[0] * in + fi->a[1] * fi->y[0];
fi->y[0] = f;
}
else if (fi->type.type)
{
fi->x[2] = fi->x[1];
fi->x[1] = fi->x[0];
fi->x[0] = in;
f = fi->b[0] * fi->x[0] + fi->b[1] * fi->x[1] + fi->b[2] * fi->x[2] - fi->a[1] * fi->y[0] - fi->a[2] * fi->y[1];
fi->y[1] = fi->y[0];
fi->y[0] = f;
}
return f;
}
使用滤波器函数滤波原始采样信号
for (int i = 0; i < NUM_POINTS; i++)
{
m_signal_filter[i] = Filter(&pre_filter_Vo, m_signal[i]); // 滤波
}
绘制信号
m_signal采样信号如下，由于包含直流分量1，所以它的幅值由正弦信号的（-1,1）沿Y轴向上偏移1位变成（0,2）
Fig8.png
滤波后的信号m_signal_filter如下所示
Fig9.png
很明显，通过中心频率为1Hz，品质因子为1的二阶带通滤波器，采样信号的直流分量被完全隔离掉，只剩下原始正弦信号。
结论和遗留问题
通过上述实验验证了我们给出的二阶带通滤波器算法的有效性，但是我们滤波后的信号其实并不是真正的正弦信号，而是一个无限逼近正弦的信号，这是由于以下几个原因造成的
1）采样时间：采样速率越高，信号失真越小
2）计算精度：计算机处理信号时在存储浮点数会出现一定误差
3）算法精度：算法本身就有精度误差
改进：
1）读者们可以自行提高采样速率，去验证采样速率对信号误差的影响，评判标准可以使用均方根误差
2）读者们还可以通过更改中心频率，品质因子，选择不同的滤波器类型测试对信号的影响
