本文目录一览：
- 1、求快速傅里叶变换的算法实现，C/C++/JAVA都行，要求适用于所有的整数，谢谢！
- 2、z7 7020中如何实现fft
- 3、java绘制wav波形
- 4、一维复数序列的快速傅里叶变换（FFT）
- 5、FFT的公式是什么和算法是怎样实现
求快速傅里叶变换的算法实现，C/C++/JAVA都行，要求适用于所有的整数，谢谢！
#includegraphics.h
#includestdio.h
#includemath.h
#includeconio.h
#includedos.h
float ar[2048],ai[2048];
float a[4100];
void testdata()
{
int i;
for(i=0;i512;i++)
{
ar[i]=50*cos(i*3.14159/32)+60*cos(i*3.14159/16)+53*cos(i*3.14159/64)+24*cos(i*3.14159/8)+10*cos(i*3.14159/128);
ai[i]=0;
}
}
void fft(int nn,int t)
{
int n1,n2,i,j,k,l,m,s,l1;
float t1,t2,x,y;
float w1,w2,u1,u2,z;
float pi=3.14159;
switch(nn)
{
case 2048: s=11; break;
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=cos(pi/k);
w2=t*sin(pi/k);
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]);
}
}
int main ()
{
int i;
int gdriver=DETECT,gmode;
initgraph(gdriver,gmode,"");
setfillstyle(SOLID_FILL,WHITE);
bar(0,0,639,479);
//模拟测试数据
testdata();
//波形显示
setcolor(BLACK);
line(10,119,550,119); // x轴
line(10,10,10,239); // y轴
setcolor(BLUE);
for(i=0;i511;i++)
line(i+10,119-ar[i],i+10+1,119-ar[i+1]);
//FFT
fft(512,-1);
//频谱显示
setcolor(BLACK);
line(10,459,550,459); // x轴
line(10,259,10,459); // y轴
setcolor(RED);
for(i=0;i256;i++)
line(2*i+10,459-a[i],2*i+10,459);
getch();
closegraph();
return 0;
}
z7 7020中如何实现fft
自己重新编程。
_ft是内建函数，不是matlab写的，看不到源代码的。下面是我写的一个fft,可以用functionxn=myfft(x)N=length(x);M=log2(N);xtmp=zeros(1,N);value=zeros(1,M);for i=0:N-1repr=i;fort=1:1:Mrepr=bitshift(i,1-t);value(t)=bitand(repr,1);endpos=0;for k=1:1:Mpos=pos+value(k)*2^(M-k);endxtmp(pos+1)=x(i+1);endfor i=1:deepth=2^(i-1);width=2^(M-i);for t=1:2^i:Nfor k=1:deepthtmp=xtmp(t+k-1);wn=width*(k-1);xtmp(t+k-1)=tmp+exp(-j*2*pi*wn/N)*xtmp(t+k+deepth-1);xtmp(t+k+deepth-1)=tmp-exp(-j*2*pi*wn/N)*xtmp(t+k+deepth-1);endendendxn=xtmp;
__FT相当于对行和列分别进行一维FFT运算。具体的实现办法如下： 先对各行逐一进行一维FFT，然后再对变换后的新矩阵的各列逐一进行一维FFT。
java绘制wav波形
本文讲述如何实时录音，以及将录音波形频谱实时显示的方法。Windows提供了一个多媒体控制接口(MCI)，用它可以录音，很方便。但是这种方法不能实时给出录音的原始数据，因此要显示波形和频谱都是不可能实现的。要达到实时的效果，就要使用Windows提供的另一套函数，即低级音频函数。其中和录音有关的是以waveIn开头的一组函数。
低级音频函数的使用比较繁琐，大致要有以下几个步骤。
用waveInOpen打开设备，并设置回调。因为要保证实时性，所以不能用查询的方式，而必须设置回调。
为设备分配足够的内存做缓冲区，动态分配或静态数组都可以。为了保证实时性，程序用了双缓冲技术，在处理一个缓冲区数据的同时另一个缓冲区用于录音。为了便于说明写成Buffer1、Buffer2。
将Buffer1关联到设备上去，waveInPrepareBuffer、waveInAddBuffer。
开始录音，waveInStart
当驱动程序填满这个缓冲区(Buffer1)时就会产生回调(消息为WIM_DATA)，这时立刻将Buffer2关联到设备上继续录音，然后处理Buffer1，当驱动程序填满Buffer2时又会产生回调，这是再将Buffer1关联到设备上，而去处理Buffer2，如此反复就使得录音能够实时的进行下去。
停止录音，waveInStop
关闭设备，waveInClose
上面就是实时录音的一个基本流程，产生回调时，其lParam就是一个指向WAVEHDR结构的指针，通过这个指针就可以得到原始数据(注意了：我这个程序只是单声道8位的PCM格式，因此正好一个字节就是一个采样点的值，对于其他格式的我就不知道了)。有了原始数据就可以很容易的画出波形了，而频谱也就是多一次FFT变换而已，没什么太难的。
这里是我用LCC写的例子，非常简单，而且相信大家都可以很容易的移植到BCB上来。另外，我偷了一下懒，FFT用了我以前封装到DLL中的函数(现在已经忘了怎么写了^_^)，可以在此下载 (不要觉得奇怪，就是这个时钟，将其中的DEMO.DLL拿来用就是了)。当然你也可以自己写一个FFT放到程序里去。
-*-*-PATCH-*-*-2002-06-11
上次那个程序还不好，原因是双缓冲还不能满足实时的要求，因此这次用了8缓冲循环技术。新的程序也是用LCC写的。
开始时依次将0－6号缓冲区关联到设备上，而7号缓冲区不关联，然后开始录音，这是驱动程序开始向0号缓冲区输出数据，当0号缓冲区数据满时产生回调，这时驱动程序会自动向1号缓冲区输出数据，但这时应将7号缓冲区关联到设备上去，然后处理0号缓冲区的数据。然后当1号缓冲区数据满时又会产生回调，这时驱动程序会自动向2号缓冲区输出数据，而将0号缓冲区关联到设备上去。这样如此反复，始终保持1个缓冲区用于驱动输出数据，1个缓冲区用于处理数据，6个缓冲区处于等待状态。
而如果只有两个缓冲区，就只能一个用于驱动输出，一个用于处理数据。这样在产生回调的时候就会出现没有处于等待状态的缓冲区的情况，这样驱动程序就不能自动向缓冲区输出数据。而要等waveInPrepareBuffer、waveInAddBuffer调用之后，这样就会出现一个小小的间隙，而导致录音不连续。
一维复数序列的快速傅里叶变换（FFT）
设x（N）为N点有限长离散序列，代入式（8－3）、式（8－4），并令 其傅里叶变换（DFT）为
地球物理数据处理基础
反变换（IDFT）为
地球物理数据处理基础
两者的差异只在于W的指数符号不同，以及差一个常数1/N，因此下面我们只讨论DFT正变换式（8－5）的运算量，其反变换式（8－6）的运算是完全相同的。
一般来说，W是复数，因此，X（j）也是复数，对于式（8－5）的傅里叶变换（DFT），计算一个X（j）值需要N次复数乘法和N－1次复数加法。而X（j）一共有N个值（j＝0，1，…，N－1），所以完成整个DFT运算总共需要N2次复数乘法和N（N－1）次复数加法。
直接计算DFT，乘法次数和加法次数都是与N2成正比的，当N很大时，运算量会很大，例如，当N＝8时，DFT需64次复数乘法；而当N＝1024时，DFT所需乘法为1048576次，即一百多万次的复数乘法运算，对运算速度要求高。所以需要改进DFT的计算方法，以减少运算次数。
分析Wjk，表面上有N2个数值，由于其周期性，实际上仅有N个不同的值W0，W1，…，WN－1。对于N＝2m时，由于其对称性，只有N/2个不同的值W0，W1，…，
地球物理数据处理基础
因此可以把长序列的DFT分解为短序列DFT，而前面已经分析DFT与N2成正比，所以N越小越有利。同时，利用ab＋ac＝a（b＋c）结合律法则，可以将同一个Wr对应的系数x（k）相加后再乘以Wr，就能大大减少运算次数。这就是快速傅里叶变换（FFT）的算法思路。
下面，我们来分析N＝2m情况下的FFT算法。
1.N＝4的FFT算法
对于m＝2，N＝4，式（8－5）傅里叶变换为
地球物理数据处理基础
将式（8－7）写成矩阵形式
地球物理数据处理基础
为了便于分析，将上式中的j，k写成二进制形式，即
地球物理数据处理基础
代入式（8－7），得
地球物理数据处理基础
分析Wjk的周期性来减少乘法次数
地球物理数据处理基础
则 代回式（8－9），整理得
地球物理数据处理基础
上式可分层计算，先计算内层，再计算外层时就利用内层计算的结果，可避免重复计算。写成分层形式
地球物理数据处理基础
则X（j1 j0）＝X2（j1 j0）。
上式表明对于N＝4的FFT，利用Wr的周期关系可分为m＝2步计算。实际上，利用Wr的对称性，仍可以对式（8－11）进行简化计算。考虑到
地球物理数据处理基础
式（8－11）可以简化为
地球物理数据处理基础
令j＝j0；k＝k0，并把上式表示为十进制，得
地球物理数据处理基础
可以看到，完成上式N＝4的FFT计算（表8－1）需要N·（m－1）/2＝2次复数乘法和N·m＝8次复数加法，比N＝4的DFT算法的N2＝16次复数乘法和N·（N－1）＝12次复数加法要少得多。
表8－1 N＝4的FFT算法计算过程
注：W0＝1；W1＝－i。
［例1］求N＝4样本序列1，3，3，1的频谱（表8－2）。
表8－2 N＝4样本序列
2.N＝8的FFT算法
类似N＝4的情况，用二进制形式表示，有
地球物理数据处理基础
写成分层计算的形式：
地球物理数据处理基础
则X（j2 j1 j0）＝X3（j2 j1 j0）。
对式（8－14）的X1（k1 k0 j0）进行展开，有
地球物理数据处理基础
还原成十进制，并令k＝2k1＋k0，即k＝0，1，2，3，有
地球物理数据处理基础
用类似的方法对式（8－14）的X2（k0 j1 j0），X3（j2 j1 j0）进行展开，整理得
地球物理数据处理基础
用式（8－16）、式（8－17）逐次计算到X3（j）＝X（j）（j＝0，1，…，7），即完成N＝23＝8的FFT计算，其详细过程见表8－3。
表8－3 N＝8的FFT算法计算过程
注：对于正变换 对于反变换 所
［例2］求N＝8样本序列（表8－4）x（k）＝1，2，1，1，3，2，1，2的频谱。
表8－4 N＝8样本序列
3.任意N＝2m的FFT算法
列出N＝4，N＝8的FFT计算公式，进行对比
地球物理数据处理基础
观察式（8－18）、式（8－19），不难看出，遵循如下规律：
（1）等式左边的下标由1递增到m，可用q＝1，2，…，m代替，则等式右边为q－1；
（2）k的上限为奇数且随q的增大而减小，至q＝m时为0，所以其取值范围为k＝0，1，2，…，（2m－q－1）；
（3）j的上限为奇数且随q的增大而增大，且q＝1时为0，其取值范围为j＝0，1，2，…，（2q－1－1）；
（4）k的系数，在等式左边为2q，等式右边为2q－1（包括W的幂指数）；
（5）等式左边序号中的常数是2的乘方形式，且幂指数比下标q小1，即2q－1；等式右边m对式子序号中的常数都是定值2m－1。
归纳上述规则，写出对于任意正整数m，N＝2m的FFT算法如下：
由X0（p）＝x（p）（p＝0，1，…，N－1）开始：
（1）对q＝1，2，…，m，执行（2）～（3）步；
（2）对k＝0，1，2，…，（2m－q－1）及j＝0，1，2，…，（2q－1－1），执行
地球物理数据处理基础
（3）j，k循环结束；
（4）q循环结束；由Xm（p）（p＝0，1，…，N－1）输出原始序列x（p）的频谱X（p）。
在计算机上很容易实现上述FFT算法程序，仅需要三个复数数组，编程步骤如下：
（1）设置复数数组X1（N－1），X2（N－1）和 （数组下界都从0开始）；
（2）把样本序列x赋给X1，即X1（k）＝x（k）（k＝0，1，…，N－1）；
（3）计算W，即正变换 反变换
（4）q＝1，2，…，m，若q为偶数，执行（6），否则执行第（5）步；
（5）k＝0，1，2，…，（2m－q－1）和j＝0，1，2，…，（2q－1－1）循环，作
X2（2qk＋j）＝X1（2q－1k＋j）＋X1（2q－1k＋j＋2m－1）
X2（2qk＋j＋2q－1）＝［X1（2q－1k＋j）－X1（2q－1k＋j＋2m－1）］W（2q－1k）
至k，j循环结束；
（6）k＝0，1，2，…，（2m－q－1）和j＝0，1，2，…，（2q－1－1）循环，作
X1（2qk＋j）＝X2（2q－1k＋j）＋X2（2q－1k＋j＋2m－1）
X1（2qk＋j＋2q－1）＝［X2（2q－1k＋j）－X2（2q－1k＋j＋2m－1）］W（2q－1k）
至k，j循环结束；
（7）q循环结束，若m为偶数，输出X1（j），否则输出X2（j）（j＝0，1，…，N－1），即为所求。
FFT的公式是什么和算法是怎样实现
二维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算法的实现。下面讨论一维FFT的算法原理。
【1D-FFT的算法实现】
设序列h(n)长度为N，将其按下标的奇偶性分成两组，即he和ho序列，它们的长度都是N/2。这样，可以将h(n)的FFT计算公式改写如下 ：
(A)
由于
所以，(A)式可以改写成下面的形式：
按照FFT的定义，上面的式子实际上是：
其中，k的取值范围是 0~N-1。
我们注意到He(k)和Ho(k)是N/2点的DFT，其周期是N/2。因此，H(k)DFT的前N/2点和后N/2点都可以用He(k)和Ho(k)来表示
