Fourier 变换的实现！

wlyz333 · 发表于 2004-11-26 05:44:10

<

> 
<DIV class=O v:shape="_x0000_s1026">
<DIV >1.Fourier 变换的实现！ </DIV>
<DIV > (1) 组题背景： </DIV>
<DIV > Fourier在工程中得到了非常广泛的应用。因此掌握Fourier的一些基本知识在今后的学习中起到非常重要的作用。 </DIV>
<DIV > (2) 解决问题： </DIV>
<DIV > 计算量的估计问题、连续Fourier变换的性质与特点、离散 Fourier 变换的计算方法，离散 Fourier 变换蝶形图、 FFT 方法的特点等;应用MATLAB来进行仿真实验. </DIV>
<DIV > (3) 参考书: ”数字信号处理”. </DIV>
<DIV >有谁知道这个题具体的怎么做吗？</DIV>
<DIV ></DIV>
<DIV ></DIV></DIV>

aiwa · 发表于 2005-3-17 19:07:24

<

>这是快速傅立叶变换的程序<

>void fft(double xr[],double xi[],int m,int n)
{
int loop1,loop2,loop3,le,le1,ip,i=1;
double wr,wi,tr,ti,ur1,ur,ui,is;
is=log(n)/log(2);
ss(xr,xi,n);
for(loop1=1;loop1<=is;loop1++)
{
ur=1.0;
ui=0.0;
le=pow(2,loop1);
le1=le/2;
wr=cos(PI*m/le1);
wi=-sin(PI*m/le1);
for(loop2=0;loop2<le1;loop2++)
{
for(loop3=loop2;loop3<n;loop3=loop3+le)
{
ip=loop3+le1;
tr=ur*xr[ip]-ui*xi[ip];
ti=ur*xi[ip]+ui*xr[ip];
xr[ip]=xr[loop3]-tr;
xi[ip]=xi[loop3]-ti;
xr[loop3]=tr+xr[loop3];
xi[loop3]=ti+xi[loop3];
}
ur1=wr*ur-ui*wi;
ui=wr*ui-ur*wi;
ur=ur1;
}
}
if(m==-1)
{
for(i=0;i<n;i++)
{
xr=xr/n;
xi=xi/n;
}
}
}

aiwa · 发表于 2005-3-17 19:10:10

<

>这是我一年前写的一个仿真实验用的程序<

># include <math.h>
# include <stdio.h>
# include <string.h>
# include <conio.h>
# include <stdlib.h>
# include <graphics.h>
# define PI 3.14159<

>void ss(float z[],float y[],int n)
{
int i=0,j,s1;
float a,bj;
for(j=1;j<n;j++)
{
 for(s1=n/2;s1<=i;s1=s1/2)
 {
i=i-s1;
 }
i=i+s1;
 if(i>j)
{
a=z;
bj=y[j];
y=y[j];
z=z[j];
z[j]=a;
y[j]=bj;
}
}
}void fft(float xr[],float xi[],int m,int n)
{
int loop1,loop2,loop3,le,le1,ip,i=1;
float wr,wi,tr,ti,ur1,ur,ui,is;
is=log(n)/log(2);
ss(xr,xi,n);
for(loop1=1;loop1<=is;loop1++)
{
ur=1.0;
ui=0.0;
le=pow(2,loop1);
le1=le/2;
wr=cos(PI*m/le1);
wi=-sin(PI*m/le1);
for(loop2=0;loop2<le1;loop2++)
{
for(loop3=loop2;loop3<n;loop3=loop3+le)
{
ip=loop3+le1;
tr=ur*xr[ip]-ui*xi[ip];
ti=ur*xi[ip]+ui*xr[ip];
xr[ip]=xr[loop3]-tr;
xi[ip]=xi[loop3]-ti;
xr[loop3]=tr+xr[loop3];
xi[loop3]=ti+xi[loop3];
}
ur1=wr*ur-ui*wi;
ui=wr*ui+ur*wi;
ur=ur1;
}
}
if(m==-1)
{
for(i=0;i<n;i++)
{
xr=xr/n;
xi=xi/n;
}
}
}void draw_image(float *x,int m,char *title1,char *title2,
 char *xdis1,char *xdis2,int dis_type)
{
int gdriver=DETECT, gmode,errorcode;
int i,scx,scy,y0,signa,signb;
int style, userpat;
int start_x=40,start_y=40,end_x=10,end_y=60;
long tlen;
float ys,xs,ym;
char dis[40];
/*initializes the graphics mode */
initgraph(&gdriver,&gmode,"");
errorcode=graphresult();
if (errorcode != grOk) {
printf("Graphics error: %s\n",grapherrormsg(errorcode));
printf("Press any key to halt!\n");
getch();
exit(1);
}
scx=getmaxx();
scy=getmaxy();
ym=1.e-90;
signa=0;
signb=0; for(i=0;i<m;i++) {
if ((*(x+i)>0)&&(*(x+i)>ym)) ym = *(x+i);
if ((*(x+i)<0)&&(- *(x+i)>ym)) ym = - *(x+i);
}
for(i=0;i<m;i++) {
if (*(x+i)>fabs(ym/20)) signa=1;
if (*(x+i)<-fabs(ym/20)) signb=1;
}
if ((signa==1)&&(signb==1)) ys=(double)((scy - start_y - end_y)>>1)/ym;
else ys=(float)((scy - start_y - end_y)/ym);
xs=(float)(scx - start_x - end_x)/m;
y0=((scy - start_y - end_y)>>1)+start_y; /* draw the frame */ setcolor(LIGHTGREEN);
rectangle(start_x-1,start_y-20,scx-end_x+1,scy-end_y+20); setcolor(DARKGRAY);
/* select the line style */
style=DASHED_LINE;
userpat = 1;
setlinestyle(style, userpat, 1);
/* a user defined line pattern */
/* binary: "0000000000000001" */
for(i=0;i<=10;i++)
line(start_x,start_y+(scy-start_y-end_y)*i/10,scx-end_x,start_y+(scy-start_y-end_y)*i/10);
for(i=0;i<=10;i++)
line(start_x+(scx-start_x-end_x)*i/10,start_y,start_x+(scx-start_x-end_x)*i/10,scy-end_y);
setcolor(GREEN);
style=SOLID_LINE;
userpat = 1;
setlinestyle(style, userpat, 1);
rectangle(start_x,start_y,scx-end_x,scy-end_y);
setcolor(YELLOW);
for(i=0;i<=10;i++)
line(start_x,start_y+(scy-start_y-end_y)*i/10,start_x+5,start_y+(scy-start_y-end_y)*i/10);
for(i=0;i<=10;i++)
line(start_x+(scx-start_x-end_x)*i/10,scy-end_y+15,start_x+(scx-start_x-end_x)*i/10,scy-end_y+20);
settextstyle(DEFAULT_FONT,HORIZ_DIR,1);
setcolor(YELLOW);
if((signa==1)&&(signb==0)) {
strcpy(dis,"0");
outtextxy(start_x+2,scy-end_y+4,dis);
gcvt(ym,5,dis);
outtextxy(start_x+1,start_y-10,dis);
outtextxy(start_x-10,scy-end_y+24,xdis1);
outtextxy(scx-2-strlen(xdis2)*8,scy-end_y+24,xdis2);
}
else if((signb==1)&&(signa==0)) {
strcpy(dis,"0");
outtextxy(start_x+2,start_y-10,dis);
gcvt(ym,5,dis);
outtextxy(start_x+2,scy-end_y+4,"-");
outtextxy(start_x+10,scy-end_y+4,dis);
outtextxy(start_x-10,scy-end_y+24,xdis1);
outtextxy(scx-2-strlen(xdis2)*8,scy-end_y+24,xdis2);
}
else {
line(start_x,y0,scx-end_x,y0);
strcpy(dis,"0");
outtextxy(start_x-10,y0,dis);
gcvt(ym,5,dis);
outtextxy(start_x+2,start_y-10,dis);
outtextxy(start_x+2,scy-end_y+4,"-");
outtextxy(start_x+10,scy-end_y+4,dis);
outtextxy(start_x-10,scy-end_y+24,xdis1);
outtextxy(scx-2-strlen(xdis2)*8,scy-end_y+24,xdis2);
}
strcpy(dis,"Press any key to continue...");
setcolor(LIGHTRED);
outtextxy((scx-28*8)>>1,scy-16,dis); settextstyle(DEFAULT_FONT,HORIZ_DIR,2);
tlen=strlen(title1);
if ((tlen<<4)<scx) {
setcolor(LIGHTGREEN);
outtextxy((start_x+scx-end_x-(tlen<<4))>>1,start_y-40,title1);
} settextstyle(DEFAULT_FONT,VERT_DIR,1);
tlen=strlen(title2);
if ((tlen<<4)<scy) {
setcolor(LIGHTGREEN);
outtextxy(start_x-20,(scy-end_y-(tlen<<3))>>1,title2);
}
/*draw the amplitude image*/
setcolor(WHITE);
if((signa==1)&&(signb==0)) y0=scy-end_y;
else if((signb==1)&&(signa==0)) y0=start_y;
if (dis_type == 0) {
for(i=0;i<m-1;i++)
 line(xs*i+start_x,y0-*(x+i)*ys,xs*(i+1)+start_x,y0-*(x+i+1)*ys);
}
else if (dis_type == 1) {
for(i=0;i<=m;i++)
 line(xs*i+start_x,y0-*(x+i)*ys,xs*i+start_x,y0);
}
getch();
closegraph();
}
void draw(float x[],int N)
{int i=0;
int graphdriver=DETECT,graphmode;
initgraph(&graphdriver,&graphmode,"");
setcolor(2);
 while(i<N)
 {line(i,10,i,10-10*x);
i++;
 }
setcolor(4);
line(0,10,640,10);
line(640,10,630,15);
line(640,10,630,5);
setlinestyle(CENTER_LINE,0,THICK_WIDTH);
rectangle(0,0,640,480);
getch();
closegraph();
}
void rec_win(float x[],int N)
{int i=0;
while(i<N)
{x=1;
i++;
}
}
void han_win(float x[],int N)
{int i=0;
while(i<N)
{x=0.5*(1+cos(PI*(i-N/2)/(N-1)));
 i++;
 }
}
void ham_win(float x[],int N)
{int i=0;
while(i<N)
{x=0.54+0.46*cos(PI*(i-N/2)/(N-1));
i++;
}
}
void bla_win(float x[],int N)
{int i=0;
while(i<N)
{x=0.42+0.5*cos(PI*(i-N/2)/(N-1))+0.08*cos(2*PI*(i-N/2)/(N-1));
 i++;
}
}
void getmode(float x[],float y[],int N)
{int i=0;
while(i<N)
{x=sqrt(x*x+y*y);
 i++;
}
}
void guiyi(float x[],int N)
{float w=x[0];
int i;
for(i=1;i<N;i++)
if(x>w)w=x;
for(i=0;i<N;i++)
x=x/w;
}
void lg(float x[],int N)
{int i=0;
while(i<N)
{x=20*log10(fabs(x)+0.00000001);
 i++;
}
}
main()
{int i=0,m=0,k=0,N,L,fftN,W,FFTN;
float h1[1024]={0};
float h2[1024]={0};
float h[1024]={0};
float x1[512]={0};
float x2[512]={0};
float y1[512]={0};
float y2[512]={0};
float f,Wd;
char d[12]="Cheng Aihua\0";
char e[12]="Cheng Aihua\0";
char p[12]="Cheng Aihua\0";
char q[12]="Cheng Aihua\0";
printf("\n\n\n ***SHI YAN 4***\n\n\n\n\nPlease input the length of the Hillbert_filter pulse_response now!\n\nWARNING:Shoud be odd_number such as 13,15,17,19,21,23,25,27 and so on\n\nlength=");
scanf("%d",&L); while(i<2*L-1)
{h1=2/((i-L+0.00000001)*PI);
 h2=0;
 i++;
 h1=0;
 h2=0;
 i++;
 }
draw_image(h1,2*L,d,e,p,q,1); printf("\n\nPlease input the sin/cos_signal_length now!\n signal_length=");
scanf("%d",&N);
printf("\n\nPlease input the FFT_points now!\nWARNING:FFT_points shoud not be less than signal_length,\nand only selected from:\n16, 32, 64, 128, 256, 512,!\n FFT_points=");
scanf("%d",&fftN);
printf("\n\n\nPlease input the sin/cos_signal frequence now! f=___? \n\n\n\n\n Made by Cheng Aihua 2003.12.01\n Thanks for using! go on! f=");
scanf("%f",&f); while(k<N)
 {x1[k]=sin(2*PI*f*k);
 x2[k]=cos(2*PI*f*k);
 k++;
}
draw_image(x1,N,d,e,p,q,1);
draw_image(x2,N,d,e,p,q,1); fft(h1,h2,1,fftN);
fft(x1,x2,1,fftN);
while(m<fftN)
{y1[m]=h1[m]*x1[m]-h2[m]*x2[m];
 y2[m]=h1[m]*x2[m]+h2[m]*x1[m];
 m++;
}
i=0;
while(i<fftN)
{h=sqrt(h1*h1+h2*h2);
 i++;
 }
guiyi(h,fftN);
draw_image(h,fftN,d,e,p,q,1);
fft(y1,y2,-1,fftN);
draw_image(y1,fftN,d,e,p,q,1);
draw_image(y2,fftN,d,e,p,q,1);
printf("\n\n ***SHI YAN 2***\n\n\n\nThe program is about to show you four kind of windows one by one:\nRectangle_window,Hanning_window,Hamming_window and blackman_window.\nInput the window_length you design now:window_length=");
scanf("%d",&W);
printf("\n\nPlease input FFT points now!\nWARNING:Never be less than window_length and selected from:\n16, 32, 64, 128, 256, 512, FFTN points=");
scanf("%d",&FFTN);
printf("\n\nPlease input the 'Wd' of the low_pass filter you design now!\nNOTICE:From 0.00001 to 3.14159 Wd=");
scanf("%f",&Wd);i=0;
while(i<512)
{h1=0.0000;
 h2=0.0000;
 i++;
}
rec_win(h1,W);
draw_image(h1,W,d,e,p,q,1);
i=0;
while(i<W)
{h1=(sin(Wd*(i-((W-1)/2)))/(PI*(i-((W-1)/2)+0.0000001)));
 i++;
 }
fft(h1,h2,1,FFTN);
getmode(h1,h2,FFTN);
guiyi(h1,FFTN);
draw_image(h1,FFTN,d,e,p,q,1);
lg(h1,FFTN);
draw(h1,FFTN);
i=0;
while(i<512)
{h1=0.0000;
 h2=0.0000;
 i++;
}
rec_win(h1,W);
fft(h1,h2,1,FFTN);
getmode(h1,h2,FFTN);
guiyi(h1,FFTN);
lg(h1,FFTN);
draw_image(h1,FFTN,d,e,p,q,1);i=0;
while(i<512)
{h1=0.0000;
 h2=0.0000;
 i++;
}
han_win(h1,W);
draw_image(h1,W,d,e,p,q,1);
i=0;
while(i<W)
{h1=h1*(sin(Wd*(i-((W-1)/2)))/(PI*(i-((W-1)/2)+0.0000001)));
 i++;
 }
 fft(h1,h2,1,FFTN);
getmode(h1,h2,FFTN);
guiyi(h1,FFTN);
draw_image(h1,FFTN,d,e,p,q,1);
lg(h1,FFTN);
draw(h1,FFTN);
i=0;
while(i<512)
{h1=0.0000;
 h2=0.0000;
 i++;
}
han_win(h1,W);
fft(h1,h2,1,FFTN);
getmode(h1,h2,FFTN);
guiyi(h1,FFTN);
lg(h1,FFTN);
draw_image(h1,FFTN,d,e,p,q,1);
i=0;
while(i<512)
{h1=0.0000;
 h2=0.0000;
 i++;
}
ham_win(h1,W);
draw_image(h1,W,d,e,p,q,1);
i=0;
while(i<W)
{h1=h1*(sin(Wd*(i-((W-1)/2)))/(PI*(i-((W-1)/2)+0.0000001)));
 i++;
 }
 fft(h1,h2,1,FFTN);
getmode(h1,h2,FFTN);
guiyi(h1,FFTN);
draw_image(h1,FFTN,d,e,p,q,1);
lg(h1,FFTN);
draw(h1,FFTN);
i=0;
while(i<512)
{h1=0.0000;
 h2=0.0000;
 i++;
}
ham_win(h1,W);
fft(h1,h2,1,FFTN);
getmode(h1,h2,FFTN);
guiyi(h1,FFTN);
lg(h1,FFTN);
draw_image(h1,FFTN,d,e,p,q,1);
i=0;
while(i<512)
{h1=0.0000;
 h2=0.0000;
 i++;
}
bla_win(h1,W);
draw_image(h1,W,d,e,p,q,1);
i=0;
while(i<W)
{h1=h1*(sin(Wd*(i-((W-1)/2)))/(PI*(i-((W-1)/2)+0.0000001)));
 i++;
 }
 fft(h1,h2,1,FFTN);
getmode(h1,h2,FFTN);
guiyi(h1,FFTN);
draw_image(h1,FFTN,d,e,p,q,1);
lg(h1,FFTN);
draw(h1,FFTN);
 i=0;
while(i<512)
{h1=0.0000;
 h2=0.0000;
 i++;
}
bla_win(h1,W);
fft(h1,h2,1,FFTN);
getmode(h1,h2,FFTN);
guiyi(h1,FFTN);
lg(h1,FFTN);
draw_image(h1,FFTN,d,e,p,q,1); printf("\n\n\n\n\n\n\n\n\n\n\n\n Made by Cheng Aihua\n\n 2003.12.06\n\n Thanks for using!\n Good Luck!");
scanf("%d",&N);
}

wlyz333 · 发表于 2005-5-14 02:25:35

十分感谢啊！

bsxbs · 发表于 2005-6-25 23:02:38

其中的graphics.h文件怎么回事，本人是学机械方面的。这位大哥写的程序很牛，本人想运行试试啊！！！[em08]

		自动登录	找回密码
密码			注-册-帐-号