fft算法的实现
Ⅰ 怎样用C语言实现FFT算法啊
1、二维FFT相当于对行和列分别进行一维FFT运算。具体的实现办法如下：
先对各行逐一进行一维FFT，然后再对变换后的新矩阵的各列逐一进行一维FFT。相应的伪代码如下所示：
for (int i=0; i<M; i++)
FFT_1D(ROW[i]，N);
for (int j=0; j<N; j++)
FFT_1D(COL[j]，M);
其中，ROW[i]表示矩阵的第i行。注意这只是一个简单的记法，并不能完全照抄。还需要通过一些语句来生成各行的数据。同理，COL[i]是对矩阵的第i列的一种简单表示方法。
所以，关键是一维FFT算法的实现。
2、例程：
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#defineN1000
/*定义复数类型*/
typedefstruct{
doublereal;
doubleimg;
}complex;
complexx[N],*W;/*输入序列,变换核*/
intsize_x=0;/*输入序列的大小，在本程序中仅限2的次幂*/
doublePI;/*圆周率*/
voidfft();/*快速傅里叶变换*/
voidinitW();/*初始化变换核*/
voidchange();/*变址*/
voidadd(complex,complex,complex*);/*复数加法*/
voidmul(complex,complex,complex*);/*复数乘法*/
voidsub(complex,complex,complex*);/*复数减法*/
voidoutput();
intmain(){
inti;/*输出结果*/
system("cls");
PI=atan(1)*4;
printf("Pleaseinputthesizeofx: ");
scanf("%d",&size_x);
printf("Pleaseinputthedatainx[N]: ");
for(i=0;i<size_x;i++)
scanf("%lf%lf",&x[i].real,&x[i].img);
initW();
fft();
output();
return0;
}
/*快速傅里叶变换*/
voidfft(){
inti=0,j=0,k=0,l=0;
complexup,down,proct;
change();
for(i=0;i<log(size_x)/log(2);i++){/*一级蝶形运算*/
l=1<<i;
for(j=0;j<size_x;j+=2*l){/*一组蝶形运算*/
for(k=0;k<l;k++){/*一个蝶形运算*/
mul(x[j+k+l],W[size_x*k/2/l],&proct);
add(x[j+k],proct,&up);
sub(x[j+k],proct,&down);
x[j+k]=up;
x[j+k+l]=down;
}
}
}
}
/*初始化变换核*/
voidinitW(){
inti;
W=(complex*)malloc(sizeof(complex)*size_x);
for(i=0;i<size_x;i++){
W[i].real=cos(2*PI/size_x*i);
W[i].img=-1*sin(2*PI/size_x*i);
}
}
/*变址计算，将x(n)码位倒置*/
voidchange(){
complextemp;
unsignedshorti=0,j=0,k=0;
doublet;
for(i=0;i<size_x;i++){
k=i;j=0;
t=(log(size_x)/log(2));
while((t--)>0){
j=j<<1;
j|=(k&1);
k=k>>1;
}
if(j>i){
temp=x[i];
x[i]=x[j];
x[j]=temp;
}
}
}
/*输出傅里叶变换的结果*/
voidoutput(){
inti;
printf("Theresultareasfollows ");
for(i=0;i<size_x;i++){
printf("%.4f",x[i].real);
if(x[i].img>=0.0001)printf("+%.4fj ",x[i].img);
elseif(fabs(x[i].img)<0.0001)printf(" ");
elseprintf("%.4fj ",x[i].img);
}
}
voidadd(complexa,complexb,complex*c){
c->real=a.real+b.real;
c->img=a.img+b.img;
}
voidmul(complexa,complexb,complex*c){
c->real=a.real*b.real-a.img*b.img;
c->img=a.real*b.img+a.img*b.real;
}
voidsub(complexa,complexb,complex*c){
c->real=a.real-b.real;
c->img=a.img-b.img;
}
Ⅱ matlab中fft（）函数是什么意思
FFT（快速傅里叶变换）是一种实现DFT（离散傅里叶变换）的快速算法，是利用复数形式的离散傅里叶变换来计算实数形式的离散傅里叶变换，matlab中的fft()函数是实现该算法的实现。
MATLAB它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。
快速傅里叶变换, 即利用计算机计算离散傅里叶变换（DFT)的高效、快速计算方法的统称，简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显着。
(2)fft算法的实现扩展阅读：
matlab优势特点：
1、高效的数值计算及符号计算功能，能使用户从繁杂的数学运算分析中解脱出来；
2、具有完备的图形处理功能，实现计算结果和编程的可视化；
3、友好的用户界面及接近数学表达式的自然化语言，使学者易于学习和掌握；
4、功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ，为用户提供了大量方便实用的处理工具。
参考资料来源：
网络-快速傅里叶变换
网络-MATLAB
Ⅲ 在DSP上实现FFT算法
void FFT( COMPLEX *Y, int N) /* input sample array, number of points */
{
COMPLEX temp1,temp2; /*temporary storage variables */蔽型
int i,j,k; /*loop counter variables */
int upper_leg, lower_leg; /*index of upper/lower butterfly leg */
int leg_diff; /宏握猜*difference between upper/lower leg */
int num_stages=0; /*number of FFT stages, or iterations */
int index, step; /*index and step between twiddle factor*/
/* log(base 2) of # of points = # of stages */
i=1;
do
{
num_stages+=1;
i = i *2 ;
} while (i!=N);
/* starting difference between upper and lower butterfly legs*/
leg_diff = N/2;
/* step between values in twiddle factor array twiddle.h */
step = 512 / N;
/* For N-point FFT */
for ( i = 0 ; i < num_stages ; i++ )
{
index = 0;
for ( j = 0; j < leg_diff ; j++ )
{
for ( upper_leg = j; upper_leg < N ; upper_leg += (2*leg_diff) )
{
lower_leg = upper_leg + leg_diff;
temp1.real=(Y[upper_leg]).real + (Y[lower_leg]).real;
temp1.imag=(Y[upper_leg]).imag + (Y[lower_leg]).imag;
temp2.real=(Y[upper_leg]).real - (Y[lower_leg]).real;
temp2.imag=(Y[upper_leg]).imag - (Y[lower_leg]).imag;
(Y[lower_leg]).real = ((long)temp2.real * (w[index]).real)/8192;
(Y[lower_leg]).real -= ((long)temp2.imag * (w[index]).imag)/8192;
(Y[lower_leg]).imag = ((long)temp2.real * (w[index]).imag)/8192;
(Y[lower_leg]).imag += ((long)temp2.imag * (w[index]).real)/8192;
(Y[upper_leg]).real = temp1.real;
(Y[upper_leg]).imag = temp1.imag;
}
index+=step;
}
leg_diff = leg_diff / 2;
step *= 2;
}
/* bit reversal for resequencing data */
j=0;
for ( i=1 ; i < (N-1) ; i++ )
{
k = N / 2;
while ( k <= j)
{
j = j - k;
k >>皮旅= 1;
}
j = j + k;
if ( i < j )
{
temp1.real = (Y[j]).real;
temp1.imag = (Y[j]).imag;
(Y[j]).real = (Y[i]).real;
(Y[j]).imag = (Y[i]).imag;
(Y[i]).real = temp1.real;
(Y[i]).imag = temp1.imag;
}
}
return;
}
参考一下的吧，这个是TI官方的在5416上实现的程序~
Ⅳ 实序列的FFT算法
在以上讨论FFT算法中，均假定序列x(l)为复的，但实际问题中的序列大多为实的。当然，我们可以把实序列处理成虚部为零的复序列。因此，就要引进许多零参加运算。这样一来，在机器运算时间和存储单元方面都将造成很大的浪费。在本段中，我们介绍对实序列x(l)应用FFT算法的一个有效方法。
1.同时计算两个实序列的FFT算法
设有N=4的两个实序列x1(l)与x2(l)。为了求得它们的谱X1(m)与X2(m)，我们用此二实序列构造成如下复序列
物探数字信号分析与处理技术
利用上一段的方法，可以求得复序列x(l)的谱X(m)。根据(7-3-1)得到
物探数字信号分析与处理技术
上式中的m用N-m代替，则得
物探数字信号分析与处理技术
将上式两端取共轭，根据对称性有
物探数字信号分析与处理技术
根据DFT的复共轭性质，对于实序列x1(l)与x2(l)，有
物探数字信号分析与处理技术
于是从(7-3-4)得到
物探数字信号分析与处理技术
联立求解(7-3-2)和(7-3-6)便得到
物探数字信号分析与处理技术
例如设有两个N=4点的实序列，
物探数字信号分析与处理技术
我们用它们构造一个N=4点的复序列
物探数字信号分析与处理技术
利用FFT算法求X(m)，m=0，1，2，3(图7-3-1)，
图7-3-1 N=4点的FFT算法流程图
于是得到
物探数字信号分析与处理技术
因此从式(7-3-7)得到
物探数字信号分析与处理技术
物探数字信号分析与处理技术
2.实序列的FFT算法
设有N点的实序列x(l)，l=0，1，2，…，N-1。按照点的奇偶编号，将它们分成N/2个点的两个子序列
物探数字信号分析与处理技术
设x1(l)的谱与x2(l)的谱分别为X1(m)与X2(m)
物探数字信号分析与处理技术
其中
于是可以将实序列x(l)的谱X(m)，用两个子序列x1(l)，x2(l)的谱X1(m)，X2(m)来表示
物探数字信号分析与处理技术
其中
物探数字信号分析与处理技术
注意，x1(l)，x2(l)与X1(m)，X2(m)均以N/2为周期，
利用x1(l)、x2(l)构成如下复序列
物探数字信号分析与处理技术
利用FFT算法可以求得复序列 的谱 。根据(7-3-7)就求得两个实子序列的谱X1(m)与X2(m)
物探数字信号分析与处理技术
有了X1(m)，X2(m)，根据(7-3-10)就可求得X(m)。以上就是用FFT算法求实序列x(l)的谱X(m)的方法。必须指出，用公式(7-3-10)求X(m)时，第一，两个实子序列的谱X1(m)，X2(m)及复序列x珓(l)的谱珘X(m)均是以N/2为周期的周期序列;第二，由于x
(l)是实序列，根据DFT的复共轭性质有X(m)=X*(N-m)，m=0，1，…，N/2，故只需求得前(N/2)+1个点的X(m)，就得到全部N个点的X(m)了
例如，有N=8点的实序列，
物探数字信号分析与处理技术
首先，按点的奇偶编号分成两个实子序列，
物探数字信号分析与处理技术
其次用它们构造如下复序列，
物探数字信号分析与处理技术
用FFT算法求此复序列的谱 (图7-3-2)
图7-3-2 N=4点的FFT算法流程图
于是得到:
根据周期性，有
物探数字信号分析与处理技术
根据(7-3-12)式，
物探数字信号分析与处理技术
根据周期性，有
物探数字信号分析与处理技术
故最终由(7-3-10)得到
物探数字信号分析与处理技术