C语言实现matlab的butter函数

//不会用的请参考我的博客:

#include “stdafx.h“
#include 
#include 
#include 




//必需包含下面两个头文件
#include “complex.h“
#include 
double pi = 3.1415926535897932384626433832795;


/*************************************************************************
* @brief   :  Butterworth filter prototype matlab函数
* @inparam :  n		阶数
*			  
* @outparam:  p		复矩阵
*			  k		
* @author  :  zone53
* @date    :  2017/10/18 15:40
* @version :  ver 1.0
*************************************************************************/
static void mybuttap（int n _C_double_complex* p double *k）
{
	int i = 0 j = 0;
	int size = （n - 1） / 2;
	int isodd = n % 2;

	_C_double_complex temp;

	for （i = 1 j = 0; i < n; j++）
	{
		temp._Val[0] = 0;
		temp._Val[1] = pi * i / （2 * n） + pi / 2;
		p[j * 2] = cexp（temp）;
		i += 2;
	}

	for （int m = 1 i = 0; i < j; i++）
	{
		p[m] = conj（p[m - 1]）;
		m += 2;
	}

	if （isodd）
	{
		p[size * 2]._Val[0] = -1.0;
		p[size * 2]._Val[1] = 0.0;
	}

	_C_double_complex a = _Cmulcc（_Cmulcr（p[0] -1） _Cmulcr（p[1] -1））;
	for （int m = 2; m < size * 2 + isodd; m++）
	{
		a = _Cmulcc（a _Cmulcr（p[m] -1））;
	}

	*k = creal（a）;
}

/*************************************************************************
* @brief   :  Characteristic polynomial or polynomial with specified roots matlab函数
* @inparam :  p		复矩阵
*			  np	复矩阵的大小
* @outparam:  d		返回复矩阵
* @author  :  zone53
* @date    :  2017/10/18 15:40
* @version :  ver 1.0
*************************************************************************/
static void mypoly（_C_double_complex* p int np _C_double_complex *d）
{
	_C_double_complex *c = （_C_double_complex *）malloc（（np + 1） * sizeof（_C_double_complex））;
	c[0]._Val[0] = 1.0;
	c[0]._Val[1] = 0.0;
	d[0]._Val[0] = 1.0;
	d[0]._Val[1] = 0.0;
	for （int i = 1; i	{
		c[i]._Val[0] = 0.0;
		c[i]._Val[1] = 0.0;
		d[i]._Val[0] = 0.0;
		d[i]._Val[1] = 0.0;
	}

	_C_double_complex temp;
	for （int i = 0; i < np; i++）
	{
		for （int j = 1; j <= i + 1; j++）
		{
			temp = _Cmulcc（p[i] d[j - 1]）;
			c[j]._Val[0] = d[j]._Val[0] - temp._Val[0];
			c[j]._Val[1] = d[j]._Val[1] - temp._Val[1];
		}
		for （int j = 1; j <= i + 1; j++）
		{
			d[j]._Val[0] = c[j]._Val[0];
			d[j]._Val[1] = c[j]._Val[1];
		}
	}
	free（c）;
}

/*************************************************************************
* @brief   :  实数矩阵相乘
* @inparam :  a		矩阵A
*			  b		矩阵B
*			  m		矩阵A与乘积矩阵C的行数
*			  n		矩阵A的行数矩阵B的列数
*			  k		矩阵B与乘积矩阵C的列数
* @outparam:  c		乘积矩阵 C=AB 
* @author  :  zone53
* @date    :  2017/10/18 15:40
* @version :  ver 1.0
*************************************************************************/
static void mytrmul（double a[] double b[] int m int n int k double c[]）
{
	int i j l u;
	for （i = 0; i <= m - 1; i++）
		for （j = 0; j <= k - 1; j++）
		{
			u = i*k + j; c[u] = 0.0;
			for （l = 0; l <= n - 1; l++）
				c[u] = c[u] + a[i*