|
发表于 2011-6-2 13:07:06
|
显示全部楼层
#include<conio.h>
#include<graphics.h>
#include<iostream.h>
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#define PI 3.1415926
// #define N 64
#define N1 64
struct rptr
{
double *r;
double *i;
}COMPLEX;
void fft(double *dr,double *di,int N);
struct rptr *dftn(double a,double h,int n,double *f,int nw,struct rptr *tran);
void draw(double *x,double *h,int lenx,int lenh);
void main()
{
int fp=3;
int i;
struct rptr *tran;
double *fw,*f;
double time=64;//time of signal
double h=0.25;
double a=-time/2;
int N=time/h;
tran->i=(double *)calloc(N+2,sizeof(double));
tran->r=(double *)calloc(N+2,sizeof(double));
fw=(double *)calloc(N+2,sizeof(double));
f=(double *)calloc(N+2,sizeof(double));
cout<<"\nPlease select the input function:\n";
cout<<"1.cos(0.125*PI*t) 2.sin(PI*t)/(PI*t)\n";
cout<<"3.sin(0.25*PI*t)+2*cos(0.5*PI*t) 4.cos(0.125*PI*t)+2*cos(0.25*PI*t)\n";
cout<<"5.sin(0.25*PI*t)\n<-----------";
lp: cin>>fp;
for(i=0;i<=N;i++)
{
double t=(a+h*i);
switch(fp)
{
case 1: f=tran->r=cos(0.125*PI*t);tran->i=0;
case 2: if(t==0) f=tran->r=1;
else
f=tran->r=sin(1.0*PI*t)/(1.0*PI*t);
tran->i=0;
break;
case 3:
f=tran->r=sin(0.25*PI*t)+2*cos(0.5*PI*t);tran->i=0;
break;
case 4:
f=tran->r=cos(0.125*PI*t)+2*cos(0.25*PI*t);tran->i=0;
break;
case 5: f=tran->r=sin(0.25*PI*t);tran->i=0;
break;
default: cout<<"Wrong!\n";
goto lp;
}
}
fft(tran->r,tran->i,N);
cout<<"\n\n";
for(int j=0;j<N;j++)
{fw[j]=pow(tran->r[j],2.0)+pow(tran->i[j],2.0);
printf("%-10.2f",fw[j]);
}getch();
draw(f,fw,N,N);
}
void fft(double *dr,double *di,int N)
{
int lh,m;
lh=N/2;
m=int(log(N)/log(2)+0.9999);
int j=lh;
int k;
for(int i=1;i<(N-1);i++)
{
if(i<j)
{
double tr,ti;
tr=dr;ti=di;
dr=dr[j];di=di[j];
dr[j]=tr;di[j]=ti;
}
k=lh;
while(j>=k) {j=j-k;k=k/2;}
j=j+k;
}
for(i=1;i<=m;i++)
{
int b=int(pow(2.0,(i-1)));
for(j=0;j<b;j++;)
{
double p;
p=j*pow(2.0,(m-i))*2.0*PI/N;
for(k=j;k<N;)
{
double tr,ti;
tr=dr[k+b]*cos(p)+di[k+b]*sin(p);
ti=di[k+b]*cos(p)-dr[k+b]*sin(p);
dr[k+b]=dr[k]-tr;
di[k+b]=di[k]-ti;
dr[k]=dr[k]+tr;
di[k]=di[k]+ti;
k=int(k+pow(2.0,i));
}
}
}
}
void draw(double *x,double *h,int lenx,int lenh)
{
int i,j,k;
double mx,mh,my;
int gdriver=DETECT,gmode,errorcode;
// errorcode=registerbgidriver(EGAVGA_driver);
initgraph(&gdriver,&gmode,"c:\\borlandc\\bgi");
errorcode=graphresult();
if(errorcode!=grOk)
{
cout<<"Graphics error: %s\n",grapherrormsg(errorcode);
cout<<"Press any key to halt!\n";
getch();
exit(1);
}
setbkcolor(WHITE);
setcolor(2);
line(80,140,560,140);
line(80,140,80,40);
line(80,370,560,370);
line(80,370,80,270);
// line(300,250,600,250);
// line(300,250,300,100);
double heigh=80;
long int a;
int sign=1;
char tmp[40],tmp1[40];
mx=x[0];j=0;
double m=x[0];k=0;
for(i=1;i<lenx;i++)
{ if(x<0) {m=(x<m) ? x:m;if(m==x) k=i;}
else {mx=(x>mx) ? x:mx;if(mx==x) j=i;}
}
if((-1)*m>mx)
{
mx=(-1)*m;
sign=-1;
j=k;
}
a=mx;
if(sign==-1) strcpy(tmp,"-");
else strcpy(tmp,"");
ltoa(a,tmp1,10);
strcat(tmp,tmp1);
strcat(tmp,".");
a=(mx-a)*1000;
itoa(a,tmp1,10);
strcat(tmp,tmp1);
if(sign==-1) outtextxy(220.0*j/lenx,140+heigh,tmp);
else outtextxy(220.0*j/lenx,130-heigh,tmp);
setcolor(1);
mx=heigh/mx;
moveto(80,140-mx*x[0]);
for(i=0;i<lenx;i++)
lineto(80+450.0*i/lenx,140-mx*x);
mh=h[0];j=0;
m=h[0];k=0;
sign=1;
for(i=1;i<lenh;i++)
{
if(h<0) {m=(h<m ? h:m;if(m==h) k=i;}
else {mh=(h>mh) ? h:mh;if(mh==h) j=i;}
}
if((-1)*m>mh)
{
mh=(-1)*m;
sign=-1;
j=k;
}
a=mh;
if(sign==-1) strcpy(tmp,"-");
else strcpy(tmp,"");
ltoa(a,tmp1,10);
strcat(tmp,tmp1);
strcat(tmp,".");
a=(mh-a)*1000;
itoa(a,tmp1,10);
strcat(tmp,tmp1);
setcolor(2);
if(sign==-1) outtextxy(220.0*j/lenh,370+heigh,tmp);
else outtextxy(220.0*j/lenh,358-heigh,tmp);
mh=heigh/mh;
setcolor(1);
moveto(80,370-mh*h[0]);
for(i=0;i<lenh;i++)
{
moveto(80+450.0*i/lenh,370);
lineto(80+450.0*i/lenh,370-mh*h);
}
settextstyle(DEFAULT_FONT,HORIZ_DIR,2);
setcolor(7);
outtextxy(75,15,"Input x(t)");
outtextxy(70,230,"Response h(t)");
// outtextxy(320,68,"Onput x(t)*h(t)");
settextstyle(DEFAULT_FONT,HORIZ_DIR,1);
setcolor(2);
outtextxy(44,38,"x(t)");
outtextxy(44,266,"h(t)");
// outtextxy(302,96,"x(t)*h(t)");
outtextxy(272,140,"t");
outtextxy(272,380,"t");
// outtextxy(602,252,"t");
outtextxy(39,150,"0");
outtextxy(39,374,"0");
// outtextxy(300,254,"0");
getch();
closegraph();
} |
|