您的位置:

包含java实现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)来表示