resample函数的C语言实现

由于matlab程序将各信号处理函数封装，无法看到原始代码，这里是其中resample函数的C语言实现源码，可实现任意整数倍上采及下采

/*
FIR滤波器的I 和 II两种类型

调用格式：Firls_LP（double * h_firls int Lendouble Fint M）;
len表示滤波器的阶数;
M = [1 1 0 0];
F给表示频率范围，格式为4*1的数组;
输出一个len+1的数组，由h返回;
*/
# include “resample_definition.h“
# include 
# include 
# include 
# include 


//# define PI 3.1415926
//# define EPSILON 0.000001

//double * Firls_LP（int double * int *）;   //计算FIR滤波器抽头系数
//bool Constant_Diff（double *）;             //判断频率范围是否合理
//double Sinc（double）;                      //辛格函数
//double * Divide_Matrix（double ** double * int）;  //矩阵左除
//void Inverse_Matrix（double ** int）;      //求逆矩阵

//int main（void）
//{
//	double * h;                         //h数组存放firls输出
//	int Len;                            //根据调用函数设置，测试时给定，可变
//	int M[4] = {1 1 0 0};            //低通滤波器幅度范围，恒定
//	double F[4] = {0 0.2 0.4 1};	    //根据调用函数设置，测试时给定，可变
//	printf（“请输入滤波器阶数：“）;
//	scanf（“%d“ &Len）;	
//	h = （double *）malloc（sizeof（double） * （Len + 1））;
//	h = Firls_LP（Len F M）;	           //len代表滤波器阶数，F为频率范围，M为幅度范围          
//
//	return 0;
//}

void Firls_LP（double *h_firlsint N double * F int * M）
{
	int i j fullband constant_weight Nodd L need_matrix s;
	constant_weight = 1;               //类型I和II滤波器，constant_weight值恒定，由w矩阵计算而来
	int W[2] = {1 1};                //类型I和II滤波器，w矩阵值恒定

	double * m * k_Odd *k_Even * b_Even * b_Odd * bNew * a;
	double ** I1 ** I2 **G;         
	double b0 m0 b1;	
	N = N + 1;                        //N阶滤波器具有N+1个系数
	Nodd = N%2;                       //判断序列长度奇偶性
	L = （N-1）/2;                      //中心对称点
	m = （double *）malloc（sizeof（double） * （L+1））;


	for （i=0; i<4; i++）
		F[i] /= 2;		              //频率归一化

	if （Constant_Diff（F））             //判断频率是否恒定
		fullband = 1;
	else
		fullband = 0;		                      

	if （!Nodd）
	{
		//序列长为奇数
		for （i=0; i<=L; i++）		
			m[i] = i + 0.5;		
	}
	else   //序列长为偶数
	{
		for （i=0; i<=L; i++）		
			m[i] = i;		
	}	
	need_matrix = （!fullband）||（!constant_weight）;    //参数

	if （need_matrix）  
	{
		//定义（L+1）*（L+1）大小的方阵，其值等于m方阵+m方阵的转置矩阵
		I1 = （double **）malloc（sizeof（double*） * （L+1））;
		for （i=0; i<=L; i++）
			*（I1 + i） = （double *）malloc（sizeof（double） * （L+1））;	

		for（i=0; i<=L; i++）
			for （j=0; j<=L; j++）
				I1[i][j] = m[i] + m[j];

		//定义（L+1）*（L+1）大小的方阵，其值等于m方阵-m方阵的转置矩阵
		I2 = （double **）malloc（sizeof（double*） * （L+1））;
		for （i=0; i<=L; i++）
			*（I2 + i） = （double *）malloc（sizeof（double） * （L+1））;

		for（i=0; i<=L; i++）
			for （j=0; j<=L; j++）
				I2[i][j] = m[i] - m[j];

		//定义（L+1）*（L+1）大小的零矩阵
		G = （double **）malloc（sizeof（double*） * （L+1））;
		for （i=0; i<=L; i++）				
			*（G + i） = （double *）malloc（sizeof（double） * （L+1））;			

		for（i=0; i<=L; i++）
			for （j=0; j<=L; j++）
				G[i][j] = 0;
	}
	k_Even = （double *）malloc（sizeof（double） * L）;       //偶阶次的k序列
	k_Odd = （double *）malloc（sizeof（double） * （L+1））;    //奇阶次序列的k序列
	bNew = （double *）malloc（sizeof（double） * （L+1））;

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     目录           0  2018-10-23 23:54  resample_1029\
     目录           0  2018-10-23 23:49  resample_1029\Debug\
     文件       47104  2018-10-23 23:49  resample_1029\Debug\resample_1029.exe
     文件      357500  2018-10-23 23:49  resample_1029\Debug\resample_1029.ilk
     文件      543744  2018-10-23 23:49  resample_1029\Debug\resample_1029.pdb
     目录           0  2018-10-23 23:50  resample_1029\resample_1029\
     文件    24268800  2018-10-23 23:54  resample_1029\resample_1029.sdf
     文件         906  2018-10-23 23:49  resample_1029\resample_1029.sln
     文件     1038774  2018-10-23 23:50  resample_1029\resample_1029\data_i.txt
     文件     1039121  2018-10-23 23:50  resample_1029\resample_1029\data_r.txt
     目录           0  2018-10-23 23:50  resample_1029\resample_1029\Debug\
     文件        2638  2018-10-23 23:50  resample_1029\resample_1029\Debug\cl.command.1.tlog
     文件       18648  2018-10-23 23:50  resample_1029\resample_1029\Debug\CL.read.1.tlog
     文件        2382  2018-10-23 23:50  resample_1029\resample_1029\Debug\CL.write.1.tlog
     文件       20694  2018-10-23 23:50  resample_1029\resample_1029\Debug\firls.obj
     文件       12095  2018-10-23 23:50  resample_1029\resample_1029\Debug\kaiser1.obj
     文件           2  2018-10-23 23:50  resample_1029\resample_1029\Debug\link-cvtres.read.1.tlog
     文件           2  2018-10-23 23:50  resample_1029\resample_1029\Debug\link-cvtres.write.1.tlog
     文件           2  2018-10-23 23:50  resample_1029\resample_1029\Debug\link.4300-cvtres.read.1.tlog
     文件           2  2018-10-23 23:50  resample_1029\resample_1029\Debug\link.4300-cvtres.write.1.tlog
     文件           2  2018-10-23 23:50  resample_1029\resample_1029\Debug\link.4300.read.1.tlog
     文件           2  2018-10-23 23:50  resample_1029\resample_1029\Debug\link.4300.write.1.tlog
     文件           2  2018-10-23 23:50  resample_1029\resample_1029\Debug\link.6936-cvtres.read.1.tlog
     文件           2  2018-10-23 23:50  resample_1029\resample_1029\Debug\link.6936-cvtres.write.1.tlog
     文件           2  2018-10-23 23:50  resample_1029\resample_1029\Debug\link.6936.read.1.tlog
     文件           2  2018-10-23 23:50  resample_1029\resample_1029\Debug\link.6936.write.1.tlog
     文件        1986  2018-10-23 23:50  resample_1029\resample_1029\Debug\link.command.1.tlog
     文件        3406  2018-10-23 23:50  resample_1029\resample_1029\Debug\link.read.1.tlog
     文件        1132  2018-10-23 23:50  resample_1029\resample_1029\Debug\link.write.1.tlog
     文件         420  2018-10-23 23:50  resample_1029\resample_1029\Debug\mt.command.1.tlog
     文件         410  2018-10-23 23:50  resample_1029\resample_1029\Debug\mt.read.1.tlog
............此处省略22个文件信息