#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;
}
评论0