rls.c

用c编的一个RLS自适应滤波功能验证程序

#include "stdio.h"
#include "math.h"
#include<stdlib.h>
#define pi 3.1415926
#define delta 0.001
#define Lambda 0.98



double uniform(double a,double b,long int *seed)
{
double t;
*seed=2045*(*seed)+1;
*seed=*seed-(*seed/1048576)*1048576;
t=(*seed)/1048576.0;
t=a+(b-a)*t;
return(t);
}


double gauss(double mean,double sigma,long int *s)
{
int i;
double x,y;
x=0.0;
for(i=0;i<12;i++)
x+=uniform(0.0,1.0,s);
x=x-6.0;;
y=mean+x*sigma;
return(y);
}
void data_save(FILE* fp,double *data,int M)
{
	int i;
	if((fp=fopen("RLS","wb"))==NULL)
	{
		printf("cannot open file\n");
		return;
	}
	for(i=0;i<M;i++)
		if(fwrite(&data[i],sizeof(double),1,fp)!=1)
			printf("file write error\n");
	fclose(fp);
}
void main()	
{

int i,j;
       //时域抽头RLS算法滤波器阶数
double *d;
double *w;
double *x;
double *error;
double *g;
double *p;
double *u;
double *z;
double *y;
long int s=13579;
FILE *fp;

int N=2048;   //输入信号抽样点数
int k=64; 

d=(double *)malloc(N*sizeof(double));
error=(double *)malloc(N*sizeof(double));
x=(double *)malloc(N*sizeof(double));
g=(double *)malloc(k*sizeof(double));
w=(double *)malloc(k*sizeof(double));
p=(double *)malloc(k*k*sizeof(double));
y=(double *)malloc(N*sizeof(double));
u=(double *)malloc(k*sizeof(double));
z=(double *)malloc(k*sizeof(double));
for(i=0;i<N;i++)
{
d[i]=sin(0.05*pi*i);
x[i]=d[i]+gauss(0.0,0.3,&s);
error[i]=0;
if(i>=k)y[i]=0;
}
for(i=0;i<k;i++)
{
for(j=0;j<k;j++)
{
	if(i==j)
p[i*k+j]=delta;
	else
		p[i*k+j]=0;
}
}
for(i=0;i<k;i++)
{
	y[i]=x[i];
    w[i]=0;
    g[i]=0;
	z[i]=0;
	u[i]=0;
}           for(i=k-1;i<=N-1;i++)
			{

int c;
double f;
double tt;

c=0;
for(j=i;j>=i-k+1;j--) 
{
	u[c]=x[j];
	c=c+1;
}

tt=0;
for(c=0;c<k;c++)
{
	
tt=tt+w[c]*u[c];
}
y[i]=tt;
    for(c=0;c<k;c=c+1)
	{
		z[c]=0;
		for(j=0;j<k;j=j+1)
		{
			z[c]=z[c]+u[j]*p[k*j+c];
		}
	}
	f=0;
    for(c=0;c<k;c++)
	{
		
    f=f+z[c]*u[c];
	}
	f=f+Lambda;
   for(c=0;c<k;c++)
   {
      g[c]=z[c]/f;
   }
error[i]=d[i]-y[i];
printf("%f\n",error[i]);
for(c=0;c<k;c++)
{
w[c]=w[c]+g[c]*error[i];
}

for(j=0;j<k;j++)
{
	for(c=0;c<k;c++)
	{

    p[j*k+c]=(p[j*k+c]-g[j]*z[c])/Lambda;
	}
}

			}

data_save(fp,error,N);			
free(d);
free(w);
free(x);
free(error);
free(g);
free(p);
free(u);
free(z);
free(y);	
}