📄 nllbfdlg.cpp
字号:
j=1;
for(i=1;i<=(nv2+1);i++) //改进时,将nm1替换为nv2
{
if(i<j)
{
mdata[0]=x[j];
x[j]=x[i];
x[i]=mdata[0];
//**改进算法部分begin
if(j<(nv2+1))
{
mdatb[0]=x[n-j+1];
x[n-j+1]=x[n-i+1];
x[n-i+1]=mdatb[0];
}
//**改进算法部分end
}
k=nv2;
sign=1; //循环标识
while(sign=1)
{
if(k>=j)
{
j=j+k;
break;//跳出while()循环
}
else
{
j=j-k;
k=k/2;
continue;//继续while()循环
}
}
}
double PI=3.1415926;
int le,le1,s,t,ip;
for(i=1;i<=m;i++)
{
le=pow(2,i);
le1=le/2;
u[0].re=1.0;
u[0].im=0.0;
w[0].re=cos(PI/le1);
w[0].im=-1.0*sin(PI/le1);
for(s=1;s<=le1;s++) //J
{
for(t=s;t<=n;t+=le) //I
{
ip=t+le1;
mdata[0]=CMul(x[ip],u[0]);
x[ip]=CSub(x[t],mdata[0]);
x[t]=CAdd(x[t],mdata[0]);
}
u[0]=CMul(u[0],w[0]);
}
}
return;
}
//**************************************************//
//******* 一维基2 FFT 算法 ***********//
//名称:fft1d
//参数说明:
//x为CNumber型的时域指针
//n为变换点数,2的整数次方,如4,8,16,32...... //
//m为变换次数
//n==2**m
//flag:1时正变换,不为1时为反变换
//***************** //
//作者:柯丹
//时间:2007.11.15
//*************************************************//
void CNllbfDlg::fft1d(CNumber* x, int n, int m, int flag)
{
//n=2**m;
int nv2,nm1,i,j,k,sign;
CNumber *mdata,*mdatb,*u,*w;
mdata=new CNumber[1];
mdatb=new CNumber[1];
u=new CNumber[1];
w=new CNumber[1];
nv2=n/2;
nm1=n-1;
j=1;
for(i=1;i<=(nv2+1);i++) //改进时,将nm1替换为nv2
{
if(i<j)
{
mdata[0]=x[j];
x[j]=x[i];
x[i]=mdata[0];
//**改进算法部分begin
if(j<(nv2+1))
{
mdatb[0]=x[n-j+1];
x[n-j+1]=x[n-i+1];
x[n-i+1]=mdatb[0];
}
//**改进算法部分end
}
k=nv2;
sign=1; //循环标识
while(sign=1)
{
if(k>=j)
{
j=j+k;
break;//跳出while()循环
}
else
{
j=j-k;
k=k/2;
continue;//继续while()循环
}
}
}
double PI=3.1415926;
int le,le1,s,t,ip;
for(i=1;i<=m;i++)
{
le=pow(2,i);
le1=le/2;
u[0].re=1.0;
u[0].im=0.0;
w[0].re=cos(PI/le1);
w[0].im=-1.0*sin(PI/le1);
for(s=1;s<=le1;s++) //J
{
for(t=s;t<=n;t+=le) //I
{
ip=t+le1;
mdata[0]=CMul(x[ip],u[0]);
x[ip]=CSub(x[t],mdata[0]);
x[t]=CAdd(x[t],mdata[0]);
}
u[0]=CMul(u[0],w[0]);
}
}
//**FFT逆变换,取共轭,再除以N // 测试成功
if(flag!=1)
{
for(i=1;i<=n;i++)
{
x[i].im=-1.0*x[i].im; //虚部取反,共轭
w[0].re=1.0/n;
w[0].im=0.0;
x[i]=CMul(x[i],w[0]);
}
}
return;
}
//******************************************************//
//******* 二维基2 FFT 算法 ***************//
//名称:fft2d
//参数说明:
//y为CNumber型的时域指针,存放二维数组
//m为二维数组的行数,为2的整数次方,如4,8,16,32...... //
//mr为行的变换次数, m=2**mr
//n为二维数组的列数
//nr为列的变换次数, n=2**nr
//flag:1时正变换,不为1时为反变换
//********************** //
//作者:柯丹
//时间:2007.11.16
//******************************************************//
void CNllbfDlg::fft2d(CNumber* y, int n, int m, int nr, int mr, int flag)
{
//数组行与列方向的数据均从1开始计数:y[1].....y[n]
int i,j,k;
CNumber *a;//,*b;
a=new CNumber[n+1];
//i=1;
//按行进行变换,调用一维FFT
for(i=1;i<=m;i++)//按行进行变换 //调试成功
{
for(j=1;j<=n;j++)
{
a[j]=y[(i-1)*n+j];
}
fft1d(a,n,nr,flag);//一维FFT
for(j=1;j<=n;j++)
{
y[(i-1)*n+j]=a[j];
}
}
///* //按列进行变换,调用一维FFT //****调试成功
for(j=1;j<=n;j++)//按列进行变换 //
{
for(i=1;i<=m;i++)
{
a[i]=y[(i-1)*n+j];
}
fft1d(a,m,mr,flag);//一维FFT
for(i=1;i<=m;i++)
{
y[(i-1)*n+j]=a[i];
}
}
return;
}
//功率谱密度
void CNllbfDlg::specd(CNumber* y, int n, int m, int nr, int mr, float length,float wn1[],float spd1[],float c[])
{
///*
//计算功率谱密度
int k,lng;//窗口长度
float tdta,tdtb,tdtc,tdtd;//temperatory data a,b,c,d
float maxr,minr;
int maxk,mink,tnum2,tnum,i,j;
maxr=0.0;
minr=1000.0;
float wn[5000],spd[5000];//wavenumber and spectrum density
//求最大小最等半径//调试通过
for(i=1;i<=m;i++)
{
for(j=1;j<=n;j++)
{
tdta=y[(i-1)*n+j].re*y[(i-1)*n+j].re;
tdtb=y[(i-1)*n+j].im*y[(i-1)*n+j].im;
tdtc=tdta+tdtb;
tdtc=sqrt(tdtc)*length/n;
//求最大等效半径
if(tdtc>=maxr)
maxr=tdtc;
//求最小等效半径
if(tdtc<=minr)
minr=tdtc;
}
}
maxk=int(maxr);//最大等效半径取整
mink=int(minr)+1;//最小等效半径取整
lng=maxk-mink+1;//数组的长度
tnum2=0;
for(k=mink;k<=maxk;k++)//对每一个J
{
tnum=0;
tdta=0.0;
tdtb=0.0;
tdtc=0.0;
tdtd=0.0;
for(i=1;i<=m;i++)
{
for(j=1;j<=n;j++)
{
tdta=y[(i-1)*n+j].re*y[(i-1)*n+j].re;
tdtb=y[(i-1)*n+j].im*y[(i-1)*n+j].im;
tdtc=tdta+tdtb;
tdtc=sqrt(tdtc)*length/n;
//介于k与(k+1)之间的等效半径
if(tdtc>k)
if(tdtc<(k+1))
{
tdtd=tdtd+tdtc*tdtc;
tnum=tnum+1;
}
}
}
//wn[k-mink+1]=(length)/(2*3.1415926*k);
//wn[k-mink+1]=(2*3.1415926*length)/(k);
wn[k-mink+1]=(2*3.1415926*k)/(length);
wn[k-mink+1]=1.0/wn[k-mink+1];//取倒数
if(tdtd>0)
{
if(tnum>0)
{
tdtc=(length*length*tdtd)/(n*n);//*tnum);
tdtc=tdtc/(n*n);
spd[k-mink+1]=tdtc;
tnum2=tnum2+1;
}
else
spd[k-mink+1]=-1.0;
}
}//对每一个J
//**********************
//剔除平均功率谱密度为零的值//调试通过
tnum=0;
for(k=mink;k<=maxk;k++)//对每一个J
{
tdtd=spd[k-mink+1];
if(tdtd>0)//功率谱不为零
{
spd1[tnum+1]=tdtd;//log(tdtd);
wn1[tnum+1]=log(wn[k-mink+1]);
tnum=tnum+1;
}
}//剔除结束
//**********************
//ercheng(&wn[1],&spd[1],&c[0],lng);//y=c[1]*x+c[2]*b
ercheng(&wn1[1],&spd1[1],&c[0],tnum);//y=c[1]*x+c[2]*b
c[2]=4.0+c[1]/2.0;//分维值,据吴小羊
c[3]=mink;//最小半径波数的取整
c[4]=maxk;//最大半径波数的取整
c[5]=float(lng);//最大半径波数与最小半径波数的间隔
c[6]=float(tnum2);//不为零的功率谱的个数
c[7]=float(tnum);//不为零的功率谱的个数
return;
}
///****数据扩充算法:将数据扩充为2的整数次幂
//名称:nexpand
//参数说明:
//num为int型的指针，存放数组
//num[0]存放原始数据的行数，横向 //
//num[1]存放原始数据的列数，纵向 //
//num[2]存放扩充后2的次幂数
//num[3]存放数据扩充后2的次幂数
//num[4]存放数据扩充后的行数，横向
//num[5]存放数据扩充后的列数，纵向
//********************** //
//作者:柯丹
//时间:2007.11.16
//******************************************************/
void CNllbfDlg::nexpand(int* num)
{
int i,se,msign,nsign;
msign=1;
nsign=1;
for(i=1;i<=11;i++)
{
se=pow(2,i);
while(msign)
{
if(num[0]<=se)
{
num[2]=i;
num[4]=pow(2,num[2]);
msign=0;
break;
}
break;
}
while(nsign)
{
if(num[1]<=se)
{
num[3]=i;
num[5]=pow(2,num[3]);
nsign=0;
break;
}
break;
}
}
return;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -