📄 c_dsp_fft.c
字号:
＃i nclude "stdio.h"
＃i nclude "graphics.h"
＃i nclude "math.h"
＃i nclude "dos.h"
#define MAX_WIDTH 640
#define MAX_HEIGHT 480
#define NUM 32
#define N 31
void conv(float x[],int M,float h[],int N1,float y[]);
void dft(float x[],int m,float rx[],float cx[],float mod[]);
void idft(float rx3[],float cx3[],float mod3[],int m,float rx[],float cx[],float mod[]);
void fft(float x[],int m,float rx[],float cx[],float mod[]);
void inition();
void figure(float x[],int m);
void first();
void second();
void third();
void fourth();
void fifth();
void sixth();
void seventh();
void eighth();
void test();
float power(float m,float n);
void xh(float x[],float h[]);
struct BITMAPFILEHEADER /*位图文件头*/
{
char bfType[2]; /*文件类型说明,各占一个字节,windows系统中为BM*/
long bfSize; /*位图文件的大小*/
int bfReserved1;/*保留字节,都为零,共两个保留*/
int bfReserved2;
long bfOffBits; /*从文件头开始到实际的图象数据间的字节偏移量*/
}BMPHeader={"BM",0,0,0,118};
/*位图信息头由位图信息头(bitmap-informationheader)和彩色表(color table)组成*/
struct BITMAPINFO
{
/*位图信息头*/
long biSize; /*位图信息头的大小,共四字节*/
long biWidth; /*图象的宽度和高度,各占四个字节*/
long biHeight;
int biPlanes; /*为目标设备说明位面数,总为1,占两个字节*/
int biBitCount; /*说明比特数/象素,其值为1、4、8、16、24、32占两字节*/
long biCompression; /*说明压缩类型,共四个字节,0为没压缩??*/
long biSizeImage; /*说明图象的大小,以字节为单位,共四个字节*/
long biXPelsPerMeter; /*说明水平分辨率,用象素/米表示,共四字节*/
long biYPelsPerMeter; /*说明垂直分辨率,用象素/米表示,共四字节*/
long biClrUsed; /*位图实际使用的彩色表中的颜色索引数,为0表使用所有调色板项*/
long biClrImportant;/*对图象显示有重要影响的颜色索引数目,为0表都重要四字节*/
}BMPInfo={40,MAX_WIDTH,MAX_HEIGHT,1,4,0,0,3780,3780,0,0};
/*color information table,set blue,green and red intensity*/
struct PALETTE_BOARD
{
int rgbBlue;
int rgbGreen;
int rgbRed;
int rgbReserved;
}Pal[16]=
{
{0,0,255,0},{255,0,0,0},{0,255,0,0},{0,255,255,0},{128,128,255,0},{255,255,0,0},
{255,255,255,0},{0,0,0,0},{128,128,128,0},{192,192,192,0},{128,0,255,0},
{128,255,128,0},{255,128,255,0},{255,128,0,0},{255,255,136,0},{0,64,128,0}
};
/*依次为红,蓝,绿,黄,浅红,青,白,黑,深灰,浅灰,洋红,浅绿,浅洋红,浅蓝,浅青,棕色*/
/*根据扫描得到的颜色与彩色信息表对照,并返回颜色索引号(即数组序号)*/
unsigned int get_index(int x,int y)
{
unsigned int index_num;
switch(getpixel(x,y))
{
case RED:index_num=0;break;
case BLUE:index_num=1;break;
case GREEN:index_num=2;break;
case YELLOW:index_num=3;break;
case LIGHTRED:index_num=4;break;
case CYAN:index_num=5;break;
case WHITE:index_num=6;break;
case BLACK:index_num=7;break;
case DARKGRAY:index_num=8;break;
case LIGHTGRAY:index_num=9;break;
case MAGENTA:index_num=10;break;
case LIGHTGREEN:index_num=11;break;
case LIGHTMAGENTA:index_num=12;break;
case LIGHTBLUE:index_num=13;break;
case LIGHTCYAN:index_num=14;break;
case BROWN:index_num=15;break;
}
return index_num;
}
void scan_screen(FILE *fp)
{
int i,j;
union COLOR /*two dots hold one byte*/
{
unsigned char clr;
struct
{
unsigned cl:4;
unsigned ch:4;
}mycolor;
}color;
for(i=MAX_HEIGHT-1;i>=0;i--)
for(j=0;j<MAX_WIDTH-1;j+=2)/*step is 2,for size of two dots is 1 byte*/
{
color.mycolor.ch=get_index(j,i); /*first dot*/
color.mycolor.cl=get_index(j+1,i); /*next dot*/
fwrite(&color.clr,1,1,fp);
}
/*sound(1500);
delay(4000);
nosound();*/
}
void Screen_Save(char *filename)
{
FILE *fp;
int i;
if((fp=fopen(filename,"wb"))==NULL)exit(0);
/*write bitmap file information*/
fwrite(BMPHeader.bfType,2,1,fp);
BMPInfo.biSizeImage=BMPInfo.biWidth*BMPInfo.biHeight*BMPInfo.biBitCount/8;
BMPHeader.bfSize=BMPInfo.biSizeImage+BMPHeader.bfOffBits;
fwrite(&BMPHeader.bfSize,4,1,fp);
fwrite(&BMPHeader.bfReserved1,2,1,fp);
fwrite(&BMPHeader.bfReserved2,2,1,fp);
fwrite(&BMPHeader.bfOffBits,4,1,fp);
/*写位图信息头和色彩表信息*/
BMPInfo.biSize=sizeof(struct BITMAPINFO);
fwrite(&BMPInfo.biSize,4,1,fp);
fwrite(&BMPInfo.biWidth,4,1,fp);
fwrite(&BMPInfo.biHeight,4,1,fp);
fwrite(&BMPInfo.biPlanes,2,1,fp);
fwrite(&BMPInfo.biBitCount,2,1,fp);
fwrite(&BMPInfo.biCompression,4,1,fp);
fwrite(&BMPInfo.biSizeImage,4,1,fp);
fwrite(&BMPInfo.biXPelsPerMeter,4,1,fp);
fwrite(&BMPInfo.biYPelsPerMeter,4,1,fp);
fwrite(&BMPInfo.biClrUsed,4,1,fp);
fwrite(&BMPInfo.biClrImportant,4,1,fp);
/*write color information table*/
for(i=0;i<16;i++)
{
fwrite(&Pal[i].rgbBlue,1,1,fp);
fwrite(&Pal[i].rgbGreen,1,1,fp);
fwrite(&Pal[i].rgbRed,1,1,fp);
fwrite(&Pal[i].rgbReserved,1,1,fp);
}
/*写位图数据*/
scan_screen(fp);
fclose(fp);
}
main()
{
int input1;
inition();
init:
gotoxy(14,2);printf("welcome!!!");
gotoxy(14,4);printf("1:plot x(n)");
gotoxy(14,6);printf("2:plot R32(n)");
gotoxy(14,8);printf("3:plot dftx");
gotoxy(14,10);printf("4:plot dfth");
gotoxy(14,12);printf("5:plot yxh");
gotoxy(14,14);printf("6:idfty");
gotoxy(14,16);printf("7:FFT16");
gotoxy(14,18);printf("8:FFT32");
gotoxy(14,20);printf("0:exit");
gotoxy(14,22);printf("please input the number:");
scanf("%d",&input1);
switch(input1)
{
case 1:clrscr();first();goto init;
case 2:clrscr();second();goto init;
case 3:clrscr();third();goto init;
case 4:clrscr();fourth();goto init;
case 5:clrscr();fifth();goto init;
case 6:clrscr();sixth();goto init;
case 7:clrscr();seventh();goto init;
case 8:clrscr();eighth();goto init;
case 0:exit(0);
default:clrscr();printf("error,please input again:");goto init;
}
getch();
}
void inition()
{
int gdriver,gmode;
gdriver=VGA;
gmode=VGAHI;
initgraph(&gdriver,&gmode,"");
}
void xh(float x[],float h[])
{
int i;
for(i=0;i<=31;i++)
{
x[i]=i+1;
h[i]=1;
}
}
float power(float m,float n)
{
int i;
float t=1;
for(i=0;i<n;i++)
t=t*m;
return t;
}
multiply(float rx1[],float cx1[],int num1,float rx2[],float cx2[],int num2,float rx3[],float cx3[],float mod3[])
{
int m;
for(m=0;m<num1;m++)
{
rx3[m]=rx1[m]*rx2[m]-cx1[m]*cx2[m];
cx3[m]=rx1[m]*cx2[m]+rx2[m]*cx1[m];
mod3[m]=sqrt(rx3[m]*rx3[m]+cx3[m]*cx3[m]);
}
}
void conv(float x[],int M,float h[],int N1,float y[])
{int m,n;
double w=0;
for(n=0;n<M+N1-1;n++)
{w=0;
for(m=0;m<M;m++)
{if((n-m)>=0&&((n-m)<N1))
w+=x[m]*h[n-m];
}
y[n]=w;
}
}
void dft(float x[],int m,float rx[],float cx[],float mod[])
{
int k,n;
for(k=0;k<=m-1;k++)
{
rx[k]=0;
cx[k]=0;
for(n=0;n<=m-1;n++)
{
rx[k]=rx[k]+x[n]*(cos(-2*3.14159265359798*n*k/m));
cx[k]=cx[k]+x[n]*(sin(-2*3.14159265359798*n*k/m));
}
mod[k]=sqrt(rx[k]*rx[k]+cx[k]*cx[k]);
}
}
void idft(float rx3[],float cx3[],float mod3[],int m,float rx[],float cx[],float mod[])
{
int k,n;
for(n=0;n<=m-1;n++)
{
rx[n]=0;
cx[n]=0;
for(k=0;k<=m-1;k++)
{
rx[n]=rx[n]+rx3[k]*cos(2*3.14159265359798*n*k/m)-cx3[k]*sin(2*3.14159265359798*n*k/m);
cx[n]=cx[n]+rx3[k]*sin(2*3.14159265359798*n*k/m)+cx3[k]*cos(2*3.14159265359798*n*k/m);
}
rx[n]=rx[n]/m;
cx[n]=cx[n]/m;
mod[n]=sqrt(rx[n]*rx[n]+cx[n]*cx[n]);
}
}
void fft(float x[],int m,float rx[],float cx[],float mod[])
{
int k,n;
float rx1[100],rx2[100],cx1[100],cx2[100],tempr[100],tempc[100],wr[100],wc[100];
for(k=0;k<m/2;k++)
{
rx1[k]=0;
cx1[k]=0;
rx2[k]=0;
cx2[k]=0;
for(n=0;n<m-1;n+=2)
{
rx1[k]=rx1[k]+x[n]*(cos(-2*3.14159265359798*(n/2)*k/(m/2)));
cx1[k]=cx1[k]+x[n]*(sin(-2*3.14159265359798*(n/2)*k/(m/2)));
} /*奇数的求x x1*/
for(n=1;n<m-1;n+=2)
{
rx2[k]=rx2[k]+x[n]*(cos(-2*3.14159265359798*((n-1)/2)*k/(m/2)));
cx2[k]=cx2[k]+x[n]*(sin(-2*3.14159265359798*((n-1)/2)*k/(m/2)));
} /*偶数的求x x2*/
wr[k]=cos(-2*3.14159265359798*k/m);
wc[k]=sin(-2*3.14159265359798*k/m); /*Wn(k) 7-11公式中*/
tempr[k]=rx2[k];
tempc[k]=cx2[k]; /*中间变量*/
rx2[k]=tempr[k]*wr[k]-tempc[k]*wc[k];
cx2[k]=tempr[k]*wc[k]+tempc[k]*wr[k]; /*变化后的x2*/
rx[k]=rx1[k]+rx2[k];
cx[k]=cx1[k]+cx2[k]; /*最后值：x*/
mod[k]=sqrt(rx[k]*rx[k]+cx[k]*cx[k]); /*first part have compute */
}
/* for(k=0;k<m/2;k++)
{
wr[k+m/2]=cos(-2*3.14159265359798*k/m);
wc[k+m/2]=sin(-2*3.14159265359798*k/m);
tempr[k+m/2]=rx2[k];
tempc[k+m/2]=cx2[k];
rx2[k+m/2]=tempr[k+m/2]*wr[k+m/2]-tempc[k+m/2]*wc[k+m/2];
cx2[k+m/2]=tempr[k+m/2]*wc[k+m/2]+tempc[k+m/2]*wr[k+m/2];
rx[k+m/2]=rx1[k]-rx2[k];
cx[k+m/2]=cx1[k]-cx2[k];
mod[k+m/2]=sqrt(rx[k+m/2]*rx[k+m/2]+cx[k+m/2]*cx[k+m/2]);
} */
for(k=m/2;k<m;k++)
{
rx[k]=rx1[k-m/2]-rx2[k-m/2];
cx[k]=cx1[k-m/2]-cx2[k-m/2];
mod[k]=sqrt(rx[k]*rx[k]+cx[k]*cx[k]);
}
/*第二部分*/
}
/* void figure(float x[],int m)
{
int i;
float mul=0;
setcolor(WHITE);
line(10,240,800,240);
line(10,10,10,600);
for(i=0;i<m;i++)
{ if(x[i]>mul)
mul=x[i];
}
mul=mul/240;
for(i=0;i<m;i++)
{setcolor(RED);
line(600/m*i+10,240-x[i]/mul,600/m*i+10,240);
}
}*/
void figure(float save[],int N1)
{
int m,i,x,y,t,x0=50,x1=549,y0=229,y1=30,y2=430,wide=50;
float max=0,min=0,addy,mul=0;
cleardevice();
/* y max =??? */
for(i=0;i<N1;i++)
{ if(save[i]>mul)
mul=save[i];
}
addy=mul/200;
/* zuo biao */
gotoxy(x0/8,15);
printf("0");
for(m=1;m<N1+1;m++)
{x=(x1-x0+1)/N1;
gotoxy((x0+x*m)/8,1);
if(m%5==0)
printf("%.1d",m);
}
for(m=-3;m<=4;m++)
{ if(m==0) continue;
gotoxy(1,15-wide*m/16);
printf("%1.2f",m*mul*wide/(y0-y1+1));
}
setcolor(9);
line(x0,y2+30,x0,y1-30);
line(x0,y1-30,x0+3,y1-30+5);
line(x0,y1-30,x0-3,y1-30+5);
line(x0,y0,x1+50, y0);
line(x1+50,y0,x1+50-5,y0+3);
line(x1+50,y0,x1+50-5,y0-3);
/* ke du xian */
for(m=-3;m<5;m++)
line(x0,y0-m*wide,x0+3,y0-m*wide);
for(m=1;m<N1+1;m++)
{ x=(x1-x0+1)/N1;
line(x0+m*x,y1-3,x0+m*x,y1-5);}
gotoxy(1,1);puts("y(n)");
gotoxy(75,16);puts("n");
for(m=0;m<N1+1;m++)
{ t=(x1-x0+1)/N1;
x=m*t+x0;
y=y0-save[m]/addy;
setcolor(10);
line(x,y0,x,y);
}
getch();
}
void first()
{
float x[50],h[50];
xh(x,h);
figure(x,NUM);
Screen_Save("dft1.bmp");
cleardevice();
}
void second()
{
float x[50],h[50];
xh(x,h);
figure(h,NUM);
Screen_Save("dft2.bmp");
cleardevice();
}
void third()
{
float x[50],h[50],rx[50],cx[50],mod[50];
xh(x,h);
dft(x,32,rx,cx,mod);
figure(mod,NUM);
Screen_Save("dft3.bmp");
cleardevice();
}
void fourth()
{
float x[50],h[50],rx[50],cx[50],mod[50];
xh(x,h);
dft(h,32,rx,cx,mod);
figure(mod,NUM);
Screen_Save("dft4.bmp");
cleardevice();
}
void fifth()
{
float x[50],h[100],y[100];
int i;
xh(x,h);
conv(x,NUM,h,NUM,y);
figure(y,63);
Screen_Save("dft5.bmp");
cleardevice();
}
void sixth()
{
float x[100]={0},h[100]={0},rx1[100]={0},cx1[100]={0},rx2[100]={0},cx2[100]={0},rx3[100]={0},cx3[100]={0},rx[100]={0},cx[100]={0},mod1[100]={0},mod2[100]={0},mod3[100]={0},mod[100]={0};
int i;
xh(x,h);
dft(x,63,rx1,cx1,mod1);
dft(h,63,rx2,cx2,mod2);
multiply(rx1,cx1,63,rx2,cx2,63,rx3,cx3,mod3);
idft(rx3,cx3,mod3,63,rx,cx,mod);
figure(mod,63);
Screen_Save("dft6.bmp");
cleardevice();
}
void seventh()
{
float x[50],h[50],rx[50],cx[50],mod[50];
xh(x,h);
fft(x,16,rx,cx,mod);
figure(mod,16);
Screen_Save("dft7.bmp");
cleardevice();
}
void eighth()
{
float x[50],h[50],rx[50],cx[50],mod[50];
xh(x,h);
fft(x,32,rx,cx,mod);
figure(mod,NUM);
Screen_Save("dft8.bmp");
cleardevice();
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -