echo_canceller.c

自适应滤波中128ms回声消除VC6.0程序

#include <stdio.h>
#include <math.h>
#define M 128
#define mu 0.005
//#define mu 0.5
#define gama 0.20
#define S16 signed short
#define ERROR 10
#define NLMS

int main()
{
	FILE * PRIMARY;
	FILE * REFERENCE;
	FILE * OUTPUT1;
       FILE * OUTPUT2;
	FILE * OUTPUT3;
	FILE * OUTPUT4;
	int i;	
	S16 u[M] = {0};
	S16 d;
	double w[M] = {0};
	double y, e;
	int n = 0;
	long double m = 0;
	long double p = 0;
	double E[M]={0};
	double abs = 0;
	double power = 0;

    PRIMARY = fopen("D:\\course\\asp\\adaptive homework 2\\d.dat", "rb");
	if(PRIMARY == NULL)
	{
		printf("Error: Unable to open d.dat!\n");
		return 1;
	}

	REFERENCE = fopen("D:\\course\\asp\\adaptive homework 2\\u.dat", "rb");
    if(REFERENCE == NULL)
	{
		printf("Error: Unable to open u.dat!\n");
		return 1;
	}

	OUTPUT1 = fopen("D:\\course\\asp\\adaptive homework 2\\output1.txt", "w");
    if(REFERENCE == NULL)
	{
		printf("Error: Unable to open output1.txt!\n");
		return 1;
	}

	OUTPUT2 = fopen("D:\\course\\asp\\adaptive homework 2\\output2.txt", "w");
    if(REFERENCE == NULL)
	{
		printf("Error: Unable to open output2.txt!\n");
		return 1;
	}

	OUTPUT3 = fopen("D:\\course\\asp\\adaptive homework 2\\output3.txt", "w");
    if(REFERENCE == NULL)
	{
		printf("Error: Unable to open output3.txt!\n");
		return 1;
	}

	OUTPUT4 = fopen("D:\\course\\asp\\adaptive homework 2\\output4.txt", "w");
    if(REFERENCE == NULL)
	{
		printf("Error: Unable to open output4.txt!\n");
		return 1;
	}	

	fprintf(OUTPUT1, "e=[");
	fprintf(OUTPUT2, "d=[");
	fprintf(OUTPUT3, "m=[");
	fprintf(OUTPUT4, "s=[");
    while(fread(&d, sizeof(S16), 1, PRIMARY) && fread(&u[M-1], sizeof(S16), 1, REFERENCE))
	{		
		y = 0;
		for(i=0;i<M;i++)
		{
			y += w[i] * u[i];
		}		
		e = (double)d - y;		
        #ifdef NLMS		  
		  abs = fabs(u[M-1]) / pow(2, 15);
		  power = (1 - gama) * power + gama * abs;
		  for(i=0;i<M;i++)
		  {
			  w[i] += (mu/pow(2, 30)) * u[i] * e / (power * power);
		  }
        #else
		  for(i=0;i<M;i++)
		  {
			  w[i] += (mu/pow(2, 30)) * u[i] * e;
		  }
        #endif
		m += e * e;
		n++;
		for(i=0;i<M-1;i++)
		{
			u[i] = u[i+1];
		}
		for(i=0;i<M-1;i++)
		{
			E[i] = E[i+1];
		}
		E[M-1] = e;
		p = 0;
		for(i=M-1;i>=0;i--)
		{
			p += E[i] * E[i];
		}
		p = p/(32768*32768);
		if(p<=0)
			p = 1e-2;
		fprintf(OUTPUT1, "%e ", e/32768);
		fprintf(OUTPUT2, "%e ", d/32768);
		fprintf(OUTPUT3, "%e ", e*e/(32768*32768));
		if(n == M)
			for(i=0;i<M;i++)
				fprintf(OUTPUT4, "%e ", 3.14 + 20 * log(sqrt(2*p/M)/4096));
		else if(n>M)
			fprintf(OUTPUT4, "%e ", 3.14 + 20 * log(sqrt(2*p/M)/4096));		
	}
	fprintf(OUTPUT1, "];\n");
	fprintf(OUTPUT2, "];\n");
	fprintf(OUTPUT3, "];\n");
	fprintf(OUTPUT4, "];\n");
	
	fclose(PRIMARY);
	fclose(REFERENCE);
	fclose(OUTPUT1);
       fclose(OUTPUT2);
	fclose(OUTPUT3);
	fclose(OUTPUT4);

	return 0;
}