169
void bit_reversal(unsigned int *, int , int);
void WTS(float *, float *, int, int);
void FFT(float *xr, float *xi, float *, float *,int,int) ;
void transpose(FILE *, int, int);
void FFT2D(FILE *, FILE *, float *, float *, unsigned int *,
int,int,int);
void main()
{
int N,N1,m,n2,i,j,ii,jj,N2,k,norder,j1,ind,yt;
unsigned int *L;
float *wr,*wi,Do,theta,nsqrt,beta,alpha,sum1,sum2,T;
FILE *fptri,*fptro,*fptr;
float *buffi,*buffo,R2,H,nsq,win,slope,gamal,gamah;
unsigned char file_name[14], ch,ch1,file_name1[14];
clrscr() ;
printf("Freq.- response can be calculated using standard\n");
printf(" functions available to you if your reply to the\n");
printf(" following question is negative.\n");
printf("Is freq. response provided in a file? (y or n)-->");
while(((ch1=getche())!='y')&&(ch1!='n'));
switch(ch1)
{
case 'n':
case 'N':
printf("\nEnter # of points to be generated (e.g. 32x32)-->");
scanf("%d%c%d",&N1,&ch,&N1);
break;
case 'y' :
case 'Y ':
printf("\nEnter name of file storing magnitude response-->");
scanf("%s",file_name1);
fptr=fopen(file_name1,"r");
fscanf(fptr,"%d %d ",&N1,&N1);
}
N=N1>>1;
yt=wherey();
again :
gotoxy(1,yt);
printf(" ");
gotoxy(1,yt);
printf("Enter file name for storing impulse response -->");
scanf("%s",file_name);
170
if(((stricmp("FFT.DAT",file_name))==0)||
((stricmp("TEMP.DAT",file_name))==0)||
((stricmp("IFFT.DAT",file_name))==0))
{
printf("This is a reserved file name. Use some other name.");
goto again;
}
gotoxy(1,yt);
printf( " ");
ind=access(file_name,0);
while(!ind)
{
gotoxy(1,yt+1);
printf ( " ");
gotoxy(1,yt+1);
printf("File exists. Wish to overwrite? (y or n)-->");
while(((ch=tolower(getch()))!='y')&&(ch!='n'));
putch(ch);
switch(ch)
{
case 'y' :
ind=1 ;
break ;
case 'n' :
gotoxy(1,yt+1);
printf ( " ");
printf ( " ");
gotoxy(1,yt);
gotoxy(1,yt);
printf("Enter file name -->");
scanf("%s",file_name);
ind=access(file_name,0);
}
}
fptri=fopen("FFT.DAT","wb+");
fptro=fopen("IFFT.DAT","wb+");
buffi=(float *)malloc((N1<<1)*sizeof(float));
switch(ch1)
{
case 'n' :
case 'N' :
/*Generating data for IFFT. */
printf
("\nEnter cut-off freq. in rad./sec. ( cut-off <= n.)->");
171
scanf("%f",&Do);
Do=(Do/pi)*(float)N;
printf("Enter choice (1,2,3,4 or 5):\n");
printf(" 1. Low-pass with abrupt transition.\n");
printf(" 2. High-pass with abrupt transition.\n");
printf(" 3. Low-pass Butterworth.\n");
printf(" 4. High-pass Butterworth.\n");
printf(" 5. Homomorphic characteristcs--->");
while(((ch=getche())!='l')&&(ch!='2')
&&(ch!=13)&&(ch!='4')&&(ch!='5'));
if(ch=='5')
{
printf("\nEnter gain at không freq.-->");
scanf("%f",&gamah);
printf("Enter gain at high frequencies.-->");
scanf("%f",&gamah);
slope=(gamah-gamal)/Do;
}
Do*=Do;
if((ch=='3')||(ch=='4'))
{
printf("\n Enter order of Butterworth-->");
scanf("%d",&norder);
for(i=0;i<norder;i++)
Do*=Do;
}
printf("\n");
for(i=0;i<N1;i++)
{
ii=(i-N)*(i-N);
for(j=0;j<N1;j++)
{
R2=(float)((j-N)*(j-N)+ii);
if((ch=='3')||(ch=='4'))
{
for(j1=0;j1<norder;j1++)
R2*=R2;
}
switch(ch)
{
case '1': /* low-pass abrupt transition.*/
if(R2<Do) H=(float)1.0;
else H=(float)0.0;
break;
172
case '2': /* high-pass abrupt transition*/
if(R2<Do) H=(float)0.0;
else H=(float)1.0;
break;
case '3': /* low-order Butterworth */
H=Do/(Do+0.414*R2);
break;
case '4': /* high-order Butterworth. */
H=R2/(R2+0.414-Do);
break;
case '5': /* homomorphic characteristics. */
if ( R2<Do)
H=slope*sqrt((double)R2)+gamal;
else
H=gamah ;
}
jj=j<<1;
buffi[jj]=H;
buffi[jj+1]=0.0;
}
fwrite(buffi,N1,2*sizeof(float),fptri);
}
break;
case 'y':
case 'Y' :
for(i=0;i<N1;i++)
{
for(j=0;j<N1;j++)
{
fscanf(fptr,"%f ",&H);
jj=j <<1 ;
buffi[jj]=H;
buffi[jj+1]=(float)0.0;
}
fwrite(buffi,N1,2*sizeof(float),fptri);
}
fclose(fptr);
}
rewind(fptri);
m=(int)(log10((double)N1)/log10((double)2));
/* Allocating memory for bit reversal LUT. */
L=(unsigned int *)malloc(N1*sizeof(unsigned int));
173
/* Generate Look-up table for bit reversal. */
bit_reversal(L,m,N1);
/* Allocating memory for twiddle factors. */
n2=N1-1;
wr=(float *)malloc(n2*sizeof(float));
wi=(float *)malloc(n2*sizeof(float));
/* Generating twiddle factor. */
WTS(wr,wi,N1,1);
clrscr() ;
FFT2D(fptri,fptro,wr,wi,L,N1,m,1);
remove("FFT.DAT");
clrscr() ;
fptri=fopen("IFFT.DAT","rb");
fptro=fopen(file_name,"wb+");
nsq=(float)(N1*N1);
buffo=(float *)malloc(N1*sizeof(float));
for(i=0;i<N1;i++)
{
fread(buffi,N1,2*sizeof(float),fptri);
for(j=0;j<N1;j++)
buffo[j]=buffi[2*j]/nsq;
fwrite(buffo,N1,sizeof(float),fptro);
}
fclose(fptri);
remove("IFFT.DAT");
rewind(fptro);
printf("Enter selection of window function:");
printf("\n 1.Rectangular.");
printf("\n 2.Hann.");
printf("\n 3.Hamming.");
printf("\n 4.Blackmann.");
printf("\n 5.Kaiser.");
printf("\nEnter (1,2 ... etc.)-->");
while(((ch=getche())!='I')&&(ch!='2')&&
(ch!='3')&&(ch!='4')&&(ch!='5'));
printf("\n");
if(ch!='1')
{
N2=N1/2 ;
theta=pi/((float)N2*sqrt((double)2.0));
for(i=0;i<N1;i++)
{