本文目录一览：
- 1、用c语言实现FFT
- 2、matlab fft2的c代码
- 3、基于FFT的算法优化 要C语言完整程序（利用旋转因子的性质），有的请留言，答谢！！！(有核心代码，望指教
- 4、怎样用C语言实现FFT算法啊？
- 5、求FFT的C语言程序……最好是1024点的……希望大家帮帮我！
- 6、C语言 1024点快速傅里叶变换（FFT）程序，最好经过优化，执行速度快
用c语言实现FFT
float ar[1024],ai[1024];/* 原始数据实部，虚部 */
float a[2050];
void fft(int nn) /* nn数据长度 */
{
int n1,n2,i,j,k,l,m,s,l1;
float t1,t2,x,y;
float w1,w2,u1,u2,z;
float fsin[10]={0.000000,1.000000,0.707107,0.3826834,0.1950903,0.09801713,0.04906767,0.02454123,0.01227154,0.00613588,};
float fcos[10]={-1.000000,0.000000,0.7071068,0.9238796,0.9807853,0.99518472,0.99879545,0.9996988,0.9999247,0.9999812,};
switch(nn)
{
case 1024: s=10; break;
case 512: s=9; break;
case 256: s=8; break;
}
n1=nn/2; n2=nn-1;
j=1;
for(i=1;i=nn;i++)
{
a[2*i]=ar[i-1];
a[2*i+1]=ai[i-1];
}
for(l=1;ln2;l++)
{
if(lj)
{
t1=a[2*j];
t2=a[2*j+1];
a[2*j]=a[2*l];
a[2*j+1]=a[2*l+1];
a[2*l]=t1;
a[2*l+1]=t2;
}
k=n1;
while (kj)
{
j=j-k;
k=k/2;
}
j=j+k;
}
for(i=1;i=s;i++)
{
u1=1;
u2=0;
m=(1i);
k=m1;
w1=fcos[i-1];
w2=-fsin[i-1];
for(j=1;j=k;j++)
{
for(l=j;lnn;l=l+m)
{
l1=l+k;
t1=a[2*l1]*u1-a[2*l1+1]*u2;
t2=a[2*l1]*u2+a[2*l1+1]*u1;
a[2*l1]=a[2*l]-t1;
a[2*l1+1]=a[2*l+1]-t2;
a[2*l]=a[2*l]+t1;
a[2*l+1]=a[2*l+1]+t2;
}
z=u1*w1-u2*w2;
u2=u1*w2+u2*w1;
u1=z;
}
}
for(i=1;i=nn/2;i++)
{
ar[i]=4*a[2*i+2]/nn; /* 实部 */
ai[i]=-4*a[2*i+3]/nn; /* 虚部 */
a[i]=4*sqrt(ar[i]*ar[i]+ai[i]*ai[i]); /* 幅值 */
}
}
(;si=2)
打字不易，如满意，望采纳。
matlab fft2的c代码
傅立叶变换的c语言源代码
128点DIT FFT函数：
/* 采样来的数据放在dataR[ ]数组中，运算前dataI[ ]数组初始化为0 */
void FFT(float dataR[],float dataI[])
{int x0,x1,x2,x3,x4,x5,x6;
int L,j,k,b,p;
float TR,TI,temp;
/********** following code invert sequence ************/
for(i=0;i128;i++)
{ x0=x1=x2=x3=x4=x5=x6=0;
x0=i0x01; x1=(i/2)0x01; x2=(i/4)0x01; x3=(i/8)0x01;x4=(i/16)0x01; x5=(i/32)0x01; x6=(i/64)0x01;
xx=x0*64+x1*32+x2*16+x3*8+x4*4+x5*2+x6;
dataI[xx]=dataR[i];
}
for(i=0;i128;i++)
{ dataR[i]=dataI[i]; dataI[i]=0; }
/************** following code FFT *******************/
for(L=1;L=7;L++) { /* for(1) */
b=1; i=L-1;
while(i0)
{b=b*2; i--；} /* b= 2^(L-1) */
for(j=0;j=b-1;j++) /* for (2) */
{ p=1; i=7-L;
while(i0) /* p=pow(2,7-L)*j; */
{p=p*2; i--；}
p=p*j;
for(k=j;k128;k=k+2*b) /* for (3) */
{ TR=dataR[k]; TI=dataI[k]; temp=dataR[k+b];
dataR[k]=dataR[k]+dataR[k+b]*cos_tab[p]+dataI[k+b]*sin_tab[p];
dataI[k]=dataI[k]-dataR[k+b]*sin_tab[p]+dataI[k+b]*cos_tab[p];
dataR[k+b]=TR-dataR[k+b]*cos_tab[p]-dataI[k+b]*sin_tab[p];
dataI[k+b]=TI+temp*sin_tab[p]-dataI[k+b]*cos_tab[p];
} /* END for (3) */
} /* END for (2) */
} /* END for (1) */
for(i=0;i32;i++){ /* 只需要32次以下的谐波进行分析 */
w[i]=sqrt(dataR[i]*dataR[i]+dataI[i]*dataI[i]);
w[i]=w[i]/64；}
w[0]=w[0]/2;
} /* END FFT */
基于FFT的算法优化 要C语言完整程序（利用旋转因子的性质），有的请留言，答谢！！！(有核心代码，望指教
实现(C描述)
#include stdio.h
#include math.h
#include stdlib.h
//#include "complex.h"
// --------------------------------------------------------------------------
#define N 8 //64
#define M 3 //6 //2^m=N
#define PI 3.1415926
// --------------------------------------------------------------------------
float twiddle[N/2] = {1.0, 0.707, 0.0, -0.707};
float x_r[N] = {1, 1, 1, 1, 0, 0, 0, 0};
float x_i[N]; //N=8
/*
float twiddle[N/2] = {1, 0.9951, 0.9808, 0.9570, 0.9239, 0.8820, 0.8317, 0.7733,
0.7075, 0.6349, 0.5561, 0.4721, 0.3835, 0.2912, 0.1961, 0.0991,
0.0000,-0.0991,-0.1961,-0.2912,-0.3835,-0.4721,-0.5561,-0.6349,
-0.7075,-0.7733, 0.8317,-0.8820,-0.9239,-0.9570,-0.9808,-0.9951}; //N=64
float x_r[N]={1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,};
float x_i[N];
*/
FILE *fp;
// ----------------------------------- func -----------------------------------
/**
* 初始化输出虚部
*/
static void fft_init( void )
{
int i;
for(i=0; iN; i++) x_i[i] = 0.0;
}
/**
* 反转算法.将时域信号重新排序.
* 这个算法有改进的空间
*/
static void bitrev( void )
{
int p=1, q, i;
int bit_rev[ N ]; //
float xx_r[ N ]; //
bit_rev[ 0 ] = 0;
while( p N )
{
for(q=0; qp; q++)
{
bit_rev[ q ] = bit_rev[ q ] * 2;
bit_rev[ q + p ] = bit_rev[ q ] + 1;
}
p *= 2;
}
for(i=0; iN; i++) xx_r[ i ] = x_r[ i ];
for(i=0; iN; i++) x_r[i] = xx_r[ bit_rev[i] ];
}
/* ------------ add by sshc625 ------------ */
static void bitrev2( void )
{
return ;
}
/* */
void display( void )
{
printf("nn");
int i;
for(i=0; iN; i++)
printf("%ft%fn", x_r[i], x_i[i]);
}
/**
*
*/
void fft1( void )
{ fp = fopen("log1.txt", "a+");
int L, i, b, j, p, k, tx1, tx2;
float TR, TI, temp; // 临时变量
float tw1, tw2;
/* 深M. 对层进行循环. L为当前层, 总层数为M. */
for(L=1; L=M; L++)
{
fprintf(fp,"----------Layer=%d----------n", L);
/* b的意义非常重大,b表示当前层的颗粒具有的输入样本点数 */
b = 1;
i = L - 1;
while(i 0)
{
b *= 2;
i--;
}
// -------------- 是否外层对颗粒循环, 内层对样本点循环逻辑性更强一些呢! --------------
/*
* outter对参与DFT的样本点进行循环
* L=1, 循环了1次(4个颗粒, 每个颗粒2个样本点)
* L=2, 循环了2次(2个颗粒, 每个颗粒4个样本点)
* L=3, 循环了4次(1个颗粒, 每个颗粒8个样本点)
*/
for(j=0; jb; j++)
{
/* 求旋转因子tw1 */
p = 1;
i = M - L; // M是为总层数, L为当前层.
while(i 0)
{
p = p*2;
i--;
}
p = p * j;
tx1 = p % N;
tx2 = tx1 + 3*N/4;
tx2 = tx2 % N;
// tw1是cos部分, 实部; tw2是sin部分, 虚数部分.
tw1 = ( tx1=N/2)? -twiddle[tx1-N/2] : twiddle[ tx1 ];
tw2 = ( tx2=N/2)? -twiddle[tx2-(N/2)] : twiddle[tx2];
/*
* inner对颗粒进行循环
* L=1, 循环了4次(4个颗粒, 每个颗粒2个输入)
* L=2, 循环了2次(2个颗粒, 每个颗粒4个输入)
* L=3, 循环了1次(1个颗粒, 每个颗粒8个输入)
*/
for(k=j; kN; k=k+2*b)
{
TR = x_r[k]; // TR就是A, x_r[k+b]就是B.
TI = x_i[k];
temp = x_r[k+b];
/*
* 如果复习一下 (a+j*b)(c+j*d)两个复数相乘后的实部虚部分别是什么
* 就能理解为什么会如下运算了, 只有在L=1时候输入才是实数, 之后层的
* 输入都是复数, 为了让所有的层的输入都是复数, 我们只好让L=1时候的
* 输入虚部为0
* x_i[k+b]*tw2是两个虚数相乘
*/
fprintf(fp, "tw1=%f, tw2=%fn", tw1, tw2);
x_r[k] = TR + x_r[k+b]*tw1 + x_i[k+b]*tw2;
x_i[k] = TI - x_r[k+b]*tw2 + x_i[k+b]*tw1;
x_r[k+b] = TR - x_r[k+b]*tw1 - x_i[k+b]*tw2;
x_i[k+b] = TI + temp*tw2 - x_i[k+b]*tw1;
fprintf(fp, "k=%d, x_r[k]=%f, x_i[k]=%fn", k, x_r[k], x_i[k]);
fprintf(fp, "k=%d, x_r[k]=%f, x_i[k]=%fn", k+b, x_r[k+b], x_i[k+b]);
} //
} //
} //
}
/**
* ------------ add by sshc625 ------------
* 该实现的流程为
* for( Layer )
* for( Granule )
* for( Sample )
*
*
*
*
*/
void fft2( void )
{ fp = fopen("log2.txt", "a+");
int cur_layer, gr_num, i, k, p;
float tmp_real, tmp_imag, temp; // 临时变量, 记录实部
float tw1, tw2;// 旋转因子,tw1为旋转因子的实部cos部分, tw2为旋转因子的虚部sin部分.
int step; // 步进
int sample_num; // 颗粒的样本总数(各层不同, 因为各层颗粒的输入不同)
/* 对层循环 */
for(cur_layer=1; cur_layer=M; cur_layer++)
{
/* 求当前层拥有多少个颗粒(gr_num) */
gr_num = 1;
i = M - cur_layer;
while(i 0)
{
i--;
gr_num *= 2;
}
/* 每个颗粒的输入样本数N' */
sample_num = (int)pow(2, cur_layer);
/* 步进. 步进是N'/2 */
step = sample_num/2;
/* */
k = 0;
/* 对颗粒进行循环 */
for(i=0; igr_num; i++)
{
/*
* 对样本点进行循环, 注意上限和步进
*/
for(p=0; psample_num/2; p++)
{
// 旋转因子, 需要优化...
tw1 = cos(2*PI*p/pow(2, cur_layer));
tw2 = -sin(2*PI*p/pow(2, cur_layer));
tmp_real = x_r[k+p];
tmp_imag = x_i[k+p];
temp = x_r[k+p+step];
/*(tw1+jtw2)(x_r[k]+jx_i[k])
*
* real : tw1*x_r[k] - tw2*x_i[k]
* imag : tw1*x_i[k] + tw2*x_r[k]
* 我想不抽象出一个
* typedef struct {
* double real; // 实部
* double imag; // 虚部
* } complex; 以及针对complex的操作
* 来简化复数运算是否是因为效率上的考虑!
*/
/* 蝶形算法 */
x_r[k+p] = tmp_real + ( tw1*x_r[k+p+step] - tw2*x_i[k+p+step] );
x_i[k+p] = tmp_imag + ( tw2*x_r[k+p+step] + tw1*x_i[k+p+step] );
/* X[k] = A(k)+WB(k)
* X[k+N/2] = A(k)-WB(k) 的性质可以优化这里*/
// 旋转因子, 需要优化...
tw1 = cos(2*PI*(p+step)/pow(2, cur_layer));
tw2 = -sin(2*PI*(p+step)/pow(2, cur_layer));
x_r[k+p+step] = tmp_real + ( tw1*temp - tw2*x_i[k+p+step] );
x_i[k+p+step] = tmp_imag + ( tw2*temp + tw1*x_i[k+p+step] );
printf("k=%d, x_r[k]=%f, x_i[k]=%fn", k+p, x_r[k+p], x_i[k+p]);
printf("k=%d, x_r[k]=%f, x_i[k]=%fn", k+p+step, x_r[k+p+step], x_i[k+p+step]);
}
/* 开跳!:) */
k += 2*step;
}
}
}
/*
* 后记:
* 究竟是颗粒在外层循环还是样本输入在外层, 好象也差不多, 复杂度完全一样.
* 但以我资质愚钝花费了不少时间才弄明白这数十行代码.
* 从中我发现一个于我非常有帮助的教训, 很久以前我写过一部分算法, 其中绝大多数都是递归.
* 将数据量减少, 减少再减少, 用归纳的方式来找出数据量加大代码的规律
* 比如FFT
* 1. 先写死LayerI的代码; 然后再把LayerI的输出作为LayerII的输入, 又写死代码; ......
* 大约3层就可以统计出规律来. 这和递归也是一样, 先写死一两层, 自然就出来了!
* 2. 有的功能可以写伪代码, 不急于求出结果, 降低复杂性, 把逻辑结果定出来后再添加.
* 比如旋转因子就可以写死, 就写1.0. 流程出来后再写旋转因子.
* 寥寥数语, 我可真是流了不少汗! Happy!
*/
void dft( void )
{
int i, n, k, tx1, tx2;
float tw1,tw2;
float xx_r[N],xx_i[N];
/*
* clear any data in Real and Imaginary result arrays prior to DFT
*/
for(k=0; k=N-1; k++)
xx_r[k] = xx_i[k] = x_i[k] = 0.0;
// caculate the DFT
for(k=0; k=(N-1); k++)
{
for(n=0; n=(N-1); n++)
{
tx1 = (n*k);
tx2 = tx1+(3*N)/4;
tx1 = tx1%(N);
tx2 = tx2%(N);
if(tx1 = (N/2))
tw1 = -twiddle[tx1-(N/2)];
else
tw1 = twiddle[tx1];
if(tx2 = (N/2))
tw2 = -twiddle[tx2-(N/2)];
else
tw2 = twiddle[tx2];
xx_r[k] = xx_r[k]+x_r[n]*tw1;
xx_i[k] = xx_i[k]+x_r[n]*tw2;
}
xx_i[k] = -xx_i[k];
}
// display
for(i=0; iN; i++)
printf("%ft%fn", xx_r[i], xx_i[i]);
}
// ---------------------------------------------------------------------------
int main( void )
{
fft_init( );
bitrev( );
// bitrev2( );
//fft1( );
fft2( );
display( );
system( "pause" );
// dft();
return 1;
}
本文来自CSDN博客，转载请标明出处：
怎样用C语言实现FFT算法啊？
1、二维FFT相当于对行和列分别进行一维FFT运算。具体的实现办法如下：
先对各行逐一进行一维FFT，然后再对变换后的新矩阵的各列逐一进行一维FFT。相应的伪代码如下所示：
for (int i=0; iM; i++)
FFT_1D(ROW[i]，N);
for (int j=0; jN; j++)
FFT_1D(COL[j]，M);
其中，ROW[i]表示矩阵的第i行。注意这只是一个简单的记法，并不能完全照抄。还需要通过一些语句来生成各行的数据。同理，COL[i]是对矩阵的第i列的一种简单表示方法。
所以，关键是一维FFT算法的实现。
2、例程：
#include stdio.h
#include math.h
#include stdlib.h
#define N 1000
/*定义复数类型*/
typedef struct{
double real;
double img;
}complex;
complex x[N], *W; /*输入序列,变换核*/
int size_x=0; /*输入序列的大小，在本程序中仅限2的次幂*/
double PI; /*圆周率*/
void fft(); /*快速傅里叶变换*/
void initW(); /*初始化变换核*/
void change(); /*变址*/
void add(complex ,complex ,complex *); /*复数加法*/
void mul(complex ,complex ,complex *); /*复数乘法*/
void sub(complex ,complex ,complex *); /*复数减法*/
void output();
int main(){
int i; /*输出结果*/
system("cls");
PI=atan(1)*4;
printf("Please input the size of x:n");
scanf("%d",size_x);
printf("Please input the data in x[N]:n");
for(i=0;isize_x;i++)
scanf("%lf%lf",x[i].real,x[i].img);
initW();
fft();
output();
return 0;
}
/*快速傅里叶变换*/
void fft(){
int i=0,j=0,k=0,l=0;
complex up,down,product;
change();
for(i=0;i log(size_x)/log(2) ;i++){ /*一级蝶形运算*/
l=1i;
for(j=0;jsize_x;j+= 2*l ){ /*一组蝶形运算*/
for(k=0;kl;k++){ /*一个蝶形运算*/
mul(x[j+k+l],W[size_x*k/2/l],product);
add(x[j+k],product,up);
sub(x[j+k],product,down);
x[j+k]=up;
x[j+k+l]=down;
}
}
}
}
/*初始化变换核*/
void initW(){
int i;
W=(complex *)malloc(sizeof(complex) * size_x);
for(i=0;isize_x;i++){
W[i].real=cos(2*PI/size_x*i);
W[i].img=-1*sin(2*PI/size_x*i);
}
}
/*变址计算，将x(n)码位倒置*/
void change(){
complex temp;
unsigned short i=0,j=0,k=0;
double t;
for(i=0;isize_x;i++){
k=i;j=0;
t=(log(size_x)/log(2));
while( (t--)0 ){
j=j1;
j|=(k 1);
k=k1;
}
if(ji){
temp=x[i];
x[i]=x[j];
x[j]=temp;
}
}
}
/*输出傅里叶变换的结果*/
void output(){
int i;
printf("The result are as followsn");
for(i=0;isize_x;i++){
printf("%.4f",x[i].real);
if(x[i].img=0.0001)printf("+%.4fjn",x[i].img);
else if(fabs(x[i].img)0.0001)printf("n");
else printf("%.4fjn",x[i].img);
}
}
void add(complex a,complex b,complex *c){
c-real=a.real+b.real;
c-img=a.img+b.img;
}
void mul(complex a,complex b,complex *c){
c-real=a.real*b.real - a.img*b.img;
c-img=a.real*b.img + a.img*b.real;
}
void sub(complex a,complex b,complex *c){
c-real=a.real-b.real;
c-img=a.img-b.img;
}
求FFT的C语言程序……最好是1024点的……希望大家帮帮我！
float ar[1024],ai[1024];/* 实部，虚部 */
float a[2050]; /* 实际值 */
void fft()
{
int n1,n2,i,j,k,l,m,s=10,nn=1024,l1;
float t1,t2,x,y;
float w1,w2,u1,u2,z;
float fsin[10]={0.000000,1.000000,0.707107,0.3826834,0.1950903,0.09801713,0.04906767,0.02454123,0.01227154,0.00613588,};
float fcos[10]={-1.000000,0.000000,0.7071068,0.9238796,0.9807853,0.99518472,0.99879545,0.9996988,0.9999247,0.9999812,};
n1=nn/2; n2=nn-1;
j=1;
for(i=1;i=nn;i++)
{
a[2*i]=ar[i-1];
a[2*i+1]=ai[i-1];
}
for(l=1;ln2;l++)
{
if(lj)
{
t1=a[2*j];
t2=a[2*j+1];
a[2*j]=a[2*l];
a[2*j+1]=a[2*l+1];
a[2*l]=t1;
a[2*l+1]=t2;
}
k=n1;
while (kj)
{
j=j-k;
k=k/2;
}
j=j+k;
}
for(i=1;i=s;i++)
{
u1=1;
u2=0;
m=(1i);
k=m1;
w1=fcos[i-1];
w2=-fsin[i-1];
for(j=1;j=k;j++)
{
for(l=j;lnn;l=l+m)
{
l1=l+k;
t1=a[2*l1]*u1-a[2*l1+1]*u2;
t2=a[2*l1]*u2+a[2*l1+1]*u1;
a[2*l1]=a[2*l]-t1;
a[2*l1+1]=a[2*l+1]-t2;
a[2*l]=a[2*l]+t1;
a[2*l+1]=a[2*l+1]+t2;
}
z=u1*w1-u2*w2;
u2=u1*w2+u2*w1;
u1=z;
}
}
for(i=1;i=nn/2;i++)
{
ar[i]=a[2*i+2]/nn;
ai[i]=-a[2*i+3]/nn;
a[i]=4*sqrt(ar[i]*ar[i]+ai[i]*ai[i]);
}
}
C语言 1024点快速傅里叶变换（FFT）程序，最好经过优化，执行速度快
void fft()
{
int nn,n1,n2,i,j,k,l,m,s,l1;
float ar[1024],ai[1024]; // 实部 虚部
float a[2050];
float t1,t2,x,y;
float w1,w2,u1,u2,z;
float fsin[10]={0.000000,1.000000,0.707107,0.3826834,0.1950903,0.09801713,0.04906767,0.02454123,0.01227154,0.00613588,};// 优化
float fcos[10]={-1.000000,0.000000,0.7071068,0.9238796,0.9807853,0.99518472,0.99879545,0.9996988,0.9999247,0.9999812,};
nn=1024;
s=10;
n1=nn/2; n2=nn-1;
j=1;
for(i=1;i=nn;i++)
{
a[2*i]=ar[i-1];
a[2*i+1]=ai[i-1];
}
for(l=1;ln2;l++)
{
if(lj)
{
t1=a[2*j];
t2=a[2*j+1];
a[2*j]=a[2*l];
a[2*j+1]=a[2*l+1];
a[2*l]=t1;
a[2*l+1]=t2;
}
k=n1;
while (kj)
{
j=j-k;
k=k/2;
}
j=j+k;
}
for(i=1;i=s;i++)
{
u1=1;
u2=0;
m=(1i);
k=m1;
w1=fcos[i-1];
w2=-fsin[i-1];
for(j=1;j=k;j++)
{
for(l=j;lnn;l=l+m)
{
l1=l+k;
t1=a[2*l1]*u1-a[2*l1+1]*u2;
t2=a[2*l1]*u2+a[2*l1+1]*u1;
a[2*l1]=a[2*l]-t1;
a[2*l1+1]=a[2*l+1]-t2;
a[2*l]=a[2*l]+t1;
a[2*l+1]=a[2*l+1]+t2;
}
z=u1*w1-u2*w2;
u2=u1*w2+u2*w1;
u1=z;
}
}
for(i=1;i=nn/2;i++)
{
ar[i]=a[2*i+2]/nn;
ai[i]=-a[2*i+3]/nn;
a[i]=4*sqrt(ar[i]*ar[i]+ai[i]*ai[i]); // 幅值
}
}