|
|
FFT 以及IFFT 在VC下验证的结果 (原文件名:FFTShow.jpg)
此FFT计算需要依赖复数运算库：
DSP_Complex.h
#ifndef __DSP_COMPLEX_H
#define __DSP_COMPLEX_H
typedef double __ComplexTyp;
typedef struct dsp_complex{
__ComplexTyp real;
__ComplexTyp image;
}DSP_Complex;
/* Add two complex number of A and B result stored in sum*/
void DSP_ComplexADD(DSP_Complex *sum,DSP_Complex *A,DSP_Complex *B);
/*Multiply two complex number of A and B result stored in Mul*/
void DSP_ComplexMUL(DSP_Complex *Mul,DSP_Complex *A,DSP_Complex *B);
#include "DSP_Complex.c"
#endif
DSP_Complex.c
extern void DSP_ComplexADD(DSP_Complex *sum,DSP_Complex *A,DSP_Complex *B)
{
sum->real=A->real+B->real;
sum->image=A->image+B->image;
}
extern void DSP_ComplexSUB(DSP_Complex *sub,DSP_Complex *A,DSP_Complex *B)
{
B->real=-B->real;
B->image=-B->image;
sub->real=A->real+B->real;
sub->image=A->image+B->image;
}
extern void DSP_ComplexMUL(DSP_Complex *Mul,DSP_Complex *A,DSP_Complex *B)
{
Mul->real=(A->real*B->real)-(A->image*B->image);
Mul->image=(A->real*B->image)+(A->image*B->real);
}
DSP_FFT.h
#ifndef __DSP_FFT_H
#define __DSP_FFT_H
typedef double __RFFTv; //real FFT input data type
typedef DSP_Complex __CFFTv; //Image TTF input data type
/*********************************************************************************
Functionaility :2-N based Real number input array FFT with result stored in FFTout
Input:
__RFFT *Data Pointer pointed to the real input array
int PointNum number of element in input array must to be the integer power of 2
DSP_Complex *FFTout Transform result
*********************************************************************************/
void DSP_FFT_1D2N_real(__RFFTv *Data,int PointNum,DSP_Complex *FFTout);
/*********************************************************************************
Functionaility :2-N based complex number input array FFT with result stored in FFTout
Input:
__CFFTv *Data Pointer pointed to the real input array
int PointNum number of element in input array must to be the integer power of 2
DSP_Complex *FFTout Transform result
*********************************************************************************/
void DSP_FFT_1D_Complex(__CFFTv *Data,int PointNum,DSP_Complex *FFTout)
/*********************************************************************************
Functionaility :2-N based complex number input array IFFT with result stored in FFTout
Input:
__CFFTv *Data Pointer pointed to the real input array
int PointNum number of element in input array must to be the integer power of 2
DSP_Complex *FFTout Transform result
*********************************************************************************/
void DSP_IFFT_1D_Complex(__CFFTv *Data,int PointNum,DSP_Complex *FFTout);
#include "DSP_FFT.c"
#endif
DSP_FFT.c
#include <math.h>
#define pi 3.1416
static double __cos(double in)
{
return cos(in);
}
static double __sin(double in)
{
return sin(in);
}
static void RDFT_2Point(__RFFTv *dataIn,DSP_Complex *DFTout)
{
DFTout->real = *dataIn+*(dataIn+1);
DFTout->image = 0;
DFTout++;
DFTout->real = *dataIn+(*(dataIn+1))*__cos(pi);
DFTout->image = *(dataIn+1)*-__sin(pi);
}
inline void RExchange(__RFFTv *A, __RFFTv *B)
{
__RFFTv buf;
buf=*A;
*A=*B;
*B=buf;
}
/*将输入的2 n次方个元素数的数组的奇数下脚标的元素与偶数下脚标的元素分开 偶数在前奇数在后*/
static void RAdjustOrder(__RFFTv *data,int Num)
{
int i,j,cycle;
if(Num!=4){
i=Num/2;
j=i/2;
RAdjustOrder(data,i);
RAdjustOrder(data+i,i);
for(cycle=0;cycle<j;cycle++){
RExchange(data+j+cycle,data+i+cycle);
}
}
else{
RExchange(data+1,data+2);
}
}
extern void DSP_FFT_1D_real(__RFFTv *Data,int PointNum,DSP_Complex *FFTout)
{
int i,cycle;
double angle;
DSP_Complex shift,buf1,buf2,*X2;
if(PointNum==2){
RDFT_2Point(Data,FFTout);
}
else{
angle=(2.0*pi)/(double)PointNum;
i=PointNum/2;
X2=FFTout+i;
RAdjustOrder(Data,PointNum);
DSP_FFT_1D_real(Data,i,FFTout);
DSP_FFT_1D_real(Data+i,i,X2);
for(cycle=0;cycle<i;cycle++,X2++,FFTout++){
shift.real=cos(angle*(double)cycle);
shift.image=-sin(angle*(double)cycle);
DSP_ComplexMUL(&buf1,&shift,X2);
DSP_ComplexADD(&buf2,FFTout,&buf1);
DSP_ComplexSUB(X2,FFTout,&buf1);
(FFTout)->image=buf2.image;
(FFTout)->real=buf2.real;
}
}
}
/****************************************************************************************************************/
static void CDFT_2Point(__CFFTv *dataIn,DSP_Complex *DFTout)
{
DSP_Complex buf,W_1;
DSP_ComplexADD(DFTout,dataIn,dataIn+1);
DFTout++;
W_1.real=__cos(pi);
W_1.image=-__sin(pi);
DSP_ComplexMUL(&buf,dataIn+1,&W_1);
DSP_ComplexADD(DFTout,dataIn,&buf);
}
inline void CExchange(__CFFTv *A, __CFFTv *B)
{
__CFFTv buf;
buf=*A;
*A=*B;
*B=buf;
}
/*将输入的2 n次方个元素数的数组的奇数下脚标的元素与偶数下脚标的元素分开 偶数在前奇数在后*/
static void CAdjustOrder(__CFFTv *data,int Num)
{
int i,j,cycle;
if(Num!=4){
i=Num/2;
j=i/2;
CAdjustOrder(data,i);
CAdjustOrder(data+i,i);
for(cycle=0;cycle<j;cycle++){
CExchange(data+j+cycle,data+i+cycle);
}
}
else{
CExchange(data+1,data+2);
}
}
extern void DSP_FFT_1D_Complex(__CFFTv *Data,int PointNum,DSP_Complex *FFTout)
{
int i,cycle;
double angle;
DSP_Complex shift,buf1,buf2,*X2;
if(PointNum==2){
CDFT_2Point(Data,FFTout);
}
else{
angle=(2.0*pi)/(double)PointNum;
i=PointNum/2;
X2=FFTout+i;
CAdjustOrder(Data,PointNum);
DSP_FFT_1D_Complex(Data,i,FFTout);
DSP_FFT_1D_Complex(Data+i,i,X2);
for(cycle=0;cycle<i;cycle++,X2++,FFTout++){
shift.real=cos(angle*(double)cycle);
shift.image=-sin(angle*(double)cycle);
DSP_ComplexMUL(&buf1,&shift,X2);
DSP_ComplexADD(&buf2,FFTout,&buf1);
DSP_ComplexSUB(X2,FFTout,&buf1);
(FFTout)->image=buf2.image;
(FFTout)->real=buf2.real;
}
}
}
extern void DSP_IFFT_1D_Complex(__CFFTv *Data,int PointNum,DSP_Complex *FFTout)
{
int i;
__CFFTv *Delegate;
Delegate=Data;
for(i=0;i<PointNum;i++,Delegate++)
Delegate->image=-Delegate->image;
DSP_FFT_1D_Complex(Data,PointNum,FFTout);
Delegate=FFTout;
for(i=0;i<PointNum;i++,Delegate++){
Delegate->image/=PointNum;
Delegate->real/=PointNum;
}
} |
阿莫论坛20周年了！感谢大家的支持与爱护！！
一只鸟敢站在脆弱的枝条上歇脚，它依仗的不是枝条不会断，而是自己有翅膀，会飞。
|
发表于 2012-2-18 19:21:18
