#include<math.h>
#include<stdio.h>
void four1(float data[],unsigned long nn,int isign)
{
unsigned long n,mmax,m,j,istep,i;
double wtemp,wr,wpr,wpi,wi,theta;
float tempr,tempi;
n=nn<<1;
j=1;
for(i=1;i<n;i+=2)
{
if(j>i)
{
tempr=data[j],data[j]=data[i],data[i]=tempr;
tempr=data[j+1],data[j+1]=data[i+1];data[i+1]=tempr;
}
m=n>>1;
while(m>=2 && j>m)
{
j-=m;
m>>=1;
}
j+=m;
}
mmax=2;
while(n>mmax)
{
istep=mmax<<1;
theta=isign*(6.28318530717959/mmax);
wtemp=sin(0.5*theta);
wpr=-2.0*wtemp*wtemp;
wpi=sin(theta);
wr=1.0;
wi=0.0;
for(m=1;m<mmax;m+=2){
for(i=m;i<=n;i+=istep){
j=i+mmax;
tempr=wr*data[j]-wi*data[j+1];
tempi=wr*data[j+1]+wi*data[j];
data[j]=data[i]-tempr;
data[j+1]=data[i+1]-tempi;
data[i]+=tempr;
data[i+1]+=tempi;
}
wr=(wtemp=wr)*wpr-wi*wpi+wr;
wi=wi*wpr+wtemp*wpi+wi;
}
mmax=istep;
}
}
void realft(float data[],unsigned long n,int isign)
{
void four1(float data[],unsigned long nn,int isign);
unsigned long i,i1,i2,i3,i4,np3;
float c1=0.5,c2,h1r,h1i,h2r,h2i;
double wr,wi,wpr,wpi,wtemp,theta;
theta=3.141592653589793/(double)(n>>1);
if(isign == 1)
{
c2=-0.5;
four1(data,n>>1,1);
}
else
{
c2=0.5;
theta=-theta;
}
wtemp=sin(0.5*theta);
wpr=-2.0*wtemp*wtemp;
wpi=sin(theta);
wr=1.0+wpr;
wi=wpi;
np3=n+3;
for(i=2;i<=(n>>2);i++)
{
i4=1+(i3=np3-(i2=1+(i1=i+i-1)));
h1r=c1*(data[i1]+data[i3]);
h1i=c1*(data[i2]-data[i4]);
h2r=-c2*(data[i2]+data[i4]);
h2i=c2*(data[i1]-data[i3]);
data[i1]=h1r+wr*h2r-wi*h2i;
data[i2]=h1i+wr*h2i+wi*h2r;
data[i3]=h1r-wr*h2r+wi*h2i;
data[i4]=-h1i+wr*h2i+wi*h2r;
wr=(wtemp=wr)*wpr-wi*wpi+wr;
wi=wi*wpr+wtemp*wpi+wi;
}
if(isign==1)
{
data[1]=(h1r=data[1])+data[2];
data[2]=h1r-data[2];
}
else
{
data[1]=c1*(h1r=data[1]+data[2]);
data[2]=c1*(h1r-data[2]);
four1(data,n>>1,-1);
}
}
main()
{
unsigned long n,m,pi;
n=0;
m=0;
pi=3.141592653589793;
float i,a[512],b[256],c[512],d[512];
FILE *fpa,*fpar,*fpai,*fpb,*fpc,*fpcr,*fpci,*fpd;
if((fpa=fopen("a1.txt","w"))==NULL){ //打开文件
printf("File open error!\n");
exit(0);
}
if((fpar=fopen("a1real.txt","w"))==NULL){ //打开文件
printf("File open error!\n");
exit(0);
}
if((fpai=fopen("a1image.txt","w"))==NULL){ //打开文件
printf("File open error!\n");
exit(0);
}
if((fpb=fopen("b1.txt","w"))==NULL){ //打开文件
printf("File open error!\n");
exit(0);
}
if((fpc=fopen("c1.txt","w"))==NULL){ //打开文件
printf("File open error!\n");
exit(0);
}
if((fpcr=fopen("c1real.txt","w"))==NULL){ //打开文件
printf("File open error!\n");
exit(0);
}
if((fpci=fopen("c1image.txt","w"))==NULL){ //打开文件
printf("File open error!\n");
exit(0);
}
if((fpd=fopen("d1.txt","w"))==NULL){ //打开文件
printf("File open error!\n");
exit(0);
}
/* if((fpf=fopen("fff.txt","w"))==NULL){ //打开文件
printf("File open error!\n");
exit(0);
}
for(n=0;n<=600;n++)
{
f[n]=2.44141*n;
fprintf(fpf,"%f\n",f[n]);
} */
for(i=0;i<=0.2048;i+=0.0004) //信号采样
{
printf("aa\n");
a[n]=0.2*sin(2*pi*50*i)+0.16*sin(2*pi*150*i)+0.05*cos(2*pi*250*i);
fprintf(fpa,"%f\n",i);
// fprintf(fpa,"%c",'\x20');
// fprintf(fpa,"%f\n",a[n]); //空格
n++;
}
realft(a,512,1); //FFT变换，偶数项为实部，奇数相为虚部
for(n=0;n<512;n++)
{
a[n]=a[n]/256;
if((n%2)==0){
fprintf(fpar,"%f\n",a[n]); //实部
// fprintf(fpar,"%c",'\x20');
}
else
{
fprintf(fpai,"%f\n",a[n]); //虚部
// fprintf(fpai,"%c",'\x20');
}
}
for(n=0;n<256;n++)
{
printf("bb\n");
b[n]=sqrt(a[2*n]*a[2*n]+a[2*n+1]*a[2*n+1]); //功率谱=开根号（实部2+虚部2）
fprintf(fpb,"%f\n",b[n]);
//fprintf(fpb,"%c",'\x20');
}
for(i=0;i<=0.2048;i+=0.0004)
{
printf("cc\n");
c[m]=0.05*cos(2*pi*250*i);
fprintf(fpc,"%f\n",c[m]);
// fprintf(fpc,"%c",'\x20');
m++;
}
realft(c,512,1); //C数列的FFT变换
for(n=0;n<512;n++)
{
c[n]=c[n]/256;
if((n%2)==0){
fprintf(fpcr,"%f\n",c[n]); //实部
// fprintf(fpar,"%c",'\x20');
}
else
{
fprintf(fpci,"%f\n",c[n]); //虚部
// fprintf(fpai,"%c",'\x20');
}
}
for(n=0;n<512;n+=2)
{
d[n]=a[n]*c[n]+a[n+1]*c[n+1]; //C数列与A数列共轭相乘
d[n+1]=-a[n]*c[n+1]+a[n+1]*c[n];
}
// realft(d,512,-1); //FFT反变换
for(n=0;n<512;n++)
{
fprintf(fpd,"%f\n",d[n]);
// fprintf(fpd,"%c",'\x20');
}
}
评论0