还是fft

lanren_me 2003-01-27 10:40:01
使用fft算法确定采集点中频谱的变化中,

是否振幅一样的波,分析出来的高度(值)是一样???

现在假设1024个模拟的点;连个频率同时输进:

double ff1= 2.0*PI*(f1)/F;//F--为采样频率;f1;f2为模拟的频率;
double ff2=2.0*PI*(f2)/F;
double ri[1024]; //fft数据源--模拟的采集点;
for (int i=0;i<1024;i++)
{
double k=sin(ff1*i);
k+=sin(ff2*i);
ri[i]=k;
}

这样对数据源进行fft;是否这两个频率确定的频谱高度是一样的???我现在的结果不是一样高,但是变换不是很大;

其中fft的代码如下:

//-----------------------------------------------------------------------------
void fft1(double *br,double *bi,
int il,const double *st,const double *ct)
{
//the size of br,bi should be equal to il.
//and size of st and ct no less than il.
int is=il>>1;
int mpx=0;
int ic=1;
int ia;

while(is)
{
is>>=1;
mpx++;
}

is=il;

for(ia=1;ia<=mpx;ia++)
{
int ib;
int ka=0;
is>>=1;
for(ib=1;ib<=ic;ib++)
{
int in=1;
int k;
for(k=1;k<=is;k++)
{
int j1=ka+k;
int j2=j1+is;
double xr=br[j1];
double xi=bi[j1];
double yr=br[j2];
double yi=bi[j2];
br[j1]=xr+yr;
bi[j1]=xi+yi;
xr-=yr;
xi-=yi;
br[j2]=xr*ct[in]-xi*st[in];
bi[j2]=xr*st[in]+xi*ct[in];
in+=ic;
}
ka+=is<<1;
}
ic<<=1;
}
}

void binrv(double *bc,int il,const int *lb)
{
int is=il-1;
int i;
for(i=2;i<=is;i++)
{
int ig=lb[i];
if(ig<=i)
continue;
double xx=bc[i];
bc[i]=bc[ig];
bc[ig]=xx;
}
}

void brtab(int *lbr,int il)
{
int is=il>>1;
int mpx=0;
while(is)
{
is>>=1;
mpx++;
}
int ln;
for(ln=1;ln<=il;ln++)
{
int j1=ln-1;
int ibord=0;
for(int k=1;k<=mpx;k++)
{
int j2=j1>>1;
ibord=ibord*2+(j1-2*j2);
j1=j2;
}
lbr[ln]=ibord+1;
}
}

void cstab(double *st,double *ct,int il,int ity)
{
double yy;
yy=-PI*2.0/il;
if(ity<0)
yy=-yy;
int l;
double ang=0.0;
for(l=1;l<=il;l++)
{
st[l]=sin(ang);
ct[l]=cos(ang);
ang+=yy;
}
}

void __fastcall FFT(double *br,double *bi,int n,int ity)
//ity>0 forward ,ity<0 backward
//n=pow(2,m);
//Here br and bi start from 0.
{
double *ct,*st;
int *il;
if(n<=0)
return;
double invsqrtn=1.0/sqrt(n);
ct=(double *)malloc(sizeof(double)*n)-1;
st=(double *)malloc(sizeof(double)*n)-1;
cstab(st,ct,n,ity);
il=(int *)malloc(sizeof(int)*n)-1;
brtab(il,n);
fft1(br-1,bi-1,n,st,ct);
binrv(br-1,n,il);
binrv(bi-1,n,il);
for(int i=0;i<n;i++)
{
br[i]*=invsqrtn;
bi[i]*=invsqrtn;
}
free(il+1);
free(ct+1);
free(st+1);
}


是否我的算法有问题???
...全文
134 7 打赏 收藏 转发到动态 举报
写回复
用AI写文章
7 条回复
切换为时间正序
请发表友善的回复…
发表回复
lanren_me 2003-01-27
  • 打赏
  • 举报
回复
to:hailulu_wang(没啥)
谢谢你的程序:不过跟我原来的运行结果是一样的啊!

理论上,同时有几个相同振幅的频率在才计数居中分析得出来的结果--也就是她么的频谱高度一样才对.

但是现在分析的结果为不一样高.

但是单个频率分析时,得出来的值一样.

比如:单个f的数据分析时得出来的高度都为100;

但是几个频率同时存在时,各个频率对应的高度就不一样了....

是不是就是不一样的???还是一样????

waiting....


hailulu_wang 2003-01-27
  • 打赏
  • 举报
回复
这是我的fft的算法,你看看:
void SimpleDig::Spectrum (float *x,int n)
{
int i;
float temp;
float *pr=new float[n];
float *pi=new float[n];
for (i=0;i<n;i++)
{
pr[i]=x[i];
pi[i]=0.0;
}
fft(pr,pi,n,0);
for (i=0;i<n;i++)
{
temp=hypot(pr[i],pi[i]);
//temp=sqrt(pr[i]*pr[i]+pi[i]*pi[i]);
x[i]=temp*2/n;
}
try{
if(pr) delete [] pr;
if(pi) delete [] pi;
pr=NULL;
pi=NULL;
}
catch (...)
{
ShowMessage("Point delete error");
}
}


void SimpleDig::fft (float *xr,float *xi,int n,int inv)
{
int j,k,ip,l,nm1,nv2,le,le1,nu;
register int i;
float tr,ti,ur,ui,wr,wi,ur1;
//float PI=3.1415926;
j=0;nm1=n-1;nv2=n/2;
nu=num(n);
if(nu==-1) return;
for(i=0;i<=nm1-1;i++)
{
if(i<j)
{
tr=xr[j];
ti=xi[j];
xr[j]=xr[i];
xi[j]=xi[i];
xr[i]=tr;
xi[i]=ti;
}
k=nv2;
while(k-1<j)
{
j=j-k;
k/=2;
}
j=j+k;
}

for(l=1;l<=nu;l++)
{
le=(int)pow(2.0,l);
le1=le/2; ur=1.0; ui=0.0;
wr=(float)cos(PI/le1);
wi=-(float)sin(PI/le1);
if(inv==1) wi=-wi;
for(j=0;j<=le1-1;j++)
{
for(i=j;i<=n-1;i+=le)
{
ip=i+le1;
tr=xr[ip]*ur-xi[ip]*ui;
ti=xr[ip]*ui+xi[ip]*ur;
xr[ip]=xr[i]-tr;
xi[ip]=xi[i]-ti;
xr[i]=xr[i]+tr;
xi[i]=xi[i]+ti;
}
ur1=ur*wr-ui*wi;
ui=ur*wi+ui*wr;
ur=ur1;
}
}
}
hdaq 2003-01-27
  • 打赏
  • 举报
回复
振幅一样的波,分析出来的高度(值)是不是一样还得看频率嘛。

建议你再看一看傅立叶变换理论方面的书籍~~~~~
lanren_me 2003-01-27
  • 打赏
  • 举报
回复
for (int i=0;i<1024;i++)
{
double k=sin(ff1*i);
k+=sin(ff2*i);
ri[i]=k;
}

再说了,从上代码看,两个频率出现的次数应该是一样的.
lanren_me 2003-01-27
  • 打赏
  • 举报
回复
但是已知的条件是高度一样的啊,
因为振幅一样啊.
hdaq 2003-01-27
  • 打赏
  • 举报
回复
几个频率同时存在时,各个频率对应的高度当然不一样了
NowCan 2003-01-27
  • 打赏
  • 举报
回复
?

13,826

社区成员

发帖
与我相关
我的任务
社区描述
C++ Builder相关内容讨论区
社区管理员
  • 基础类社区
加入社区
  • 近7日
  • 近30日
  • 至今
社区公告
暂无公告

试试用AI创作助手写篇文章吧