IIR
最近在写的一个程序需要用到
以前用到的
~fisher/mkfilter/trad.html
一直都偷懒直接用这个网站的结果，所以手上也没积累什么有用的代码。这次就需要自己从头做起。
我面临的问题有两个：
1.根据滤波器的参数（滤波器的类型、截止频率、滤波器的阶数等），计算出滤波器对应的差分方程的系数。
2.利用得到的差分方程的系数构造一个可以工作的滤波器。
其中第一个问题，对于不同类型的滤波器，比如
这里就先写写第二个问题。IIR滤波器对应的差分方程为：
相应的系统函数为：
这里默认
按照奥本海默写的《离散时间信号处理》上面的介绍，
直接I型class IIR_I{private:double *m_pNum;double *m_pDen;double *m_px;double *m_py;int m_num_order;int m_den_order;public:IIR_I();void reset();void setPara(double num[], int num_order, double den[], int den_order);void resp(double data_in[], int m, double data_out[], int n);double filter(double data);void filter(double data[], int len);void filter(double data_in[], double data_out[], int len);};
其中
三个
下面是实现代码：
/** \brief 将滤波器的内部状态清零，滤波器的系数保留 * \return */void IIR_I::reset(){for(int i = 0; i <= m_num_order; i++){m_pNum[i] = 0.0;}for(int i = 0; i <= m_den_order; i++){m_pDen[i] = 0.0;}}IIR_I::IIR_I(){m_pNum = NULL;m_pDen = NULL;m_px = NULL;m_py = NULL;m_num_order = -1;m_den_order = -1;};/** \brief * * \param num 分子多项式的系数，升序排列,num[0] 为常数项 * \param m 分子多项式的阶数 * \param den 分母多项式的系数，升序排列,den[0] 为常数项 * \param m 分母多项式的阶数 * \return */void IIR_I::setPara(double num[], int num_order, double den[], int den_order){delete[] m_pNum;delete[] m_pDen;delete[] m_px;delete[] m_py;m_pNum = new double[num_order + 1];m_pDen = new double[den_order + 1];m_num_order = num_order;m_den_order = den_order;m_px = new double[num_order + 1];m_py = new double[den_order + 1];for(int i = 0; i <= m_num_order; i++){m_pNum[i] = num[i];m_px[i] = 0.0;}for(int i = 0; i <= m_den_order; i++){m_pDen[i] = den[i];m_py[i] = 0.0;}}/** \brief 滤波函数，采用直接I型结构 * * \param data 传入输入数据 * \return 滤波后的结果 */double IIR_I::filter(double data){m_py[0] = 0.0; // 存放滤波后的结果m_px[0] = data;for(int i = 0; i <= m_num_order; i++){m_py[0] = m_py[0] + m_pNum[i] * m_px[i];}for(int i = 1; i <= m_den_order; i++){m_py[0] = m_py[0] – m_pDen[i] * m_py[i];}for(int i = m_num_order; i >= 1; i–){m_px[i] = m_px[i-1];}for(int i = m_den_order; i >= 1; i–){m_py[i] = m_py[i-1];}return m_py[0];}/** \brief 滤波函数，采用直接I型结构 * * \param data[] 传入输入数据，返回时给出滤波后的结果 * \param len data[] 数组的长度 * \return */void IIR_I::filter(double data[], int len){int i, k;for(k = 0; k < len; k++){m_px[0] = data[k];data[k] = 0.0;for(i = 0; i <= m_num_order; i++){data[k] = data[k] + m_pNum[i] * m_px[i];}for(i = 1; i <= m_den_order; i++){data[k] = data[k] – m_pDen[i] * m_py[i];}// we get the y value nowm_py[0] = data[k];for(i = m_num_order; i >= 1; i–){m_px[i] = m_px[i-1];}for(i = m_den_order; i >= 1; i–){m_py[i] = m_py[i-1];}}}/** \brief 滤波函数，采用直接I型结构 * * \param data_in[] 输入数据 * \param data_out[] 保存滤波后的数据 * \param len 数组的长度 * \return */void IIR_I::filter(double data_in[], double data_out[], int len){int i, k;for(k = 0; k < len; k++){m_px[0] = data_in[k];m_py[0] = 0.0;for(i = 0; i <= m_num_order; i++){m_py[0] = m_py[0] + m_pNum[i] * m_px[i];}for(i = 1; i <= m_den_order; i++){m_py[0] = m_py[0] – m_pDen[i] * m_py[i];}for(i = m_num_order; i >= 1; i–){m_px[i] = m_px[i-1];}for(i = m_den_order; i >= 1; i–){m_py[i] = m_py[i-1];}data_out[k] = m_py[0];}}
除此之外，还提供了个
/** \brief 计算 IIR 滤波器的时域响应，不影响滤波器的内部状态 * \param data_in 为滤波器的输入，0 时刻之前的输入默认为 0，data_in[M] 及之后的输入默认为data_in[M-1] * \param data_out 滤波器的输出 * \param M 输入数据的长度 * \param N 输出数据的长度 * \return */void IIR_I::resp(double data_in[], int M, double data_out[], int N){int i, k, il;for(k = 0; k < N; k++){data_out[k] = 0.0;for(i = 0; i <= m_num_order; i++){if( k – i >= 0){il = ((k – i) < M) ? (k – i) : (M – 1);data_out[k] = data_out[k] + m_pNum[i] * data_in[il];}}for(i = 1; i <= m_den_order; i++){if( k – i >= 0){data_out[k] = data_out[k] – m_pDen[i] * data_out[k – i];}}}}直接II型/**< IIR 滤波器直接II型实现 */class IIR_II{public:IIR_II();void reset();void setPara(double num[], int num_order, double den[], int den_order);void resp(double data_in[], int m, double data_out[], int n);double filter(double data);void filter(double data[], int len);void filter(double data_in[], double data_out[], int len);protected:private:double *m_pNum;double *m_pDen;double *m_pW;int m_num_order;int m_den_order;int m_N;};class IIR_BODE{private:double *m_pNum;double *m_pDen;int m_num_order;int m_den_order;complex<double> poly_val(double p[], int order, double omega);public:IIR_BODE();void setPara(double num[], int num_order, double den[], int den_order);complex<double> bode(double omega);void bode(double omega[], int n, complex<double> resp[]);};
直接II型实现中所需的存储单元少了很多。另外，这两种实现的外部接口是完全相同的。
细数门前落叶，倾听窗外雨声，
