|
- //x--双精度实型一维数组,长度的n。开始时存在要变换数据的实部,最后
- //存放变换变换结果的实部
- //y--双精度实型一维数组,长度的n。开始时存在要变换数据的虚部,最后
- //存放变换变换结果的虚部
- //n--整型变量,数据长度如果是2的整数次幂,即n=(2的m次方);
- //数据长度如果不是2的整数次幂,即(2的m次方)<n<(2的m+1次方),
- //取前(2的m次方)个数据进行fft
- //sign--整型变量。sign=1时,子函数fft()计算离散傅立叶正变换(dft);
- //sign=-1时,子函数fft()计算离散傅立叶反变换(idft)
- void time2fft(float x[],float y[],int n,int sign)//时间抽样基2fft
- {
- int i,j,m,k,l,n1,n2,q;
- float ui1,ur1,c1,e,ui,ur,s1,t,tr,ti;
- q=n;
- for(j=1,i=1;i<n;i++)//取m值
- {
- m=i;
- j=2*j;
- if(j<n&&2*j>n)
- {
- n=j;
- break;
- }
- if(j==n)
- break;
- }
- n1=n-1;
- for(j=0,i=0;i<n1;i++)//变换
- {
- if (i<j)
- {
- tr=x[j];
- ti=y[j];
- x[j]=x[i];
- y[j]=y[i];
- x[i]=tr;
- y[i]=ti;
- }
- k=n/2;
- while(k<=j)
- {
- j=j-k;
- k=k/2;
- }
- j=j+k;
- }
- n1=1;
- for(l=1;l<=m;l++)
- {
- n1=2*n1;
- n2=n1/2;
-
- ur=1;
- ui=0;
- c1=cos(3.1415926/n2);
- s1=-sign*sin(3.1415926/n2);
- for(j=0;j<n2;j++)
- {
- for(i=j;i<n;i=i+n1)
- {
- k=i+n2;
- tr=ur*x[k]-ui*y[k];
- ti=ur*y[k]+ui*x[k];
- x[k]=x[i]-tr;
- y[k]=y[i]-ti;
- x[i]=x[i]+tr;
- y[i]=y[i]+ti;
- }
-
- ur1=ur*c1-ui*s1;
- ui1=ur*s1+ui*c1;
- ur=ur1;
- ui=ui1;
- }
- }
-
- if(sign==-1)
- {
- for(i=0;i<n;i++)
- {
- x[i]/=n;
- y[i]/=n;
- }
- }
复制代码
自己以前花了不少时间编出来的,可以正确运行 |
|