?? fsample.c
字號:
/**********************************************************************
FSAMPLE.C - 取樣混迭演示程序
dft 離散傅里葉變換函數
idft IDFT 函數
draw_image 繪圖子程序
***********************************************************************/
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <conio.h>
#include <alloc.h>
#include <graphics.h>
/* COMPLEX STRUCTURE */
typedef struct {
double real, imag;
} COMPLEX;
#define PI (4.0*atan(1.0))
void idft(COMPLEX *Datain, COMPLEX *Dataout, int N);
void dft(COMPLEX *Datain, COMPLEX *Dataout, int N);
void draw_image(double *x,int m,char *title1,char *title2,
char *xdis1,char *xdis2,int dis_type);
/********************************************************/
void main(void)
{
int i,j,k,length,m,type;
char title1[80],title2[80],tmp1[20],tmp2[20];
double far *amp;
COMPLEX far *fsamp, far *tsamp, far *ssamp;
length = 256; /* length is the display sample size */
amp = (double *) farcalloc(length+1,sizeof(double));
ssamp = (COMPLEX *) farcalloc(length+1,sizeof(COMPLEX));
fsamp = (COMPLEX *) farcalloc(length+1,sizeof(COMPLEX));
tsamp = (COMPLEX *) farcalloc(length+1,sizeof(COMPLEX));
if(!tsamp) {
printf("\nUnable to allocate complex array for calculation\n");
exit(1);
}
/* Input frequency sampling data for processing */
for (i = 0; i < length; i++) {
if(i<length/14) fsamp[i].real = 24.0-24.0*14.0*i/length;
else if(i>length-length/14) fsamp[i].real = 24.0*14.0*i/length-24.0*14.0+24.0;
else fsamp[i].real = 0;
fsamp[i].imag=0;
}
printf("Generate the Signal. Waitting for the calculation...\n");
idft(fsamp, ssamp, length);
for(i = 0; i< length; i++)
amp[i]=ssamp[i].real;
strcpy(title1,"The Signal Data Sequency");
strcpy(title2,"The Magnitude of Signal");
strcpy(tmp1,"0");
gcvt((float)length/2/70,5,tmp2);
strcat(tmp2,"(s) t");
draw_image(amp,length,title1,title2,tmp1,tmp2,0);
/* Resample with fs = 10 Hz */
for (i = 0 ; i < length/7 ; i++) {
tsamp[i].real = ssamp[i*7].real;
tsamp[i].imag = ssamp[i*7].imag;
}
printf("Waitting for the calculation...\n");
dft(tsamp, fsamp, (int)(length/7));
for(i = 0; i< length/7; i++)
amp[i]=fsamp[i].real;
strcpy(title1,"The Signal Spectrum fs=10Hz");
strcpy(title2,"The Spectrum Magnitude");
strcpy(tmp1,"0");
strcpy(tmp2,"2*PI");
draw_image(amp,(int)(length/7),title1,title2,tmp1,tmp2,0);
/* Resample with fs = 7 Hz */
for (i = 0 ; i < length/10 ; i++) {
tsamp[i].real = ssamp[i*10].real;
tsamp[i].imag = ssamp[i*10].imag;
}
printf("Waitting for the calculation...\n");
dft(tsamp, fsamp, (int)(length/10));
for(i = 0; i< length/10; i++)
amp[i]=fsamp[i].real;
strcpy(title1,"The Signal Spectrum fs=7Hz");
strcpy(title2,"The Spectrum Magnitude");
strcpy(tmp1,"0");
strcpy(tmp2,"2*PI");
draw_image(amp,(int)(length/10),title1,title2,tmp1,tmp2,0);
/* Resample with fs = 14 Hz */
for (i = 0 ; i < length/5 ; i++) {
tsamp[i].real = ssamp[i*5].real;
tsamp[i].imag = ssamp[i*5].imag;
}
printf("Waitting for the calculation...\n");
dft(tsamp, fsamp, (int)(length/5));
for(i = 0; i< length/5; i++)
amp[i]=fsamp[i].real;
strcpy(title1,"The Signal Spectrum fs=14Hz");
strcpy(title2,"The Spectrum Magnitude");
strcpy(tmp1,"0");
strcpy(tmp2,"2*PI");
draw_image(amp,(int)(length/5),title1,title2,tmp1,tmp2,0);
farfree(fsamp);
farfree(ssamp);
farfree(tsamp);
farfree(amp);
}
/***********************************************************************
dft - 離散傅里葉正變換子程序
輸入參數:
COMPLEX *Datain : 輸入數據區指針;
COMPLEX *Dataout: 輸出數據區指針;
int N : 數據長度;
輸出參數:
輸出數據存放在 Dataout 所指的數據區;
無輸出參數.
void dft(COMPLEX *Datain, COMPLEX *Dataout, int N)
************************************************************************/
void dft(COMPLEX *Datain, COMPLEX *Dataout, int N)
{
int i,k,n,p;
static int nstore = 0; /* store N for future use */
static COMPLEX *cf; /* coefficient storage */
COMPLEX *cfptr,*Dinptr;
double arg;
/* Create the coefficients if N has changed */
if(N != nstore) {
if(nstore != 0) free((char *) cf); /* free previous */
cf = (COMPLEX *) calloc(N, sizeof(COMPLEX));
if (!cf) {
printf("\nUnable to allocate memory for coefficients.\n");
exit(1);
}
arg = 8.0*atan(1.0)/N;
for (i=0 ; i<N ; i++) {
cf[i].real = (float)cos(arg*i);
cf[i].imag = -(float)sin(arg*i);
}
}
/* Perform the DFT calculation */
printf("\n");
for (k=0 ; k<N ; k++) {
Dinptr = Datain;
Dataout->real = Dinptr->real;
Dataout->imag = Dinptr->imag;
Dinptr++;
for (n=1; n<N; n++) {
p = (int)((long)n*k % N);
cfptr = cf + p; /* pointer to cf modulo N */
Dataout->real += Dinptr->real * cfptr->real
- Dinptr->imag * cfptr->imag;
Dataout->imag += Dinptr->real * cfptr->imag
+ Dinptr->imag * cfptr->real;
Dinptr++;
}
if (k % 32 == 31) printf("*");
Dataout++; /* next output */
}
printf("\n");
}
/***********************************************************************
idft - 離散傅里葉反變換子程序
輸入參數:
COMPLEX *Datain : 輸入數據區指針;
COMPLEX *Dataout: 輸出數據區指針;
int N : 數據長度;
輸出參數:
輸出數據存放在 Dataout 所指的數據區;
無輸出參數.
void idft(COMPLEX *Datain, COMPLEX *Dataout, int N)
************************************************************************/
void idft(COMPLEX *Datain, COMPLEX *Dataout, int N)
{
int i,k,n,p;
static int nstore = 0; /* store N for future use */
static COMPLEX *cf; /* coefficient storage */
COMPLEX *cfptr,*Dinptr;
double arg;
/* Create the coefficients if N has changed */
if(N != nstore) {
if(nstore != 0) free((char *) cf); /* free previous */
cf = (COMPLEX *) calloc(N, sizeof(COMPLEX));
if (cf == 0) {
printf("\nUnable to allocate memory for coefficients.\n");
exit(1);
}
/* scale stored values by 1/N */
arg = 8.0*atan(1.0)/N;
for (i=0 ; i<N ; i++) {
cf[i].real = (float)(cos(arg*i)/(double)N);
cf[i].imag = (float)(sin(arg*i)/(double)N);
}
}
/* Perform the DFT calculation */
printf("\n");
for (k=0 ; k<N ; k++) {
Dinptr = Datain;
Dataout->real = Dinptr->real * cf[0].real;
Dataout->imag = Dinptr->imag * cf[0].real;
Dinptr++;
for (n=1; n<N; n++) {
p = (int)((long)n*k % N);
cfptr = cf + p; /* pointer to cf modulo N */
Dataout->real += Dinptr->real * cfptr->real
- Dinptr->imag * cfptr->imag;
Dataout->imag += Dinptr->real * cfptr->imag
+ Dinptr->imag * cfptr->real;
Dinptr++;
}
if (k % 32 == 31) printf("*");
Dataout++; /* next output */
}
printf("\n");
}
/************************************************************************
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 + -