fourier.c

来自「利用不均一分割的FDTD解析程序」· C语言 代码 · 共 56 行

/*
 *		Fourier Transform
 */

void fourier(void)
{
	double vz[nt+1];
	double cur[nt+1];
	long_double_complex an1[501],an2[501];

	long_double_complex Z;
	int a,n,i;
	double freq,loss,F,omega,current,voltage;

	for(i=0;i<=500;i++){
		an1[i] = (0.0,0.0);
		an2[i] = (0.0,0.0);
	}

	/*****    Open input file "v-i.dat"     *****/
	if ((fin=fopen("v-i.dat","r"))==NULL){
		printf("file %s can't open\n","v-i.dat");
		exit(-1);
	}

	/*****    Open output file "z-loss.dat"     *****/
	if((fout=fopen("z-loss.dat","wb"))==NULL){
		printf("file %s can't open\n","z-loss.dat");
		exit(-1);
	}

	/*****    Read data from input file    *****/
	while((fscanf(fin,"%d\t%lf\t%lf\n",&a,&voltage,&current))!=EOF){
		vz[a] = voltage;
		cur[a] = current;
	}
	fclose(fin);

	for(i=0;i<=500;i++){
		freq = (1000.0e+6)/500*i;
		omega = 2.0*PI*freq;
		for(n=0;n<=nt;n++){
			an2[i] += vz[n]*exp(-J*omega*n*ts);
			an1[i] += cur[n]*exp(-J*omega*n*ts);
		}
	}

	for(i=0;i<=500;i++){
		Z = an2[i]/an1[i]*2.0;                   // Input Impedance
		F = (double)(i)/500;
		loss = (double)20.0*log10(abs((Z-50)/(Z+50))); // Return Loss
		fprintf(fout,"%lf\t%lf\t%lf\t%lf\n",F,Z.real(),Z.imag(),loss);
	}
	fclose(fout);
}