C语言IIR双向滤波

设计一个0.5~1Hz的IIR滤波器，用巴特沃斯或者契比雪夫2，看零极点图是稳定的。

设计如下：

function Hd = iir_highpass_05_1_buter
%IIR_HIGHPASS_05_1_BUTER Returns a discrete-time filter object.

% MATLAB Code
% Generated by MATLAB(R) 9.5 and DSP System Toolbox 9.7.
% Generated on: 27-Nov-2023 20:11:44

% Butterworth Highpass filter designed using FDESIGN.HIGHPASS.

% All frequency values are in Hz.
Fs = 250;  % Sampling Frequency

Fstop = 0.5;         % Stopband Frequency
Fpass = 1;           % Passband Frequency
Astop = 30;          % Stopband Attenuation (dB)
Apass = 1;           % Passband Ripple (dB)
match = 'stopband';  % Band to match exactly

% Construct an FDESIGN object and call its BUTTER method.
h  = fdesign.highpass(Fstop, Fpass, Astop, Apass, Fs);
Hd = design(h, 'butter', 'MatchExactly', match);
% % 
% % Get the transfer function values.
% [b, a] = tf(Hd);
% b 
% a
% close all
% figure
% freqz(b,a);
% figure
% zplane(b,a);
% [H,w]=freqz(b,a);
% % [H,w]=freqz(b,a,5000);
% N = length(H);
% Hf=abs(H);  %取幅度值实部
% Hx=angle(H);  %取相位值对应相位角
% % clf
% figure
% % plot(w,20*log10(Hf))  %幅值变换为分贝单位
% plot([0:1/N:1-1/N]*Fs/2,20*log10(Hf))  %幅值变换为分贝单位
% title('离散系统幅频特性曲线')
% figure
% % plot(w,Hx)
% plot([0:1/N:1-1/N]*Fs/2,Hx);
% title('离散系统相频特性曲线')

% [EOF]

或：

function Hd = iir_highpass_05_1_cheby2
%IIR_HIGHPASS_05_1_BUTER Returns a discrete-time filter object.

% MATLAB Code
% Generated by MATLAB(R) 9.5 and DSP System Toolbox 9.7.
% Generated on: 27-Nov-2023 20:10:02

% Chebyshev Type II Highpass filter designed using FDESIGN.HIGHPASS.

% All frequency values are in Hz.
Fs = 250;  % Sampling Frequency

Fstop = 0.5;         % Stopband Frequency
Fpass = 1;           % Passband Frequency
Astop = 30;          % Stopband Attenuation (dB)
Apass = 1;           % Passband Ripple (dB)
match = 'stopband';  % Band to match exactly

% Construct an FDESIGN object and call its CHEBY2 method.
h  = fdesign.highpass(Fstop, Fpass, Astop, Apass, Fs);
Hd = design(h, 'cheby2', 'MatchExactly', match);
% 
% % Get the transfer function values.
% [b, a] = tf(Hd);
% b 
% a
% close all
% figure
% freqz(b,a);
% figure
% zplane(b,a);
% % [H,w]=freqz(b,a);
% [H,w]=freqz(b,a,5000);
% N = length(H);
% Hf=abs(H);  %取幅度值实部
% Hx=angle(H);  %取相位值对应相位角
% % clf
% figure
% % plot(w,20*log10(Hf))  %幅值变换为分贝单位
% plot([0:1/N:1-1/N]*Fs/2,20*log10(Hf))  %幅值变换为分贝单位
% title('离散系统幅频特性曲线')
% figure
% % plot(w,Hx)
% plot([0:1/N:1-1/N]*Fs/2,Hx);
% title('离散系统相频特性曲线')

% [EOF]

可以得到滤波器系数a,b.

设计一个通用的双向滤波程序，包括滤波函数和数组逆序函数,包含可以和matlab对数据的写数据代码。主要有：

main.c

#include"stdio.h"
#include"stdlib.h"
#include"math.h"


#include"highpass_iir_05_1.h"


#define M_PI 3.1415
#define N 200
#define F1 3
#define F2 40

//float fir_coef_b[5] = { 0.979965406517890, - 3.91970687664644,	5.87948294331178, - 3.91970687664644,	0.979965406517890 };
//float fir_coef_a[5] = { 1, - 3.95936787226356,	5.87908160514721, - 3.88004583425791,	0.960332197971774 };

float fir_coef_b[5] = { 1,2,1,3,1 };
float fir_coef_a[5] = { 1,1,1,1,1 };

int  main()
{
    // 生成正弦信号
    float a[N] = { 0 }, b[N] = { 0 };
    for (int i = 0; i < N; i++)
    {
        a[i] = (float)(sin(2 * M_PI * i * F1 / N) + 0.25 * sin(2 * M_PI * i * F2 / N));
    }

    FILE* fp;

    fp = fopen("text.txt", "w");
    if (feof(fp))
    {
        printf("NULL");
        exit(0);//表示如果读取为空文件就正常退出
    }

    for (int i = 0; i < N; i++)
        fprintf(fp, "%.6f\n", a[i]);
    fclose(fp);
    

    int m;

    // 双向滤波
    m = iir_filter_zhh(a, b, N, fir_coef_a, 5, fir_coef_b,5);
    reverse_arry_float(b, N);
    memcpy(a, b, N * sizeof(float));
    m = iir_filter_zhh(a, b, N, fir_coef_a, 5, fir_coef_b, 5);
    reverse_arry_float(b, N);

    fp = fopen("text_fir.txt", "w");
    if (feof(fp))
    {
        printf("NULL");
        exit(0);//表示如果读取为空文件就正常退出
    }

    for (int i = 0; i < N; i++)
        fprintf(fp, "%.6f\n", b[i]);
    fclose(fp);


    return 0;
}

highpass_iir_05_1.c

#include"stdio.h"
#include"stdlib.h"
#include"string.h"

#include "assert.h" 


// MATLAB函数
//function y = iir_filter_zhh(b, a, x)
//
//len_x = length(x);
//y = zeros(size(x));
//
//len_b = length(b);
//len_a = length(a);
//x_buf = [zeros(len_b - 1, 1); x];
//y_buf = [zeros(len_a - 1, 1); y];
//
//b_fir = fliplr(b);
//a_fir = fliplr(a);
//
//for i = 1:len_x
//y(i) = 1 / a_fir(end) * (b_fir * x_buf(i:i + len_b - 1) - a_fir(1:end - 1) * y_buf(i:i + len_a - 2));
//y_buf(i + len_a - 1) = y(i);
//end
//
//end

int iir_filter_zhh(float* sig_in, float* sig_out, int sig_len, float* a, int a_len, float* b, int b_len)
{
	// 参数个数检查,至少3个至多5个参数.系统自行检查，参数个数不对及类型不对会自动报错。
	assert(sig_len > 1);
	assert(a_len > 1);
	assert(b_len > 1);


	int i, j;
	float tmp1, tmp2;
	float* sig_in_buf = (float*)malloc((size_t)(sig_len + b_len - 1) * sizeof(float));
	if (sig_in_buf == NULL) {//判空
		printf("sig_in_buf malloc error");//打印错误信息
		return 1;
	}
	memset(sig_in_buf, 0, (size_t)(b_len - 1) * sizeof(float));
	memcpy(sig_in_buf + b_len - 1, sig_in, (size_t)sig_len * sizeof(float));

	float* sig_out_buf = (float*)malloc((size_t)(sig_len + a_len - 1) * sizeof(float));
	if (sig_out_buf == NULL) {//判空
		printf("sig_out_buf malloc error");//打印错误信息
		return 1;
	}
	memset(sig_out_buf, 0, (size_t)(sig_len + a_len - 1) * sizeof(float));

	for (i = 0; i < sig_len; i++) {
		tmp1 = 0; tmp2 = 0;
		for (j = 0; j < b_len; j++) {
			tmp1 += sig_in_buf[i - j + b_len - 1] * b[j];
		}
		for (j = 1; j < a_len; j++) {
			tmp2 += sig_out_buf[i - j + a_len - 1] * a[j];
		}
		sig_out[i] = (tmp1 - tmp2) / a[0];
		sig_out_buf[i + a_len - 1] = sig_out[i];
	}

	free(sig_in_buf);
	free(sig_out_buf);

	return 0;
}


// 数组逆序 float版
void reverse_arry_float(float* arry, unsigned arry_len) // arry_len可以为1，不报错，此时arry只有一个元素，调用后不变
{
	int i = 0;  //循环变量1, i的值为数组第一个元素的下标
	int j = arry_len - 1;  //循环变量2, j的值为数组最后一个元素的下标
	float fir_idx_buf;  //互换时的中间存储变量
	for (; i < j; ++i, --j)  /*因为i和j已经初始化过了, 所以表达式1可以省略, 但表达式1后面的分号不能省。*/
	{
		fir_idx_buf = arry[i];
		arry[i] = arry[j];
		arry[j] = fir_idx_buf;
	}
}

highpass_iir_05_1.h

#pragma once


extern int iir_filter_zhh(float* sig_in, float* sig_out, int sig_len, float* a, int a_len, float* b, int b_len);
extern void reverse_arry_float(float* arry, unsigned arry_len);

matlab对数据代码：

close all,clear,clc

% fs = 200;
% f1 = 3;f2 = 40;
% t = 0:1/fs:1-1/fs;
% x = sin(2*pi*t*f1)+0.25*sin(2*pi*t*f2);

fid=fopen('D:\zhh\work\VS projects\highpass_iir_05_1\text.txt'); %D:\zhh\work\VS projects\filter_int_v3
x=fscanf(fid,'%f');
fclose(fid);
filename = 'data.txt';
% dlmwrite(filename, x);
fid = fopen(filename, 'w');
for i=1:length(x)
    fprintf(fid, '%.6f\n', x(i));
end
fclose(fid);

% Hd = iir_highpass_05_1_cheby2;
% % 直接matlab滤波
% [b, a] = tf(Hd);
b=[1,2,1,3,1 ];
a=[1,1,1,1,1 ];
y1 = iir_filter_zhh(b,a,x);
y2 = flipud(  filter( b,a,flipud(y1) )  );

figure
plot(x);hold on
plot(y1);
figure
plot(x);hold on
plot(y2);

filename = 'data_fir.txt';
fid = fopen(filename, 'w');
for i=1:length(x)
    fprintf(fid, '%.6f\n', y2(i));
end
fclose(fid);


zhh = 1;