📄 fft842.cpp
字号:
#include "stdafx.h"
#include <math.h>
#include <conio.h>
#include<stdio.h>
#include<stdlib.h>
#include<dos.h>
#include<math.h>
#define PI 4.0*atan(1.0)
#define PI2 8.0*atan(1.0)
#define P7 1.0/sqrt(2.0)
void R2TX(unsigned int NTHPO,float *X,float *Y);
void R4TX(unsigned int NTHPO,float *X,float *Y);
void R8TX(unsigned int NXTLT,unsigned int NTHPO,unsigned int LENGT,float *X,float *Y);
int FFT842(bool bFFFT,unsigned int N,float *X,float *Y);
void RockWavelet(float dt,float mf,int lb,float *wb);
void ArrayFold(int nx,int nb,float *x,float *b,float *s);
int main(int argc, char* argv[])
{
int Nr=50;
int DT=4;/*ms*/
int Fdown=5;/*HZ*/
int Fup=65;/*HZ*/
int Vdown=10;/*中值下限*/
int Vup=39;/*中值上限*/
int Nw=1+(Fup-Fdown)/DT;/*子波长度*/
int i;
if(!(Nw%2))
Nw--;
int Ns=Nr+Nw-1;
float Wbuf[500],Wbuf1[500];
float Rbuf[500],Rbuf1[500];
float Sbuf[1000],Sbuf1[1000];
for(i=0;i<Ns;i++)
Wbuf[i]=Wbuf1[i]=Rbuf[i]=Rbuf1[i]=0.0;
/*求Rock子波*/
RockWavelet((float)DT,(float)(Vup-Vdown),Nw,Wbuf);
printf("输出Rock子波数据 Wbuf:\n");
printf("-------------------------------------------------------\n");
for(i=0;i<Ns;i++)
{
printf("\n%4d %9f %9f",i,Wbuf[i],Wbuf1[i]);
}
printf("\n");
/*反射系统*/
Rbuf[Nr/5]=(float)0.35;
Rbuf[Nr/3]=(float)0.67;
Rbuf[Nr/2]=(float)0.20;
Rbuf[Nr-1-Nr/3]=(float)-0.43;
Rbuf[Nr-1-Nr/4]=(float)-0.23;
Rbuf[Nr-1-Nr/5]=(float)0.89;
printf("输出反射系统数据 Rbuf:\n");
printf("-------------------------------------------------------\n");
for(i=0;i<Ns;i++)
{
printf("\n%4d %9f %9f",i,Rbuf[i],Rbuf1[i]);
}
printf("\n");
getch();
/*输出时域卷积记录*/
for(i=0;i<Ns;i++)
Sbuf[i]=0;
ArrayFold(Nr,Nw,Rbuf,Wbuf,Sbuf);
printf("输出时域卷积记录：\n");
printf("\n-------------------------------------------------------\n");
for(i=0;i<Ns;i++)
{
printf("\n%4d %9f",i,Sbuf[i]);
}
printf("\n");
getch();
/*对Rock子波进行FFT变换*/
FFT842(true,Ns,Wbuf,Wbuf1);
printf("After FFT Wbuf 输出Rock子波FFT变换后的结果:\n");
printf("-------------------------------------------------------\n");
for(i=0;i<Ns;i++)
{
printf("\n%4d %9f %9f",i,Wbuf[i],Wbuf1[i]);
}
printf("\n");
/*对反射系统进行FFT变换*/
FFT842(true,Ns,Rbuf,Rbuf1);
printf("After FFT Rbuf 输出反射系统FFT变换后的结果:\n");
printf("-------------------------------------------------------\n");
for(i=0;i<Ns;i++)
{
printf("\n%4d %9f %9f",i,Rbuf[i],Rbuf1[i]);
}
printf("\n");
getch();
/*对FFT变换后的信号进行复数乘法运算*/
/*Y(K)=X(K)*Y(k)*/
for(i=0;i<Ns;i++)
{
Sbuf[i]=Wbuf[i]*Rbuf[i]-Wbuf1[i]*Rbuf1[i];
Sbuf1[i]=Wbuf[i]*Rbuf1[i]+Wbuf1[i]*Rbuf[i];
}
printf("\n");
printf(" 输出两波信号傅立叶变换后复乘结果 Sbuf:\n");
printf("\n-------------------------------------------------------\n");
for(i=0;i<Ns;i++)
{
printf("\n%4d %9f %9f",i,Sbuf[i],Sbuf1[i]);
}
printf("\n");
/*对乘积结果进行IFFT变换*/
/*IFFT for Y(K)*/
FFT842(false,Ns,Sbuf,Sbuf1);
/*Output y(n)*/
printf("\n");
printf("After iFFT Sbuf 输出两波信号的傅立叶变换乘积的逆变换结果:\n");
printf("\n-------------------------------------------------------\n");
for(i=0;i<Ns;i++)
{
printf("\n%4d %9f %9f",i,Sbuf[i],Sbuf1[i]);
}
printf("\n");
getch();
return(0);
}
/*在时间域内求卷积*/
void ArrayFold(int nx,int nb,float *x,float *b,float *s)
{
int i,j;
for(i=0;i<nx;i++)
for(j=0;j<nb;j++)
{
s[i+j]=s[i+j]+(float)((double)x[i]*(double)b[j]); /* 卷积 */
}
}
/*--产生Rock子波---*/
void RockWavelet(float dt,float mf,int lb,float *wb)
{
int l1,ll1,ll2,ll3,mlb;
double tau,uu,wtt,wtt1,wtt2,wmax;
for(l1=0;l1<lb;l1++)
wb[l1]=0;
mlb=(lb-1)/2;
wmax=sqrt(3.1415927)/4.0;
for(ll1=0;ll1<mlb;ll1++)
{
tau=(double)(ll1+1)*(double)dt/1000.0;
uu=2.0*sqrt(6.0)*tau*(double)mf;
wtt=uu*uu/4.0;
wtt1=(wtt-0.5)*sqrt(3.1415927)/2.0;
wtt2=exp(-wtt);
ll2=mlb+1+ll1;
ll3=mlb-1-ll1;
wb[ll2]=-(float)((wtt1*wtt2)/wmax);
wb[ll3]=wb[ll2];
}
wb[mlb]=1;
}
/*************************************************************************
*If bFFFT=TRUE is FFT 正变换 *
*If bFFFT=FALSE is IFFT 逆变换 *
*N is the length Of a^N *
*X is Real 实部 *
*Y is Image 虚部 *
*If return 0 NORMAL *
*If retuen -1 is not a^N *
*************************************************************************/
int FFT842(bool bFFFT,unsigned int N,float *X,float *Y)
{
int M;
unsigned int L[15];
unsigned int N2POW,N8POW,LENGT,NXTLT,NTHPO;
unsigned int I,J,J1,J2,J3,J4,J5,J6,J7,J8,J9,J10,J11,J12,J13,J14;
float R;
for(J=1; J<=15; J++)
{
N2POW = J;
NTHPO = 1<<J;
if(NTHPO == N)
break;
}
if(NTHPO != N)
return(-1);
if(bFFFT)
{
for(J=0; J<NTHPO; J++)
Y[J] = -Y[J];
}
N8POW=N2POW/3;
if(N8POW)
{
for(J=1; J<=N8POW; J++)
{
NXTLT = 1<<(N2POW-3*J);
LENGT = 8*NXTLT;
R8TX(NXTLT,NTHPO,LENGT,X,Y);
}
}
M = N2POW-3*N8POW-1;
if(M == 0)
R2TX(NTHPO,X,Y);
else
{
if(M > 0)
R4TX(NTHPO,X,Y);
}
for(J=0;J<15;J++)
{
L[J] = 1;
M = N2POW-J;
if(M >= 1)
L[J] = 1<<M;
}
//
I=1;
for(J1=1;J1<=L[14];J1++)
for(J2=J1;J2<=L[13];J2+=L[14])
for(J3=J2;J3<=L[12];J3+=L[13])
for(J4=J3;J4<=L[11];J4+=L[12])
for(J5=J4;J5<=L[10];J5+=L[11])
for(J6=J5;J6<=L[9];J6+=L[10])
for(J7=J6;J7<=L[8];J7+=L[9])
for(J8=J7;J8<=L[7];J8+=L[8])
for(J9=J8;J9<=L[6];J9+=L[7])
for(J10=J9;J10<=L[5];J10+=L[6])
for(J11=J10;J11<=L[4];J11+=L[5])
for(J12=J11;J12<=L[3];J12+=L[4])
for(J13=J12;J13<=L[2];J13+=L[3])
for(J14=J13;J14<=L[1];J14+=L[2])
for(J=J14;J<=L[0];J+=L[1],I++)
{
M = I-J;
if(M < 0)
{
R=X[I-1];
X[I-1]=X[J-1];
X[J-1]=R;
R=Y[I-1];
Y[I-1]=Y[J-1];
Y[J-1]=R;
}
}
if(bFFFT)
{
for(J=0; J<NTHPO; J++)
Y[J] = -Y[J];
}
else
{
for(J=0; J<NTHPO; J++)
{
X[J] = (float)((double)X[J]/(double)NTHPO);
Y[J] = (float)((double)Y[J]/(double)NTHPO);
}
}
return(0);
}
void R2TX(unsigned int NTHPO,float *X,float *Y)
{
unsigned int i,j;
double R;
for(i=0;i<NTHPO;i+=2)
{
j = i+1;
R = (double)X[i]+(double)X[j];
X[j] = (float)((double)X[i]-(double)X[j]);
X[i] = (float)R;
R = (double)Y[i]+(double)Y[j];
Y[j] = (float)((double)Y[i]-(double)Y[j]);
Y[i] = (float)R;
}
}
//
void R4TX(unsigned int NTHPO,float *X,float *Y)
{
unsigned int i,j,k,l;
double R1,R2,R3,R4;
double FI1,FI2,FI3,FI4;
for(i=0;i<NTHPO;i+=4)
{
j = i+1;
k = j+1;
l = k+1;
R1 = (double)X[i]+(double)X[k];
R2 = (double)X[i]-(double)X[k];
R3 = (double)X[j]+(double)X[l];
R4 = (double)X[j]-(double)X[l];
FI1 = (double)Y[i]+(double)Y[k];
FI2 = (double)Y[i]-(double)Y[k];
FI3 = (double)Y[j]+(double)Y[l];
FI4 = (double)Y[j]-(double)Y[l];
X[i] = (float)(R1+R3);
Y[i] = (float)(FI1+FI3);
X[j] = (float)(R1-R3);
Y[j] = (float)(FI1-FI3);
X[k] = (float)(R2-FI4);
Y[k] = (float)(FI2+R4);
X[l] = (float)(R2+FI4);
Y[l] = (float)(FI2-R4);
}
}
//
void R8TX(unsigned int NXTLT,unsigned int NTHPO,unsigned int LENGT,float *X,float *Y)
{
unsigned int j,k0,k1,k2,k3,k4,k5,k6,k7;
double ar0,ar1,ar2,ar3,ar4,ar5,ar6,ar7;
double ai0,ai1,ai2,ai3,ai4,ai5,ai6,ai7;
double br0,br1,br2,br3,br4,br5,br6,br7;
double bi0,bi1,bi2,bi3,bi4,bi5,bi6,bi7;
double scale,arg,tr,ti;
double c1,c2,c3,c4,c5,c6,c7;
double s1,s2,s3,s4,s5,s6,s7;
scale = PI2/(double)LENGT;
for(j=0; j<NXTLT; j++)
{
arg = (double)j*scale;
c1 = cos(arg);
s1 = sin(arg);
c2 = c1*c1-s1*s1;
s2 = c1*s1+c1*s1;
c3 = c1*c2-s1*s2;
s3 = c2*s1+s2*c1;
c4 = c2*c2-s2*s2;
s4 = c2*s2+c2*s2;
c5 = c2*c3-s2*s3;
s5 = c3*s2+s3*c2;
c6 = c3*c3-s3*s3;
s6 = c3*s3+c3*s3;
c7 = c3*c4-s3*s4;
s7 = c4*s3+s4*c3;
for(k0=j; k0<NTHPO; k0+=LENGT)
{
k1 = k0+NXTLT;
k2 = k1+NXTLT;
k3 = k2+NXTLT;
k4 = k3+NXTLT;
k5 = k4+NXTLT;
k6 = k5+NXTLT;
k7 = k6+NXTLT;
//
ar0 = (double)X[k0]+(double)X[k4];
ar1 = (double)X[k1]+(double)X[k5];
ar2 = (double)X[k2]+(double)X[k6];
ar3 = (double)X[k3]+(double)X[k7];
ar4 = (double)X[k0]-(double)X[k4];
ar5 = (double)X[k1]-(double)X[k5];
ar6 = (double)X[k2]-(double)X[k6];
ar7 = (double)X[k3]-(double)X[k7];
//
ai0 = (double)Y[k0]+(double)Y[k4];
ai1 = (double)Y[k1]+(double)Y[k5];
ai2 = (double)Y[k2]+(double)Y[k6];
ai3 = (double)Y[k3]+(double)Y[k7];
ai4 = (double)Y[k0]-(double)Y[k4];
ai5 = (double)Y[k1]-(double)Y[k5];
ai6 = (double)Y[k2]-(double)Y[k6];
ai7 = (double)Y[k3]-(double)Y[k7];
//
br0 = ar0+ar2;
br1 = ar1+ar3;
br2 = ar0-ar2;
br3 = ar1-ar3;
br4 = ar4-ai6;
br5 = ar5-ai7;
br6 = ar4+ai6;
br7 = ar5+ai7;
bi0 = ai0+ai2;
bi1 = ai1+ai3;
bi2 = ai0-ai2;
bi3 = ai1-ai3;
bi4 = ai4+ar6;
bi5 = ai5+ar7;
bi6 = ai4-ar6;
bi7 = ai5-ar7;
//
X[k0] = (float)(br0+br1);
Y[k0] = (float)(bi0+bi1);
if(j)
{
X[k1]=(float)(c4*(br0-br1)-s4*(bi0-bi1));
Y[k1]=(float)(c4*(bi0-bi1)+s4*(br0-br1));
X[k2]=(float)(c2*(br2-bi3)-s2*(bi2+br3));
Y[k2]=(float)(c2*(bi2+br3)+s2*(br2-bi3));
X[k3]=(float)(c6*(br2+bi3)-s6*(bi2-br3));
Y[k3]=(float)(c6*(bi2-br3)+s6*(br2+bi3));
tr=P7*(br5-bi5);
ti=P7*(br5+bi5);
X[k4]=(float)(c1*(br4+tr)-s1*(bi4+ti));
Y[k4]=(float)(c1*(bi4+ti)+s1*(br4+tr));
X[k5]=(float)(c5*(br4-tr)-s5*(bi4-ti));
Y[k5]=(float)(c5*(bi4-ti)+s5*(br4-tr));
tr=0.0-P7*(br7+bi7);
ti=P7*(br7-bi7);
X[k6]=(float)(c3*(br6+tr)-s3*(bi6+ti));
Y[k6]=(float)(c3*(bi6+ti)+s3*(br6+tr));
X[k7]=(float)(c7*(br6-tr)-s7*(bi6-ti));
Y[k7]=(float)(c7*(bi6-ti)+s7*(br6-tr));
}
else
{
X[k1]=(float)(br0-br1);
Y[k1]=(float)(bi0-bi1);
X[k2]=(float)(br2-bi3);
Y[k2]=(float)(bi2+br3);
X[k3]=(float)(br2+bi3);
Y[k3]=(float)(bi2-br3);
tr=P7*(br5-bi5);
ti=P7*(br5+bi5);
X[k4]=(float)(br4+tr);
Y[k4]=(float)(bi4+ti);
X[k5]=(float)(br4-tr);
Y[k5]=(float)(bi4-ti);
tr=0.0-P7*(br7+bi7);
ti=P7*(br7-bi7);
X[k6]=(float)(br6+tr);
Y[k6]=(float)(bi6+ti);
X[k7]=(float)(br6-tr);
Y[k7]=(float)(bi6-ti);
}
}
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -