120 lines
2.6 KiB
C
120 lines
2.6 KiB
C
|
#include "math.h"
|
|||
|
|
#include "fft.h"
|
|||
|
|
|
|||
|
|
void conjugate_complex(int n,complex in[],complex out[])
|
|||
|
|
{
|
|||
|
|
int i = 0;
|
|||
|
|
for(i=0;i<n;i++)
|
|||
|
|
{
|
|||
|
|
out[i].imag = -in[i].imag;
|
|||
|
|
out[i].real = in[i].real;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
void c_abs(complex f[],float out[],int n)
|
|||
|
|
{
|
|||
|
|
int i = 0;
|
|||
|
|
float t;
|
|||
|
|
for(i=0;i<n;i++)
|
|||
|
|
{
|
|||
|
|
t = f[i].real * f[i].real + f[i].imag * f[i].imag;
|
|||
|
|
out[i] = iot_sqrt(t);
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
void c_plus(complex a,complex b,complex *c)
|
|||
|
|
{
|
|||
|
|
c->real = a.real + b.real;
|
|||
|
|
c->imag = a.imag + b.imag;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
void c_sub(complex a,complex b,complex *c)
|
|||
|
|
{
|
|||
|
|
c->real = a.real - b.real;
|
|||
|
|
c->imag = a.imag - b.imag;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
void c_mul(complex a,complex b,complex *c)
|
|||
|
|
{
|
|||
|
|
c->real = a.real * b.real - a.imag * b.imag;
|
|||
|
|
c->imag = a.real * b.imag + a.imag * b.real;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
void c_div(complex a,complex b,complex *c)
|
|||
|
|
{
|
|||
|
|
c->real = (a.real * b.real + a.imag * b.imag)/(b.real * b.real +b.imag * b.imag);
|
|||
|
|
c->imag = (a.imag * b.real - a.real * b.imag)/(b.real * b.real +b.imag * b.imag);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
|
|||
|
|
|
|||
|
|
void Wn_i(int n,int i,complex *Wn,char flag)
|
|||
|
|
{
|
|||
|
|
Wn->real = iot_cos(2*PI*i/n);
|
|||
|
|
|
|||
|
|
if(flag == 1)
|
|||
|
|
Wn->imag = -iot_sin(2*PI*i/n);
|
|||
|
|
else if(flag == 0)
|
|||
|
|
Wn->imag = -iot_sin(2*PI*i/n);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
void fft(int N,complex f[])
|
|||
|
|
{
|
|||
|
|
complex t,wn;
|
|||
|
|
int i,j,k,m,n,l,r,M;
|
|||
|
|
int la,lb,lc;
|
|||
|
|
/*----M=log2(N)----*/
|
|||
|
|
for(i=N,M=1;(i=i/2)!=1;M++);
|
|||
|
|
|
|||
|
|
for(i=1,j=N/2;i<=N-2;i++)
|
|||
|
|
{
|
|||
|
|
if(i<j)
|
|||
|
|
{
|
|||
|
|
t=f[j];
|
|||
|
|
f[j]=f[i];
|
|||
|
|
f[i]=t;
|
|||
|
|
}
|
|||
|
|
k=N/2;
|
|||
|
|
while(k<=j)
|
|||
|
|
{
|
|||
|
|
j=j-k;
|
|||
|
|
k=k/2;
|
|||
|
|
}
|
|||
|
|
j=j+k;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/*----FFT----*/
|
|||
|
|
for(m=1;m<=M;m++)
|
|||
|
|
{
|
|||
|
|
la=iot_pow(2,m); //la=2^m<><6D><EFBFBD><EFBFBD><EFBFBD><EFBFBD>m<EFBFBD><6D>ÿ<EFBFBD><C3BF><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><DAB5><EFBFBD>
|
|||
|
|
lb=la/2; //lb<6C><62><EFBFBD><EFBFBD><EFBFBD><EFBFBD>m<EFBFBD><6D>ÿ<EFBFBD><C3BF><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ε<EFBFBD>Ԫ<EFBFBD><D4AA>
|
|||
|
|
//ͬʱ<CDAC><CAB1>Ҳ<EFBFBD><D2B2>ʾÿ<CABE><C3BF><EFBFBD><EFBFBD><EFBFBD>ε<EFBFBD>Ԫ<EFBFBD><D4AA><EFBFBD>½ڵ<C2BD>֮<EFBFBD><D6AE><EFBFBD>ľ<EFBFBD><C4BE><EFBFBD>
|
|||
|
|
/*----<2D><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>----*/
|
|||
|
|
for(l=1;l<=lb;l++)
|
|||
|
|
{
|
|||
|
|
r=(l-1)*iot_pow(2,M-m);
|
|||
|
|
for(n=l-1;n<N-1;n=n+la) //<2F><><EFBFBD><EFBFBD>ÿ<EFBFBD><C3BF><EFBFBD><EFBFBD><EFBFBD>飬<EFBFBD><E9A3AC><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ΪN/la
|
|||
|
|
{
|
|||
|
|
lc=n+lb; //n,lc<6C>ֱ<EFBFBD><D6B1><EFBFBD><EFBFBD><EFBFBD>һ<EFBFBD><D2BB><EFBFBD><EFBFBD><EFBFBD>ε<EFBFBD>Ԫ<EFBFBD><D4AA><EFBFBD>ϡ<EFBFBD><CFA1>½ڵ<C2BD><DAB5><EFBFBD><EFBFBD><EFBFBD>
|
|||
|
|
Wn_i(N,r,&wn,1);//wn=Wnr
|
|||
|
|
c_mul(f[lc],wn,&t);//t = f[lc] * wn<77><6E><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
|||
|
|
c_sub(f[n],t,&(f[lc]));//f[lc] = f[n] - f[lc] * Wnr
|
|||
|
|
c_plus(f[n],t,&(f[n]));//f[n] = f[n] + f[lc] * Wnr
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
//<2F><><EFBFBD><EFBFBD>Ҷ<EFBFBD><D2B6><EFBFBD>任
|
|||
|
|
void ifft(int N,complex f[])
|
|||
|
|
{
|
|||
|
|
int i=0;
|
|||
|
|
conjugate_complex(N,f,f);
|
|||
|
|
fft(N,f);
|
|||
|
|
conjugate_complex(N,f,f);
|
|||
|
|
for(i=0;i<N;i++)
|
|||
|
|
{
|
|||
|
|
f[i].imag = (f[i].imag)/N;
|
|||
|
|
f[i].real = (f[i].real)/N;
|
|||
|
|
}
|
|||
|
|
}
|