fourier_trans.c

傅立叶变换的C语言实现

/* 
************************************************************************
*  fourier.c: Discrete Fourier Transformation                          *
*								       *
*  taken from: "Projects in Computational Physics" by Landau and Paez  * 
*	       copyrighted by John Wiley and Sons, New York            *      
*                                                                      *
*  written by: students in PH465/565, Computational Physics,           *
*	       at Oregon State University                              *
*              code copyrighted by RH Landau                           *
*  supported by: US National Science Foundation, Northwest Alliance    *
*                for Computational Science and Engineering (NACSE),    *
*                US Department of Energy 	                       *
*								       *
*  UNIX (DEC OSF, IBM AIX): cc fourier.c -lm    		       *
*  								       *
*  comment: The program reads its input data from a file in the        *
* 	    same directory called input.dat. This file has to contain  *
*           only y(t) values separated by whitespaces which are real.  *
*	    The output is the direct output from the algorithm which   *
*	    will probably look very different than what you are used   *
*           to.	The output has the form 			       *
*	    	frequency index \t real part \t imaginary part         *
*  related programs: invfour.c					       *
************************************************************************
*/
#include <stdio.h>
#include <math.h>

#define max 1000			/* max number of input data */
#define PI 3.1415926535897932385E0

main()
{
   double imag, real,input[max+1];
   int i=0,j,k;
 
   FILE *data;
   FILE *output;				
   data=fopen("input.dat", "r");	/* read data from input.dat */
   output=fopen("fourier.dat", "w");	/* save data in fourier.dat */

   while ((fscanf(data, "%lf", &input[i]) !=EOF) && (i<max))
   { 				
      i++;				/* reading input data */
   }
 
   for (j=0; j<i; j++)			/* loop for frequency index */
   {
      real=imag=0.0;			/* clear variables */
      for (k=0; k<i; k++)		/* loop for sums */
      {
         real+=input[k]*cos((2*PI*k*j)/i);
         imag+=input[k]*sin((2*PI*k*j)/i);
      }
      fprintf(output, "%d\t%f\t%f\n", j, real/i, imag/i );
   }
   printf("data stored in fourier.dat.\n");
   fclose(data);
   fclose(output);
}