//不會用的請參考我的博客:

#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?	{
		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?	{
		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?	{
		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?	{
		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*