C版本- #include <stdlib.h>
- #include <stdio.h>
- #include <string.h>
- #include <math.h>
- #define MAXLINESIZE 128
- /*Calculates the Fourier transform over Nf points from frequency f1 to f2,of N(N>1) samples of
- a signal taken at sps sample per second.The values of cos(n*Omega)&sin(n*Omega) are calculated
- using Chebyshev polynomials.*/
- void hrft(int N,double *q,double Omega,double *Qr,double *Qi)
- {
- int n;
- double C0,S0;
- double T0,T1,T2;
- double U0,U1,U2;
- double Cn,Sn;
- C0=cos(Omega);
- S0=sin(Omega);
- T0=1.0;
- T1=C0;
- U0=1.0;
- U1=2.0*C0;
- *Qr=q[0]+q[1]*T1;
- *Qi=q[1]*S0*U0;
- C0 *=2.0;
- for(n=2;n<N;++n)
- {
- T2=C0*T1-T0;
- T0=T1;
- T1=T2;
- U2=C0*U1-U0;
- U0=U1;
- U1=U2;
- Cn=T2;
- Sn=S0*U0;
- *Qr +=q[n]*Cn;
- *Qi +=q[n]*Sn;
- }
- }
- int main(int argc,char *argv[])
- {
- double f1,f2;
- int Nf;
- double df;
- int N;
- int i,n;
- double sps;
- double Omega;
- double dOmega;
- FILE *fp;
- char linebuff[MAXLINESIZE];
- double *q;
- double Qr,Qi;
- if (argc !=8)
- {
- printf ("\nhrft calculates the Fourier transform at Nf equally\n");
- printf ("spaced points in the frequency interval [f1,f2]\n");
- printf ("usage: hrft f1 f2 Nf N sps infile outfile\n");
- printf (" f1=starting frequency [Hz]\n");
- printf (" f2=ending frequency [Hz]\n");
- printf (" Nf=number of frequengcy points to calculate FT\n");
- printf (" N=number of samples to read from the input file\n");
- printf (" sps=samples per second\n");
- printf (" infile=input file\n");
- printf ("outfile=output file\n");
- exit (-1);
- }
- f1=atof(argv[1]);
- f2=atof(argv[2]);
- Nf=atoi(argv[3]);
- N=atoi(argv[4]);
- sps=atof(argv[5]);
- if (N<2)
- {
- printf("N must be greater than 1\n");
- exit (-1);
- }
- if ((fp=fopen(argv[6],"r"))==NULL)
- {
- perror("Error opening input file");
- exit (-1);
- }
- printf("Reading data file: %s\n\n",argv[6]);
- for(n=0;n<N;++n)
- {
- fgets (linebuff,MAXLINESIZE,fp);
- if (linebuff[0] != '#') break;
- printf("%s",linebuff);
- }
- q=(double *)malloc(N*sizeof(double));
- sscanf (linebuff,"%lf",&q[0]);
- for (n=0;n<N;++n)
- {
- fscanf(fp,"%lf",&q[n]);
- }
- fclose(fp);
- if ((fp=fopen(argv[7],"w"))==NULL)
- {
- perror("Error opening output file");
- exit (-1);
- }
- printf ("Opening output file:%s\n",argv[7]);
- df=(f2-f1)/((double)(Nf-1));
- fprintf(fp,"# This file produced by program htft\n");
- fprintf(fp,"# %20f :f1,starting frequency [Hz]\n",f1);
- fprintf(fp,"# %20f :f2,ending frequency [Hz]\n",f2);
- fprintf(fp,"# %20f :df,frequency step size [Hz]\n",df);
- fprintf(fp,"# %20d :Nf,number of frequency points at which FT was calculated\n",Nf);
- fprintf(fp,"# %20d :N,number of samp;es read from the input file,N>1\n",N);
- fprintf(fp,"# %20f :sps,samples per second\n",sps);
- fprintf(fp,"# %s :input file name\n",argv[6]);
- fprintf(fp,"# Left column:real part of FT,Right column:imaginary part\n");
- Omega=2.0*M_PI*f1/sps;
- dOmega=2.0*M_PI*df/sps;
- for (i=0;i<Nf;++i)
- {
- hrft (N,q,Omega,&Qr,&Qi);
- fprintf (fp,"%lf %lf\n",Qr,Qi);
- Omega +=dOmega;
- }
- fclose(fp);
- free(q);
- return(0);
- }
复制代码 |