matlab butter函数 C实现

参考文档

主要参考资料来源数字信号处理 基于计算机的方法 4版（中文版）和Matlab函数butter(巴特沃斯IIR滤波器) 两部分。

IIR滤波器和分类

• IIR滤波器
  即无限冲激响应的滤波器
  缺点 是 相位响应非线性（相位畸变）且在高阶时不稳定（过度或发散响应）
  优点 是相较于FIR而言高效、实时、频率选择等。

• 分类
  如今主要的IIR滤波器有以下几类：

butterworth(巴特沃斯)

chebyshev1(切比雪夫Ⅰ型)

chebyshev2(切比雪夫Ⅱ型)

elliptic(椭圆)

bessel(贝塞尔)

（这篇文章主要介绍巴特沃斯IIR滤波器的C实现）

C代码

void butter(int n, double Wn[], int type, int analog, double *ab, double *bb)
{
   
   
	double fs = 2;
	double u[2] = {
   
   0.0, 0.0};
	// step 1: get analog, pre-warped frequencies
	if (!analog)
	{
   
   
		if (type == 1 || type == 2)
		{
   
   
			fs = 2;
			u[0] = 2 * fs * tan(PI * Wn[0] / fs);
		}
		else
		{
   
   
			fs = 2;
			u[0] = 2 * fs * tan(PI * Wn[0] / fs);
			u[1] = 2 * fs * tan(PI * Wn[1] / fs);
		}
	}
	else if (type == 3 || type == 4)
	{
   
   
		if (type == 1 || type == 2)
		{
   
   
			u[0] = Wn[0];
		}
		else
		{
   
   
			u[1] = Wn[1];
		}
	}

	// step 2: convert to low-pass prototype estimate
	double Bw = 0.0;
	if (type == 1 || type == 2)
	{
   
   
		Wn = u;
	}
	else if (type == 3 || type == 4)
	{
   
   
		Bw = u[1] - u[0];
		Wn[0] = sqrt(u[0] * u[1]); // center
		Wn[1] = 0.0;
	}

	// step 3: Get N-th order Butterworth analog lowpass prototype
	Complex *z = (Complex *)malloc(0 * sizeof(Complex));
	Complex *p = (Complex *)malloc(n * sizeof(Complex));
	double k = 0;
	buttap(n, p, &k);
	
	// Transform to state-space
	int a_size = n;
	double *a = (double *)malloc(sizeof(double) * n * n);
	double *b = (double *)malloc(sizeof(double) * n);
	double *c = (double *)malloc(sizeof