[讨论]关于fft算法问题
本帖最后由 wdhd 于 2016-3-14 14:51 编辑傅立叶变换的程序贴出如下:
#include <stdio.h>
#include <math.h>
#define N 64
#define m 6 //2^m=N
/*
float twiddle={1.0,0.707,0.0,-0.707,};
float x_r={1,1,1,1,0,0,0,0,};
float x_i; //N=8
*/
float twiddle={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={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;
void fft_init()
{
int i;
for(i=0;i<N;i++)
x_i=0.0;
}
void bitrev()
{
int p=1,q,i;
int bit_rev;
float xx_r;
bit_rev=0;
while(p<N)
{
for(q=0;q<p;q++)
{
bit_rev=bit_rev*2;
bit_rev=bit_rev+1;
}
p=p*2;
}
for(i=0;i<N;i++)
{
xx_r=x_r;
}
for(i=0;i<N;i++)
{
x_r=xx_r];
}
}
void display()
{
int i;
for(i=0;i<N;i++)
printf("%f\t%f\n",x_r,x_i);
}
void fft1()
{
int L,i,b,j,p,k,tx1,tx2;
float TI,TR,temp;
float tw1,tw2;
for(L=1;L<=m;L++)
{ /* for(1) */
b=1; i=L-1;
while(i>0)
{b=b*2;i--;} /* b= 2^(L-1) */
for(j=0;j<=b-1;j++) /* for (2) */
{ p=1; i=m-L;
while(i>0) /* p=pow(2,7-L)*j; */
{p=p*2; i--;}
p=p*j;
tx1=p%(N);
tx2=tx1+(3*N)/4;
tx2=tx2%(N);
if (tx1>=(N/2))
tw1=-twiddle;
else
tw1=twiddle;
if (tx2>=(N/2))
tw2=-twiddle;
else
tw2=twiddle;
for(k=j;k<N;k=k+2*b) /* for (3) */
{TR=x_r; TI=x_i; temp=x_r;
x_r=x_r+x_r*tw1+x_i*tw2;
x_i=x_i-x_r*tw2+x_i*tw1;
x_r=TR-x_r*tw1-x_i*tw2;
x_i=TI+temp*tw2-x_i*tw1;
}
}
}
}
void dft()
{
int i,n,k,tx1,tx2;
float tw1,tw2;
float xx_r,xx_i;
//clear any data in Real and Imaginary result arrays prior to DFT
for(k=0;k<=N-1;k++)
xx_r=xx_i=x_i=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;
else
tw1=twiddle;
if (tx2>=(N/2))
tw2=-twiddle;
else
tw2=twiddle;
xx_r=xx_r+x_r*tw1;
xx_i=xx_i+x_r*tw2;
}
xx_i=-xx_i;
}
//display
for(i=0;i<N;i++)
printf("%f\t%f\n",xx_r,xx_i);
}
void main()
{
/*
fft_init();
bitrev();
fft1();
display(); //FFT
*/
dft(); //DFT
}
回复：（simon21）[讨论]关于fft算法问题
本帖最后由 wdhd 于 2016-3-14 14:51 编辑其中分别使用fft和dft做变换
源程序如上: 其中函数fft1()代表fft变换,函数dft()代表dft变换.
当取N=8的时候,fft和dft算出来的数据是一样的.
fft和dft计算结果如下:
4.000000 0.000000
1.000000 -2.414000
0.000000 0.000000
1.000000 -0.414000
0.000000 0.000000
1.000000 0.414000
0.000000 0.000000
1.000000 2.414000
Press any key to continue
但是,当N=64的时候,fft和dft算出的数据x_r[ i ],x_ i就不相同:
fft算出的结果如下:
32.000000 0.000000
0.968460 -20.365383
0.000000 0.000000
0.993620 -6.749197
0.000000 0.000000
3.875856 -1.388196
0.000000 0.000000
... ... ... ... ;这里只取i=0~6的x_r[ i ],x_ i显示
dft算出的结果如下:
32.000000 0.000000
2.663399 -18.705200
-0.000000 0.000001
2.663400 -8.407199
0.000000 0.000000
2.663400 -2.332000
0.000000 0.000000
... ... ... ... ;这里只取i=0~6的x_r[ i ],x_ i显示
书上说,dft和fft的计算结果应该是完全一样的.这是怎么会事?是不是程序中有问题,谢谢!
回复：（simon21）[讨论]关于fft算法问题
本帖最后由 wdhd 于 2016-8-31 15:29 编辑fft的源程序不是现成的吗?看看这个网站，上面有世界上最优秀的fft算法，而且是免费的，matlab的fft就是用的他的算法。
http://www.fftw.org/
回复：（Anonymous）回复：（simon21）[讨论]关于fft算...
在某些场合也不完全通用吧，需要自己更具情况写回复：（simon21）回复：（Anonymous）回复：（simon2...
本帖最后由 wdhd 于 2016-3-14 14:51 编辑以下是引用simon21在2005-9-23 9:00:07的发言：
在某些场合也不完全通用吧，需要自己更具情况写
有啥不能通用的地方说来听听，长长见识啦，偶没有自己编过，
//一向都觉得只是调用个子程序罢了
[此贴子已经被作者于2005-9-23 15:45:38编辑过]
如果只是用在信号处理里这些程序基本够了，不过我见过一篇文章将FFT和其他算法结合起来做数值分析的，用通用程序就不行了
回复：（simon21）[讨论]关于fft算法问题
个人感觉自己写一遍对FFT的理解程度能提高一个层次<BR>要用FFT的还是要自己写一下，不管最后用不用 我也正在编程中. 好主意！
页:
[1]