?? dfstest.c
字號:
/**********************************************************************
DFSTEST.C - 離散傅里葉級數合成連續周期信號
gen_freq 產生信號的頻域數據
cal_time 計算信號的時域波形
draw_image 輸入數據繪圖子程序
***********************************************************************/
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <conio.h>
#include <graphics.h>
/* function prototypes and structures for fft functions */
/* COMPLEX STRUCTURE */
typedef struct {
double real, imag;
} COMPLEX;
#define PI (4.0*atan(1.0))
void gen_freq(COMPLEX *x, int m ,double a, double w0,int type);
void cal_time(COMPLEX *x, int m, 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 *amp;
double a,w0;
COMPLEX *fsamp,*tsamp;
a = 1; /* a : 信號幅度 */
w0 = 2.0*PI*1; /* w0: 信號的頻率(弧度) */
m = 10; /* m : DFS 長度 */
length = 512; /* length: 顯示樣點數據長度 */
type = 3; /* type: 對應于表2.2 中的信號類型(1,2,...,6) */
amp = (double *) calloc(length+1,sizeof(double));
fsamp = (COMPLEX *) calloc(m*2+2,sizeof(COMPLEX));
tsamp = (COMPLEX *) calloc(length+1,sizeof(COMPLEX));
if(!tsamp) {
printf("\nUnable to allocate complex array for calculation\n");
exit(1);
}
/* Input frequency DFS sampling data for processing */
gen_freq(fsamp,m,a,w0,type);
for (i = 0; i < 2*m+1; i++) {
amp[i]= sqrt(fsamp[i].real*fsamp[i].real+fsamp[i].imag*fsamp[i].imag);
}
strcpy(title1,"The DFS Data Sequency N=");
strcat(title1,itoa(m,tmp1,10));
strcpy(title2,"The Magnitude of DFS");
strcpy(tmp1,"-N");
strcpy(tmp2,"N");
draw_image(amp,2*m,title1,title2,tmp1,tmp2,1);
/* Calculate the time domain signal wave of the input DFS*/
printf("Waitting for the calculation...\n");
cal_time(fsamp,m,tsamp,length);
for (i = 0 ; i < length ; i++) {
amp[i] = tsamp[i].real;
}
strcpy(title1,"The Signal In Time Domain N=");
strcat(title1,itoa(m,tmp1,10));
strcpy(title2,"The Time Signal Magnitude");
strcpy(tmp1,"0");
gcvt(4.0*PI/w0,6,tmp2);
strcat(tmp2,"(s)");
draw_image(amp,length,title1,title2,tmp1,tmp2,0);
free(fsamp);
free(tsamp);
free(amp);
}
/**************************************************************************
gen_freq - 信號的頻率域數據產生子程序
輸入參數:
x : 復數指針,作為輸出頻域信號的存放區
m : 整數,輸出頻域信號的長度(length=2^m+1)
a : 信號的幅度, 雙精度浮點數
w0: 信號基波的頻率(弧度),雙精度浮點數
type: 產生信號的類型, 對應于表2.2(1,2,...,6)
輸出參數:
產生的信號存放在 x 所指定的區域
void gen_freq(COMPLEX *x, int m ,double a, double w0,int type)
*************************************************************************/
void gen_freq(COMPLEX *x, int m, double a, double w0,int type)
{
int i,j,le;
le = 2*m + 1;
/* calculate the F-Domain X(nw0) values recursively */
switch(type) {
case 1: i= -m;
for (j=0; j<le; j++) {
if(abs(i%2) == 0) {
(x+j)->real = 2.0*a/(PI*(1.0-i*i));
(x+j)->imag = 0;
}
else {
(x+j)->real = 0;
(x+j)->imag = 0;
}
i++;
}
break;
case 2: i= -m;
for (j=0; j<le; j++) {
if(abs(i%2) == 0) {
(x+j)->real = a/(PI*(1.0-i*i));
(x+j)->imag = 0;
}
else if(abs(i) != 1) {
(x+j)->real = 0;
(x+j)->imag = 0;
}
else {
(x+j)->real = 0;
(x+j)->imag = -1.0*i*a/4.0;
}
i++;
}
break;
case 3: i= -m;
for (j=0; j<le; j++) {
if(abs(i)%2 == 0) {
(x+j)->real = 0;
(x+j)->imag = 0;
}
else if( abs(i>>1)%2 == 0 ) {
(x+j)->real = 2.0*a/i/PI;
(x+j)->imag = 0;
}
else {
(x+j)->real = -2.0*a/i/PI;
(x+j)->imag = 0;
}
i++;
}
break;
case 4: i= -m;
for (j=0; j<le; j++) {
if(abs(i%2) == 0) {
(x+j)->real = 0;
(x+j)->imag = 0;
}
else {
(x+j)->real = 4.0*a*cos(i*PI/2)/(PI*PI*i*i);
(x+j)->imag = -4.0*a*sin(i*PI/2)/(PI*PI*i*i);
}
i++;
}
break;
case 5: i= -m;
for (j=0; j<le; j++) {
if(i == 0) {
(x+j)->real = 0;
(x+j)->imag = 0;
}
else {
(x+j)->real = a*pow(-1,i+1)*cos(PI/2)/(i*PI);
(x+j)->imag = -a*pow(-1,i+1)*sin(PI/2)/(i*PI);
}
i++;
}
break;
case 6: for (j=0; j<le; j++) {
(x+j)->real = w0/2.0/PI;
(x+j)->imag = 0;
}
break;
default: for (j=0; j<le; j++) {
(x+j)->real = 0;
(x+j)->imag = 0;
}
break;
}
}
/***********************************************************************
cal_time - 計算信號的時間域波形子程序
輸入參數:
x : 輸入頻率域的數據指針,復數指針;
m : 輸入數據的長度, 整數;
Dataout: 輸出數據指針, 復數指針;
N : 輸出數據長度;
輸出參數:
輸出數據在 Dataout 中, 長度為 N;
無輸出參數.
void cal_time(COMPLEX *x, int m, COMPLEX *Dataout, int N)
*************************************************************************/
void cal_time(COMPLEX *x, int m, COMPLEX *Dataout, int N)
{
int i,j,k;
COMPLEX *Dinptr,*Doutptr;
double arg;
/* Perform the DFT calculation */
arg = 4*PI/N;
printf("\n");
Doutptr=Dataout;
for (i=0; i<N; i++) {
Doutptr->real = 0;
Doutptr->imag = 0;
Doutptr++;
}
Dinptr = x;
for (i=-m ; i<=m ; i++) {
Doutptr=Dataout;
for (j=0; j<N; j++) {
Doutptr->real += Dinptr->real*cos(arg*i*j)-Dinptr->imag*sin(arg*i*j);
Doutptr->imag += Dinptr->imag*cos(arg*i*j)+Dinptr->real*sin(arg*i*j);
Doutptr++;
}
printf("*");
Dinptr++;
}
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 + -