求c 或c++ 实现fft变换求频谱

发布网友 发布时间：2023-04-13 08:08

我来回答

共1个回答

热心网友 时间：2023-10-31 06:53

#include<stdio.h>
#include<math.h>
#include<string.h>

void kkfft(double pr[], double pi[], int n, int k, double fr[], double fi[], int l, int il)
{
int it,m,is,i,j,nv,l0;
double p,q,s,vr,vi,poddr,poddi;

for (it=0; it<=n-1; it++)
{
m = it;
is = 0;
for(i=0; i<=k-1; i++)
{
j = m/2;
is = 2*is+(m-2*j);
m = j;
}
fr[it] = pr[is];
fi[it] = pi[is];
}
/*/----------------------------*/
pr[0] = 1.0;
pi[0] = 0.0;
p = 6.283185306/(1.0*n);
pr[1] = cos(p);
pi[1] = -sin(p);

if (l!=0)
pi[1]=-pi[1];

for (i=2; i<=n-1; i++)
{
p = pr[i-1]*pr[1];
q = pi[i-1]*pi[1];
s = (pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
pr[i] = p-q;
pi[i] = s-p-q;
}

for (it=0; it<=n-2; it=it+2)
{
vr = fr[it];
vi = fi[it];
fr[it] = vr+fr[it+1];
fi[it] = vi+fi[it+1];
fr[it+1] = vr-fr[it+1];
fi[it+1] = vi-fi[it+1];
}
m = n/2;
nv = 2;

for (l0=k-2; l0>=0; l0--)
{
m = m/2;
nv = 2*nv;
for(it=0; it<=(m-1)*nv; it=it+nv)
for (j=0; j<=(nv/2)-1; j++)
{
p = pr[m*j]*fr[it+j+nv/2];
q = pi[m*j]*fi[it+j+nv/2];
s = pr[m*j]+pi[m*j];
s = s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
poddr = p-q;
poddi = s-p-q;
fr[it+j+nv/2] = fr[it+j]-poddr;
fi[it+j+nv/2] = fi[it+j]-poddi;
fr[it+j] = fr[it+j]+poddr;
fi[it+j] = fi[it+j]+poddi;
}
}

if(l!=0)
for(i=0; i<=n-1; i++)
{
fr[i] = fr[i]/(1.0*n);
fi[i] = fi[i]/(1.0*n);
}

if(il!=0)
for(i=0; i<=n-1; i++)
{
pr[i] = sqrt(fr[i]*fr[i]+fi[i]*fi[i]);
if(fabs(fr[i])<0.000001*fabs(fi[i]))
{
if ((fi[i]*fr[i])>0)
pi[i] = 90.0;
else
pi[i] = -90.0;
}
else
pi[i] = atan(fi[i]/fr[i])*360.0/6.283185306;
}

return;
}

void main()
{
FILE *fp;
int ii;
if((fp=fopen("aaa.wav","rb"))==NULL)
{ printf("Cannot open the file\n");exit(0);}//如果文件打不开，则关闭所有文件，显示出错
//--------------------------
double td[256],tm[256]={0}; //原文件虚部为0
double ffr[256],ffi[256] ;
fread(td,8,256,fp); //试着从fp所指文件中读取256个double型数据，读取一次，寸入td数组中

kkfft(td,tm,256,8,ffr,ffi,0,1);
fclose(fp);
for(ii=0;ii<256;ii++)
printf("%f ",ffr[ii]);
getchar();
}