794ba8da0bd777ab47e8622078940722.rar_CZT_FFT实现CZT_czt C语言_czt实现_

#include <stdio.h> #include <stdlib.h> #include <math.h> #include "msp.h" void mcztfft(complex* x1,int n,int m,int maxnm,float dltomg, float omg0,float fs,int *ierror) { /*---------------------------------------------------------------------- Routine CZTFFT:To perform a CHIRP-Z Transform. input parameters: N : numbers of input signal to be transformed by CZT, x(i),i=0,..,N-1; M : numbers of output signal transformed by CZT, x(i),i=0,..,M-1; ! Note !: the dimension of x must be MAXNM,otherwise the result will be error! dltomg: increasing angle of frequncy domain, in Hz; omg0 : start angle of CZT, in Hz; fs : sampling frequency in Hz. output parameters: x : MAXNM dimensioned complex array. the CZT result is stored in x(0) to x(M-1). Note : In this program the largest length to do CZT MAXCZT=1024, if you want to increase the length to do CZT please change MAXCZT and the dimension of X,WORK1 and WORK2 . in Chapter 5 --------------------------------------------------------------------*/ complex a,z; complex x[1024],work1[1024],work2[1024]; float pi,ph,alpha,ang; int nfft,i,ntest,maxczt,m1; m1=max(n,m); maxczt=1024; *ierror=1; if(maxnm<m1) return; /*-------------------------------- TO COMPUTE W**(-N**2/2)-----*/ nfft=n+m; for(i=0;i<=10;i++) {ntest=pow(2,i); if(ntest>=nfft) break; } nfft=ntest; *ierror=2; if(nfft>maxczt) return; pi=4.*atan(1.); ph=-2.*pi*dltomg/fs; for(i=0;i<maxnm;i++) {alpha=-ph*i*i/2.; work1[i].real=cos(alpha); work1[i].imag=sin(alpha); work2[i].real=work1[i].real; work2[i].imag=work1[i].imag; } for(i=1;i<n;i++) {work1[nfft-i].real=work2[i].real; work1[nfft-i].imag=work2[i].imag; } if(n+m!=nfft) for(i=m;i<=nfft-n;i++) {work1[i].real=0.; work1[i].imag=0.; } mcmpfft(work1,nfft,-1); /*------------------------ TO COMPUTE A**(-N)*W**(N**2/2)-----*/ ang=-omg0*pi*2./fs; for(i=0;i<n;i++) {a.real=cos(ang*i); a.imag=sin(ang*i); z.real=work2[i].real*a.real+work2[i].imag*a.imag; z.imag=work2[i].real*a.imag-work2[i].imag*a.real; x[i].real=z.real*x1[i].real-z.imag*x1[i].imag; x[i].imag=z.real*x1[i].imag+z.imag*x1[i].real; } for(i=n;i<nfft;i++) {x[i].real=0.; x[i].imag=0.; } mcmpfft(x,nfft,-1); /*------------------------ TO COMPUTE g(n)*h(n)---------------*/ for(i=0;i<nfft;i++) {z.real=x[i].real; z.imag=x[i].imag; x[i].real=z.real*work1[i].real-z.imag*work1[i].imag; x[i].imag=z.real*work1[i].imag+z.imag*work1[i].real; } mcmpfft(x,nfft,1); /*-------------------------- TO MULTIPLY BY W**(K**2/2)--------*/ for(i=0;i<m;i++) {z.real=x[i].real; z.imag=x[i].imag; a.real=work2[i].real; a.imag=-work2[i].imag; x1[i].real=z.real*a.real-z.imag*a.imag; x1[i].imag=z.imag*a.real+z.real*a.imag; } return; }