24p.c

24点傅立叶算法

#include <stdio.h>
#include <conio.h>
#include <math.h>

#define PI	3.1416
#define DU	(PI/180)

int vsin,vcos,kkff;

void fourier_table(int points)
{
	int i,k;
	float f;
	float ks1[100],kc1[100];
	FILE *file;
	f=PI*2.0/(float)points;
	k=points;
	for (i=0;i<=k;i++)
	{
		ks1[i]=sin(i*f);
		kc1[i]=cos(i*f);
	}
	file=fopen("fourier.txt","wa");
	fprintf(file,"%d points\n",points);
	fprintf(file,"sin:\t\t\t\tcos:\n");
	for(i=0;i<=k;i++)
	{
		fprintf(file,"#define ks%d\t%f\t\t\t#define kc%d\t%f\n",i,ks1[i],i,kc1[i]);
	}
}
void kf(int a,int b)
{
	int c,d;
	if(a<0) a=~a+1;
	if(b<0) b=~b+1;
	if(a<b)
	{
	 c=a;
	 a=b;
	 b=c;
	}
	c=a+a+a+b+1;
	d=b+b+b+b+b;
 asm  {
	mov ax,d;
//  cwd
	imul b;
	idiv c;
	mov bx,1000;
	imul bx;
	mov bx,3000;
	idiv bx;
	add a,ax
//  mov ax,a;
			}
 kkff=a;
}

void f24(int * ptr,int i)//24point
{
//sin :
// k1 ( I0 + I10 - I12 - I22 )
// k2 ( I1 + I9  - I13 - I21 )
// k3 ( I2 + I8  - I14 - I20 )
// k4 ( I3 + I7  - I15 - I19 )
// k5 ( I4 + I6  - I16 - I18 )
// k6 ( I5       - I17       )
//cos :
// k1 ( I4 - I6  - I16 + I18 )
// k2 ( I3 - I7  - I15 + I19 )
// k3 ( I2 - I8  - I14 + I20 )
// k4 ( I1 - I9  - I13 + I21 )
// k5 ( I0 - I10 - I12 + I22 )
// k6 (    - I11       + I23 )

// the input is address of be calculated date
// the output date is register of:
//	AX	the value if cos
//	CX	the value of sin
	_SI=(unsigned int)ptr;
	asm{
	mov	ax,i;
	shl	ax,1;
	add si,ax;

// k1 ( I0 + I10 - I12 - I22 )
	mov	ax,[si];
	add	ax,[si-20];
	sub	ax,[si-24];
	sub	ax,[si-44];
	//k1=2588,k2=5000,k3=7071,k4=8660;k5=9659;k6=10000
	mov cx,2588
	sal cx,1
	imul cx
	mov bx,ax
	mov di,dx

// k2 ( I1 + I9  - I13 - I21 )
	mov	ax,[si-2];
	add	ax,[si-18];
	sub	ax,[si-26];
	sub	ax,[si-42];
	//k1=2588,k2=5000,k3=8090,k4=8660;k5=9659
	mov cx,5000
	sal cx,1
	imul cx
	clc
	add bx,ax
	adc di,dx

// k3 ( I2 + I8  - I14 - I20 )
	mov	ax,[si-4];
	add	ax,[si-16];
	sub	ax,[si-28];
	sub	ax,[si-40];
	//k1=2588,k2=5000,k3=7071,k4=8660;k5=9659
	mov cx,7071
	sal cx,1
	imul cx
	clc
	add bx,ax
	adc di,dx

// k4 ( I3 + I7  - I15 - I19 )
	mov	ax,[si-6];
	add	ax,[si-14];
	sub	ax,[si-30];
	sub	ax,[si-38];
	//k1=2588,k2=5000,k3=7071,k4=8660;k5=9659
	mov cx,8660
	sal cx,1
	imul cx
	clc
	add bx,ax
	adc di,dx

// k5 ( I4 + I6  - I16 - I18 )
	mov	ax,[si-8];
	add	ax,[si-12];
	sub	ax,[si-32];
	sub	ax,[si-36];
	//k1=2588,k2=5000,k3=7071,k4=8660;k5=9659
	mov cx,9659
	sal cx,1
	imul cx
	clc
	add bx,ax
	adc di,dx

// k6 ( I5       - I17       )
	mov	ax,[si-10];
	sub	ax,[si-34];
	//k1=2588,k2=5000,k3=7071,k4=8660;k5=9659
	mov cx,10000
	sal cx,1
	imul cx
	clc
	add ax,bx
	adc dx,di

	}
	vsin=_DX;
// end of sin calculation.

// Begin cos calculation:

//cos :
	asm{

// k1 ( I4 - I6  - I16 + I18 )
	mov	ax,[si-8];
	sub	ax,[si-12];
	sub	ax,[si-32];
	add	ax,[si-36];
	//k1=2588,k2=5000,k3=7071,k4=8660;k5=9659
	mov cx,2588
	sal cx,1
	imul cx
	mov bx,ax
	mov di,dx

// k2 ( I3 - I7  - I15 + I19 )
	mov	ax,[si-6];
	sub	ax,[si-14];
	sub	ax,[si-30];
	add	ax,[si-38];
	//k1=2588,k2=5000,k3=7071,k4=8660;k5=9659
	mov cx,5000
	sal cx,1
	imul cx
	clc
	add bx,ax
	adc di,dx

// k3 ( I2 - I8  - I14 + I20 )
	mov	ax,[si-4];
	sub	ax,[si-16];
	sub	ax,[si-28];
	add	ax,[si-40];
	//k1=2588,k2=5000,k3=7071,k4=8660;k5=9659
	mov cx,7071
	sal cx,1
	imul cx
	clc
	add bx,ax
	adc di,dx

// k4 ( I1 - I9  - I13 + I21 )
	mov	ax,[si-2];
	sub	ax,[si-18];
	sub	ax,[si-26];
	add	ax,[si-42];
	//k1=2588,k2=5000,k3=7071,k4=8660;k5=9659
	mov cx,8660
	sal cx,1
	imul cx
	clc
	add bx,ax
	adc di,dx

// k5 ( I0 - I10 - I12 + I22 )
	mov	ax,[si];
	sub	ax,[si-20];
	sub	ax,[si-24];
	add	ax,[si-44];
	//k1=2588,k2=5000,k3=7071,k4=8660;k5=9659
	mov cx,9659
	sal cx,1
	imul cx
	clc
	add bx,ax
	adc di,dx

// k6 (    - I11       + I23 )
	mov	ax,[si-46];
	sub	ax,[si-22];
	//k1=2588,k2=5000,k3=7071,k4=8660;k5=9659
	mov cx,10000
	sal cx,1
	imul cx
	clc
	add ax,bx
	adc dx,di

	}
	vcos=_DX;
}

void main(void) {

	int i,j;
	int x[1000];
	int *p;
	long l,ls,lc;
	
	fourier_table(24);
	getch();

	for(i=0;i<50;i++) {
		x[i]=(int)( sin(15*DU*i)*2200 );
	}

	for(i=0;i<30;i++)
		printf("x[%d]= %d\t",i,x[i]);


	for(j=0;j<10;j++) {

		p=&( x[j] );
		f24(p,24);	          //2020

		printf("\n vsin=%d , vcos=%d\t\t",vsin,vcos);

		kf(vsin,vcos);

		printf("%d-> %d\t",j,kkff);

	}
	puts("\n\n");
	getch();

}