DSP滤波程序C语言实现

#include <stdio.h>
#include <malloc.h>
//z变换转成差分方程然后迭代
//Talen @UESTC 2015.6.5

//打印数组
void print(float *x,int len);
//取后第1-n个
float* aminusone(float* a,int ar);
//差分中间值
float gtcet(float *t,float *x,int xr,int s,int d);
//主要滤波函数
float* gfilter(float *b,int br,float *a,int ar,float *x,int xr);

void print(float *x,int len)
{
for(int i=0;i<len;i++)
{printf("%9.3f ",x[i]);
if((i%7)==6)
printf("\n");}
printf("\n");
}

float* aminusone(float* a,int ar){
float* am = (float *)malloc((ar-1) * sizeof(float));
for(int i=0;i<ar-1;i++)
am[i]=a[i+1];
return am;
}
float gtcet(float *t,float *x,int xr,int s,int d)
//x[]数据源 r数据源长度 s开始点(0开始) d数据长度
//t可以是a[] b[] 长度与d一样
{
float sum=0;
float* w = (float *)malloc(d * sizeof(float));
for(int i=0;i<d;i++)
{
int temp=s-i;
if(temp>=0&&temp<xr)
w[i]=x[temp];
else
w[i]=0;
//核心语句
sum=sum+w[i]*t[i];
}
return sum;
}
float* gfilter(float *b,int br,float *a,int ar,float *x,int xr)
//b a x br b的长度 ar a的长度 xr数据长度
//use本程序定义的gtcet()和aminusone()
//输入数据，b，a，x得出滤波后的y
{
float* y = (float *)malloc(xr* sizeof(float));
for(int ti=0;ti<xr;ti++)
y[ti]=0;
for(int i=0;i<xr;i++)
{
//最核心语句，此处迭代
y[i]=gtcet(b,x,xr,i,br)-gtcet(aminusone(a,ar),y,xr,i-1,ar-1);
}
return y;
}

int main()
{
//input 源数据 可修改
////////////////////////////////////////////////////////////////
float x[120]={-4, -2, 0, -4, -6, -4, -2, -4, -6, -6, -4, -4,\
-6, -6, -2, 6, 12, 8, 0, -16, -38, -60, -84,\
-90, -66, -32, -4, 2, -4, 8, 12, 12, 10, 6,\
6, 6, 4, 0, 0, 0, 0, 0, -2, -4, 0, 0, 0, -2,\
-2, 0, 0, -2, -2, -2, -2, 0, -2, -4, -2, 0,\
-2, -4, -4, 2, 0, 0, -2, -4, -2, 0, 0, -2,\
-4, -2, 0, 0, -4, -4, -2, -2, -4, -6, -6,\
-4, -4, 8, -10, -8, -6, -6, -8, -12, -10,\
-8, -8, -10, -12, -10, -8, -8, -10, -10,\
-8, -6, -6, -8, -8, -4, -2, -4, -4, -4,\
0, 0, -2, -4, -2, -2, 0, -4};
//参数
float b[7]={ 0.0007,0.0044,0.0111,0.0148,0.0111,0.0044,0.0007};
float a[7]={ 1.0000,-3.1836,4.6222,-3.7795,1.8136,-0.4800,0.0544};
//////////////////////////////////////////////////////////////////
//不要修改
int xr=sizeof(x)/sizeof(float);
int br=sizeof(b)/sizeof(float);
int ar=sizeof(a)/sizeof(float);
//output 滤波后数据
float *y;

printf("参数b：\n");print(b,br);
printf("参数a：\n");print(a,ar);
printf("原数据：\n");print(x,xr);

y=gfilter(b,br,a,ar,x,xr);

printf("滤波后数据：\n");print(y,xr);
getchar();
return 0;
}