四阶IIR滤波器C语言代码实现

Oliver-H

已于 2024-09-19 17:52:19 修改

阅读量674

点赞数 3

分类专栏： Matlab数字信号处理文章标签： matlab 算法开发语言 c++

于 2024-09-19 09:40:13 首次发布

本文链接：https://blog.csdn.net/yishuihanq/article/details/142351392

版权

Matlab数字信号处理专栏收录该内容

91 篇文章

订阅专栏

用matlab的filterDesigner设计一个四阶高通滤波器

% %
% % IIR butterworth滤波器系数% 高通（直接II型4阶、Fs = 20KHz、Fc = 1500Hz）PLC=0.0025PLC
% % 分子、分母系数（对累加和进行滤波）
b1 = [0.5365885311137,   -2.146354124455,    3.219531186682,   -2.146354124455, 0.5365885311137];
a1 = [1,   -2.773682317549,    3.019038699424,   -1.504765051425, 0.2879304294211];
% % 阶跃响应
% figure
% sys = tf(b1,a1);
% step(sys)

cur_flt1 = filter(b1,a1,cur_diffsum); % 滤波处理
vol_flt1 = filter(b1,a1,vol_diffsum); % 滤波处理
figure("Name", "IIR cur_4jie1500")
plot(cur_flt1,'*');
figure("Name", "IIR vol_4jie1500")
plot(vol_flt1,'*');

filter的实际就是对X数组做离散运算得到滤波后数组；

使用的离散序列为：

a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + b(3)*x(n-2) + ... + b(nb+1)*x(n-nb)
- a(2)*y(n-1) - a(3)*y(n-2) - ... - a(na+1)*y(n-na)


%%%===================================%%%
x = cur_diffsum;
y(1)=b1(1)*x(1);
y(2)=b1(1)*x(2)+b1(2)*x(1)-a1(2)*y(1);
y(3)=b1(1)*x(3)+b1(2)*x(2)+b1(3)*x(1)-a1(2)*y(2)-a1(3)*y(1);
y(4)=b1(1)*x(4)+b1(2)*x(3)+b1(3)*x(2)+b1(4)*x(1)-a1(2)*y(3)-a1(3)*y(2)-a1(4)*y(1);
for i=5:len
    y(i) = b1(1)*x(i)+b1(2)*x(i-1)+b1(3)*x(i-2)+b1(4)*x(i-3)+b1(5)*x(i-4) - a1(2)*y(i-1)-a1(3)*y(i-2)-a1(4)*y(i-3)-a1(5)*y(i-4);
end
figname = figure("Name", "cur_fft_y(i)");
signal_fft(y, Fs, 256, 1, figname);
%%%===================================%%%

%%%+++++++++++++++++++++++++++++++++++%%%
SumY = 0;
SumX = 0;
FltX = [0,0,0,0,0];
FltY = [0,0,0,0,0];
for i=1:len
    FltX(5) = FltX(4);
    FltX(4) = FltX(3);
    FltX(3) = FltX(2);
    FltX(2) = FltX(1);
    FltX(1) = x(i);

    for p=1:5
        SumX = SumX + b1(p)*FltX(p);
    end
    for q=2:5
        SumY = SumY + a1(q)*FltY(q-1);
    end
    FltY(4) = FltY(3);
    FltY(3) = FltY(2);
    FltY(2) = FltY(1);

    FltY(1) = (SumX - SumY);
    Yn(i) = FltY(1);
    SumY = 0;
    SumX = 0;
end
figname = figure("Name", "cur_fft_Yn");
signal_fft(Yn, Fs, 256, 1, figname);

其C语言实现形式如下：

//% % IIR butterworth滤波器系数% 高通（直接II型4阶、Fs = 20KHz、Fc = 1500Hz）
//% % 分子、分母系数（对累加和进行滤波）
double b1[5] = { 0.5365885311137, -2.146354124455, 3.219531186682, -2.146354124455, 0.5365885311137 };
double a1[5] = { 1, -2.773682317549, 3.019038699424, -1.504765051425, 0.2879304294211 };

/**
* 功能:高通滤波
*      
* 参数:
*      Xn: 传入的 diffsum 值
* 返回值:
*      Yn: 滤波后的值
*/
double iir_filter_diffsum(double Xn)
{
	// 初始化滤波参数
	static double FltX[5] = { 0, 0, 0, 0, 0 };
	static double FltY[5] = { 0, 0, 0, 0, 0 };
	double SumX = 0, SumY = 0, Yn = 0;
	int i = 0;

    FltX[4] = FltX[3];
    FltX[3] = FltX[2];
    FltX[2] = FltX[1];
    FltX[1] = FltX[0];
    FltX[0] = Xn;

    for (i = 0; i < 5; i++)
    {
        SumX = SumX + b1[i] * FltX[i];
    }
    for (i = 1; i < 5; i++)
    {
        SumY = SumY + a1[i] * FltY[i - 1];
    }
    FltY[3] = FltY[2];
    FltY[2] = FltY[1];
    FltY[1] = FltY[0];
    FltY[0] = (SumX - SumY);
   
    Yn = FltY[0];

    return Yn;
}

滑动均值滤波：

/**
* 功能:滑动求平均
*
* 输入参数:
*      data: 传入的数据
*      d_size: data数据长度
*      w_size: 滑动窗口大小
* 输出参数:
*      result: 滤波后的值
*/
void moving_average(double *data, double* result, int d_size, int w_size)
{
	for (int i = 0; i < d_size - w_size + 1; i++)
	{
		double sum = 0;
		for (int j = 0; j < w_size;j++)
		{
			sum += data[i + j];
		}
		//result[i + (w_size - 1) / 2] = sum / w_size;
		result[i] = sum / w_size;
	}

	//将末尾的数据直接添加到平均值数组中
	for (int i = d_size - w_size + 1; i < d_size;i++)
	{
		result[i] = data[i];
	}
}