?? dif_fft.c
字號:
/**********************************************************************
DIF_FFT.C - DIF FFT 測試程序
fft 基2 DIF FFT 子程序
ifft 基2 DIF IFFT 子程序
draw_image 繪圖子程序
***********************************************************************/
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <conio.h>
#include <graphics.h>
/* COMPLEX STRUCTURE */
typedef struct {
float real, imag;
} COMPLEX;
#define PI (4.0*atan(1.0))
void fft(COMPLEX *,int);
void ifft(COMPLEX *,int);
void draw_image(double *x,int m,char *title1,char *title2,
char *xdis1,char *xdis2,int dis_type);
/********************************************************/
void main(void)
{
int i,length,m,j;
char title[80],tmp[20];
double *amp;
double a,tempflt;
COMPLEX *samp;
m=8;
length = 1<<m;
amp = (double *) calloc(length+1,sizeof(double));
samp = (COMPLEX *) calloc(length+1,sizeof(COMPLEX));
if(!samp) {
printf("\nUnable to allocate complex array for fft\n");
exit(1);
}
/* Input sampling data for processing */
printf("Waitting for sampling data...");
for (i = 0; i < length; i++) {
amp[i] = cos(2.0*PI*(double)i*2.0/length);
printf("*");
}
for (i = 0; i < length; i++) {
samp[i].real = amp[i]*cos(16.0*PI*(double)i*2.0/length)*amp[i];
amp[i]= samp[i].real;
printf("*");
}
strcpy(title,"The Sampling Signal Data");
draw_image(amp,length,title,"The Magnitude","0",itoa(length,tmp,10),0);
/* Find the spectrum of the data */
printf("Waitting for the FFT calculation...\n");
fft(samp,m);
/* calculate magnitude */
tempflt=0;
for (i = 0 ; i < length ; i++) {
tempflt = samp[i].real * samp[i].real + samp[i].imag * samp[i].imag;
amp[i] = tempflt;
}
strcpy(title,"The Signal Freq Magnitude Result");
draw_image(amp,length,title,"The Magnitude","0",itoa(length,tmp,10),0);
printf("Waitting for the IFFT calculation...\n");
ifft(samp,m);
strcpy(title,"The IFFT Result");
for (i=0;i<length;i++) amp[i]=samp[i].real;
draw_image(amp,length,title,"The Magnitude","0",itoa(length,tmp,10),0);
free(samp);
free(amp);
}
/************************************************************************
fft - 基2 DIF FFT 子程序
輸入參數:
COMPLEX *x : FFT 輸入和輸出數據區指針;
int m : FFT 長度 ( length = 2^m );
輸出參數:
輸出數據放在 x 所指的輸入數據區.
無輸出參數.
void fft(COMPLEX *x, int m)
*************************************************************************/
void fft(COMPLEX *x,int m)
{
static COMPLEX *w; /* used to store the w complex array */
static int mstore = 0; /* stores m for future reference */
static int n = 1; /* length of fft stored for future */
COMPLEX u,temp,tm;
COMPLEX *xi,*xip,*xj,*wptr;
int i,j,k,l,le,windex;
double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;
if(m != mstore) {
/* free previously allocated storage and set new m */
if(mstore != 0) free(w);
mstore = m;
if(m == 0) return; /* if m=0 then done */
/* n = 2**m = fft length */
n = 1 << m;
le = n/2;
/* allocate the storage for w */
w = (COMPLEX *) calloc(le-1,sizeof(COMPLEX));
if(!w) {
printf("\nUnable to allocate complex W array\n");
exit(1);
}
/* calculate the w values recursively */
arg = 4.0*atan(1.0)/le; /* PI/le calculation */
wrecur_real = w_real = cos(arg);
wrecur_imag = w_imag = -sin(arg);
xj = w;
for (j = 1 ; j < le ; j++) {
xj->real = (float)wrecur_real;
xj->imag = (float)wrecur_imag;
xj++;
wtemp_real = wrecur_real*w_real - wrecur_imag*w_imag;
wrecur_imag = wrecur_real*w_imag + wrecur_imag*w_real;
wrecur_real = wtemp_real;
}
}
/* start fft */
le = n;
windex = 1;
for (l = 0 ; l < m ; l++) {
le = le/2;
/* first iteration with no multiplies */
for(i = 0 ; i < n ; i = i + 2*le) {
xi = x + i;
xip = xi + le;
temp.real = xi->real + xip->real;
temp.imag = xi->imag + xip->imag;
xip->real = xi->real - xip->real;
xip->imag = xi->imag - xip->imag;
*xi = temp;
}
/* remaining iterations use stored w */
wptr = w + windex - 1;
for (j = 1 ; j < le ; j++) {
u = *wptr;
for (i = j ; i < n ; i = i + 2*le) {
xi = x + i;
xip = xi + le;
temp.real = xi->real + xip->real;
temp.imag = xi->imag + xip->imag;
tm.real = xi->real - xip->real;
tm.imag = xi->imag - xip->imag;
xip->real = tm.real*u.real - tm.imag*u.imag;
xip->imag = tm.real*u.imag + tm.imag*u.real;
*xi = temp;
}
wptr = wptr + windex;
}
windex = 2*windex;
}
/* rearrange data by bit reversing */
j = 0;
for (i = 1 ; i < (n-1) ; i++) {
k = n/2;
while(k <= j) {
j = j - k;
k = k/2;
}
j = j + k;
if (i < j) {
xi = x + i;
xj = x + j;
temp = *xj;
*xj = *xi;
*xi = temp;
}
}
}
/************************************************************************
ifft - 基2 DIF IFFT 子程序
輸入參數:
COMPLEX *x : IFFT 輸入和輸出數據區指針;
int m : IFFT 長度 ( length = 2^m );
輸出參數:
輸出數據放在 x 所指的輸入數據區.
無輸出參數.
void ifft(COMPLEX *x, int m)
*************************************************************************/
void ifft(x,m)
COMPLEX *x;
int m;
{
static COMPLEX *w; /* used to store the w complex array */
static int mstore = 0; /* stores m for future reference */
static int n = 1; /* length of ifft stored for future */
COMPLEX u,temp,tm;
COMPLEX *xi,*xip,*xj,*wptr;
int i,j,k,l,le,windex;
double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;
float scale;
if(m != mstore) {
/* free previously allocated storage and set new m */
if(mstore != 0) free(w);
mstore = m;
if(m == 0) return; /* if m=0 then done */
/* n = 2**m = inverse fft length */
n = 1 << m;
le = n/2;
/* allocate the storage for w */
w = (COMPLEX *) calloc(le-1,sizeof(COMPLEX));
if(!w) {
printf("\nUnable to allocate complex W array\n");
exit(1);
}
/* calculate the w values recursively */
arg = 4.0*atan(1.0)/le; /* PI/le calculation */
wrecur_real = w_real = cos(arg);
wrecur_imag = w_imag = sin(arg); /* opposite sign from fft */
xj = w;
for (j = 1 ; j < le ; j++) {
xj->real = (float)wrecur_real;
xj->imag = (float)wrecur_imag;
xj++;
wtemp_real = wrecur_real*w_real - wrecur_imag*w_imag;
wrecur_imag = wrecur_real*w_imag + wrecur_imag*w_real;
wrecur_real = wtemp_real;
}
}
/* start inverse fft */
le = n;
windex = 1;
for (l = 0 ; l < m ; l++) {
le = le/2;
/* first iteration with no multiplies */
for(i = 0 ; i < n ; i = i + 2*le) {
xi = x + i;
xip = xi + le;
temp.real = xi->real + xip->real;
temp.imag = xi->imag + xip->imag;
xip->real = xi->real - xip->real;
xip->imag = xi->imag - xip->imag;
*xi = temp;
}
/* remaining iterations use stored w */
wptr = w + windex - 1;
for (j = 1 ; j < le ; j++) {
u = *wptr;
for (i = j ; i < n ; i = i + 2*le) {
xi = x + i;
xip = xi + le;
temp.real = xi->real + xip->real;
temp.imag = xi->imag + xip->imag;
tm.real = xi->real - xip->real;
tm.imag = xi->imag - xip->imag;
xip->real = tm.real*u.real - tm.imag*u.imag;
xip->imag = tm.real*u.imag + tm.imag*u.real;
*xi = temp;
}
wptr = wptr + windex;
}
windex = 2*windex;
}
/* rearrange data by bit reversing */
j = 0;
for (i = 1 ; i < (n-1) ; i++) {
k = n/2;
while(k <= j) {
j = j - k;
k = k/2;
}
j = j + k;
if (i < j) {
xi = x + i;
xj = x + j;
temp = *xj;
*xj = *xi;
*xi = temp;
}
}
/* scale all results by 1/n */
scale = (float)(1.0/n);
for(i = 0 ; i < n ; i++) {
x->real = scale*x->real;
x->imag = scale*x->imag;
x++;
}
}
/************************************************************************
draw_image - 將輸入數據的幅度畫出圖形。該函數可自動調整顯示的比例, 使圖形
充滿整個屏幕。
輸入參數: double *x - 輸入數據序列的指針;
int m - 輸入數據序列的長度;
char *title1 - 顯示圖形的上標題字符串指針;
char *xdis1 - X 坐標左邊顯示標題字符串指針.
char *title2 - 顯示圖形的左標題字符串指針.
char *xdis2 - X 坐標右邊顯示標題字符串指針.
int dis_type - 顯示類型, 0:連線 1:直線.
輸出參數: 無
*************************************************************************/
void draw_image(double *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;
double 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=(double)((scy - start_y - end_y)/ym);
xs=(double)(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();
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -