-
</pre><p>快速离散傅立叶变换FFT利用DFT计算的对称性实现的,具体的介绍网上一大堆。这次自己写了个定点FFT头文件，直接用C语言写的很容易移植。</p><p></p><p><pre name=
"code" class=
"cpp">
/*
-
快速离散傅立叶算法V1.0
-
含有：FFT,IFFT
-
made by xyt
-
2015/7/8
-
C语言
-
*/
-
#ifndef _FFT_H
-
#define _FFT_H
-
#include<math.h>
-
-
#define FFT_N
8
//点数,要求2的次方
-
-
#define PI
3.14159265354
-
struct fft_complex{
-
double r,i;
-
};
-
int fft_fi(double in){
//四舍五入转整
-
if((in-(
int)in)>
0.5)
return (
int)in+
1;
-
else
return (
int)in;
-
}
-
int fft_fac2(int n){
//返回log2(n)
-
int count=
0;
-
while(n/
2!=
0){
-
n/=
2;count++;
-
}
-
return count;
-
}
-
int fft_turndat(int n,int num){
-
int g=n,m,r=
0;
-
int count=
0;
-
while(num/
2!=
0){
-
num/=
2;count++;
-
}
-
while(count>=
0){
-
m=g%
2;
-
r+=m*
pow(
2,--count);
-
g/=
2;
-
}
-
return r;
-
}
-
fft_complex fft_t(fft_complex a){
-
fft_complex tmp;
-
tmp.i=
-1*a.i;
-
tmp.r=
-1*a.r;
-
return tmp;
-
}
-
fft_complex fft_multi(fft_complex a,fft_complex b){
-
fft_complex tmp;
-
tmp.r=a.r*b.r-a.i*b.i;
-
tmp.i=a.r*b.i+a.i*b.r;
-
return tmp;
-
}
-
fft_complex fft_add(fft_complex a,fft_complex b){
-
fft_complex tmp;
-
tmp.r=a.r+b.r;
-
tmp.i=a.i+b.i;
-
return tmp;
-
}
-
/* 定点FFT，返回两个double型数组分别是实部和虚部 */
-
void FFT(int *in,double *outr,double *outi)
-
{
-
int i,j;
-
int deep;
-
const
int N=FFT_N/
2;
-
fft_complex W[N];
-
W[
0].r=
1;W[
0].i=
0;
-
W[
1].r=
cos(
2.0*PI/FFT_N);
-
W[
1].i=
sin(
2.0*PI/FFT_N);
-
for(i=
2;i<N;i++){
-
W[i]=fft_multi(W[
1],W[i
-1]);
-
}
-
deep=fft_fac2(FFT_N);
-
int g=
1;
-
int ne=
0,ge=
0;
-
-
fft_complex left[FFT_N];
-
fft_complex right[FFT_N];
-
for(i=
0;i<FFT_N;i++){
-
left[i].r=in[fft_turndat(i,FFT_N)];
-
left[i].i=
0;
-
}
-
fft_complex tpp;
-
int mggtmp;
-
while(
1)
-
{
-
if(deep==
0)
break;
-
int adt=
pow(
2,deep
-1);
-
mggtmp=
pow(
2,g);
-
for(i=
0;i<FFT_N;i+=mggtmp){
-
ne=
0;ge=
0;
-
for(j=
0;j<mggtmp;j++){
-
if(j<mggtmp/
2){
-
tpp=fft_multi(left[i+j+mggtmp/
2],W[ne]);
-
right[i+j]=fft_add(left[i+j],tpp);ne+=adt;
-
}
else{
-
tpp=
fft_t(fft_multi(left[i+j],W[ge]));
-
right[i+j]=fft_add(left[i+j-mggtmp/
2],tpp);ge+=adt;
-
}
-
}
-
}
-
for(i=
0;i<FFT_N;i++){
-
left[i]=right[i];
-
}
-
deep--;g++;
-
}
-
for(i=
0;i<FFT_N;i++){
-
outr[i]=left[i].r;
-
outi[i]=
-1*left[i].i;
-
}
-
}
-
/* 定点IFFT，输入实部和虚部,返回时域int类型 */
-
void IFFT(double *inr,double *ini,int *out)
-
{
-
int i,j;
-
int deep;
-
const
int N=FFT_N/
2;
-
fft_complex W[N];
-
W[
0].r=
1;W[
0].i=
0;
-
W[
1].r=
cos(
2.0*PI/FFT_N);
-
W[
1].i=
-1*
sin(
2.0*PI/FFT_N);
-
for(i=
2;i<N;i++){
-
W[i]=fft_multi(W[
1],W[i
-1]);
-
}
-
deep=fft_fac2(FFT_N);
-
int g=
1;
-
int ne=
0,ge=
0;
-
-
fft_complex left[FFT_N];
-
fft_complex right[FFT_N];
-
for(i=
0;i<FFT_N;i++){
-
left[i].r=inr[fft_turndat(i,FFT_N)];
-
left[i].i=
-1*ini[fft_turndat(i,FFT_N)];
-
}
-
fft_complex tpp;
-
int mggtmp;
-
while(
1)
-
{
-
if(deep==
0)
break;
-
int adt=
pow(
2,deep
-1);
-
mggtmp=
pow(
2,g);
-
for(i=
0;i<FFT_N;i+=mggtmp){
-
ne=
0;ge=
0;
-
for(j=
0;j<mggtmp;j++){
-
if(j<mggtmp/
2){
-
tpp=fft_multi(left[i+j+mggtmp/
2],W[ne]);
-
right[i+j]=fft_add(left[i+j],tpp);ne+=adt;
-
}
else{
-
tpp=
fft_t(fft_multi(left[i+j],W[ge]));
-
right[i+j]=fft_add(left[i+j-mggtmp/
2],tpp);ge+=adt;
-
}
-
}
-
}
-
for(i=
0;i<FFT_N;i++){
-
left[i]=right[i];
-
}
-
deep--;g++;
-
}
-
for(i=
0;i<FFT_N;i++){
-
out[i]=fft_fi(left[i].r/FFT_N);
-
}
-
}
-
#endif
FFT快速傅立叶算法纯C语言版本（转）
最新推荐文章于 2021-11-11 04:45:39 发布
7214
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
