加载中…胡广书第四章第9题第二问splfft
(2022-11-22 19:56:46) | 分类： VC |
代码：
complex xt;
float
es,e,a,a3,cc1,ss1,cc3,ss3,r1,r2,s1,s2,s3;
int m,n2,k,n4,j,is,id,i1,i2,i3,i0,n1,i,nn;
for(m=1;m<=16;m++)
{
nn=pow(2,m);
if(n==nn)break;
}
{
printf(" N is not a power of 2 ! \n");
return;
}
n2=n*2;
es=-isign*atan(1.0)*8.0;
for(k=1;k
{
n2=n2/2;
n4=n2/4;
e=es/n2;
a=0.0;
for(j=0;j
{
a3=3*a;
cc1=cos(a);
ss1=sin(a);
cc3=cos(a3);
ss3=sin(a3);
a=(j+1)*e;
is=j;
id=2*n2;
do
{
for(i0=is;i0
{
i1=i0+n4;
i2=i1+n4;
i3=i2+n4;
r1=x[i0].real-x[i2].real;
s1=x[i0].imag-x[i2].imag;
r2=x[i1].real-x[i3].real;
s2=x[i1].imag-x[i3].imag;
x[i0].real+=x[i2].real;x[i0].imag+=x[i2].imag;
x[i1].real+=x[i3].real;x[i1].imag+=x[i3].imag;
if(isign!=1)
{
s3=r1-s2;
r1+=s2;
s2=r2-s1;
r2+=s1;
}
else
{
s3=r1+s2;
r1=r1-s2;
s2=-r2-s1;
r2=-r2+s1;
}
x[i2].real=r1*cc1-s2*ss1;
x[i2].imag=-s2*cc1-r1*ss1;
x[i3].real=s3*cc3+r2*ss3;
x[i3].imag=r2*cc3-s3*ss3;
}
is=2*id-n2+j;
id=4*id;
}while(is
}
}
is=0;
id=4;
do
{
for(i0=is;i0
{i1=i0+1;
xt.real=x[i0].real;
xt.imag=x[i0].imag;
x[i0].real=xt.real+x[i1].real;
x[i0].imag=xt.imag+x[i1].imag;
x[i1].real=xt.real-x[i1].real;
x[i1].imag=xt.imag-x[i1].imag;
}
is=2*id-2;
id=4*id;
}while(is
j=1;
n1=n-1;
for(i=1;i<=n1;i++)
{
if(i
{
xt.real=x[j-1].real;
xt.imag=x[j-1].imag;
x[j-1].real=x[i-1].real;
x[j-1].imag=x[i-1].imag;
x[i-1].real=xt.real;
x[i-1].imag=xt.imag;
}
k=n/2;
do
{
if(k>=j)
break;
j-=k;
k/=2;
}while(1);
j+=k;
}
if(isign==-1)
return;
for(i=0;i
{
x[i].real/=(float)n;
x[i].imag/=(float)(n);
}
}
#include "stdafx.h"
#include "MSP.h"
#include "stdio.h"
#include
#include
#define N 4096
#define PI
3.14159265
complex x[N];
void msplfft(complex x[],int n,int isign)
{
if(m>16)
int main(int argc, char* argv[])
{
FILE *fp,*fp1;
int i=0;
if((fp=fopen("a.txt","r"))==NULL)
{
printf("cannot open this file\n");
exit(0);
}
for(i=0;i
{
fscanf(fp,"%f\n",&x[i].real);
printf("x[%d].real=:%f\n",i,x[i].real);
x[i].imag=0;
}
msplfft(x,N,-1);
fclose(fp);
if((fp1=fopen("d.txt","w"))==NULL)
{
printf("cannot open this file\n");
exit(0);
}
for(i=0;i
{
fprintf(fp1,"%f\n",sqrt(x[i].real*x[i].real+x[i].imag*x[i].imag));
}
fclose(fp1);
return 0;
}
幅度谱：
喜欢
0
赠金笔