二阶IIR滤波器C代码实现

#define PI 3.14159265358979323846
typedef double sample_t;

enum {
	BIQUAD_LOWPASS,
	BIQUAD_HIGHPASS,
	BIQUAD_BANDPASS_PEAK,
};

typedef struct
{
	sample_t a0, a1, a2, a3, a4;
	sample_t x1, x2, y1, y2;
}biquad_state;

// 计算滤波器的系数，2nd-order Butterworth: q_value=0.7071 ，2nd-order Chebyshev (ripple 1 dB): q_value=0.9565，2nd-order Thomson-Bessel: 
// q_value=0.5773。带通滤波器的q_value = sqrt(f1*f2)/(f2-f1)，fs是指信号的采样频率，fc为cutoff frequence 或者是centern frequenc(带通滤波器)
void filter_init(biquad_state *state, int type, int fs, float fc, float q_value)
{
	sample_t w0, sin_w0, cos_w0, alpha;
	sample_t b0, b1, b2, a0, a1, a2;

	w0 = 2 * PI * fc / fs;
	sin_w0 = sin(w0);
	cos_w0 = cos(w0);
	alpha = sin_w0 / (2.0 * q_value);

	switch(type)
	{
	case BIQUAD_LOWPASS:
		b0 = (1.0 - cos_w0) / 2.0;
		b1 = b0 * 2;
		b2 = b0;
		a0 = 1.0 + alpha;
		a1 = -2.0 * cos_w0;
		a2 = 1.0 - alpha;
		break;
	case BIQUAD_HIGHPASS:
		b0 = (1.0 + cos_w0) / 2.0;
		b1 = -b0 * 2;
		b2 = b0;
		a0 = 1.0 + alpha;
		a1 = -2.0 * cos_w0;
		a2 = 1.0 - alpha;
		break;
	case BIQUAD_BANDPASS_PEAK:
		b0 = alpha;
		b1 = 0.0;
		b2 = -alpha;
		a0 = 1.0 + alpha;
		a1 = -2.0 * cos_w0;
		a2 = 1.0 - alpha;
		break;
	}
	state->a0 = b0 / a0;
	state->a1 = b1 / a0;
	state->a2 = b2 / a0;
	state->a3 = a1 / a0;
	state->a4 = a2 / a0;
	state->x1 = state->x2 = 0.0;
	state->y1 = state->y2 = 0.0;
}

//biquard滤波
sample_t biquad(biquad_state *state, sample_t data)
{
	sample_t result = 0;
	result = state->a0 * data + state->a1 * state->x1 + state->a2 * state->x2 -  state->a3 * state->y1 - state->a4 * state->y2;
	state->x2 = state->x1;
	state->x1 = data;
	state->y2 = state->y1;
	state->y1 = result;
	return result;

}

int main(int argc, char **argv) {
    biquad_state *stat = (biquad_state *) malloc(sizeof(biquad_state));
    int numbers = 1024;
    float fs = 1024;
    sample_t data[numbers], filter_data[numbers];
    float fc = 50;
    float q_value = 0.7071;
    int type = BIQUAD_LOWPASS;
    filter_init(stat, type, fs, fc, q_value);
    for (int i = 0; i <numbers; ++i)
    {　
        data[i] = 0.5 * (sin(2 * PI * 40 * i / fs) + cos(2 * PI * 150 * i / fs + PI / 4));
    }
　　　　
    for (int i = 0; i <numbers; ++i)
    {
        filter_data[i] = biquad(stat, data[i]);
    }
    free(stat);
}