Listing for DFMSPRT2.C
/* DFMSPRT2.C */
/* This program attempts finds the Fourier spectrum for a given
set of parameters of the DFM algorithm.
*/
#include <math.h>
#include <stdio.h>
#include <malloc.h>
/* DOS Extender includes */
/*
#include <graphics.h>
*/
/* End of DOS Extender includes */
/* Start of IRIS includes */
#include <gl.h>
#include </usr/people/postgd/gansl/dsp/tek4014.h>
/* End of IRIS includes */
/* Start of fabs() for DOS Extender, as fabs() is not defined in its library */
/*
double fabs(double value)
{
if (value < 0.0)
{
value = -value;
}
return(value);
}
*/
/* End of fabs() for DOS Extender */
void main()
{
#define BesselOrderMin -14
#define BesselOrderMax 14
#define BesselIndexMin 0.0
#define BesselIndexMax 15.0
#define BesselIndexInc 0.01
/* Start of IRIS definitions */
/* WindowX = 0 to 1279, WindowY = 0 to 1023 */
#define WindowXMin 10
#define WindowXMax 1000
#define WindowYMin 10
#define WindowYMax 800
/* TEKWindowX = 0 to 4095, TEKWindowY = 0 to 2500 */
#define TekWindowXMin 0
#define TekWindowXMax 4095
#define TekWindowYMin 0
#define TekWindowYMax 3119
#define TopPixel 2500
/* End of IRIS definition */
/* Start of DOS Extender definitions */
/* 320x200 resolution, 256 colors (0 to 255) */
/*
#define XMin 0
#define XMax 319
#define YMin 0
#define YMax 199
#define TopPixel 180
*/
/* End of DOS Extender definitions */
#define FINDMIN(a,b) ((a < b) ? a:b)
#define FINDMAX(a,b) ((a > b) ? a:b)
#define MinEvenN -12
#define MaxEvenN 12
#define MaxFreqAllowed 50000
#define MaxFreqToPlot 5000
#define Freq1 440
#define Freq2 660
#define MinIndex1 0.5
#define MaxIndex1 5.0
#define Index1Inc 0.1
#define MinIndex2 0.5
#define MaxIndex2 3.5
#define Index2Inc 0.02
double j0I2,j0I1,j,temp0,temp1,temp2,temp3,temp4;
double index1,index2;
double ampA,ampB,ampC,ampD,ampE,ampF,maxAmp,oldMaxAmp=1.0;
long freq,freqA,freqB,freqC,freqD,freqE,freqF,minFreq,maxFreq;
long freq1,freq2,oddFreq1,oddFreq2,evenFreq1,evenFreq2,topFreq;
int evenN,oddN,minOddN,maxOddN,minEvenN,maxEvenN;
FILE *inFile1,*outFile1;
int file1Open=0;
float *calcHar,*oldHar;
float temp1Float,temp2Float;
double tempDouble;
double coordinate[4];
int x1,y1,x2,y2,xMax,yMax,dosColor;
float bessel1[29][1501];
int besselOrder;
float besselIndex;
/*
outFile1 = fopen("spp02f11.dfm","w");
*/
if ((inFile1 = fopen("bessel1.bin","rb")) == NULL)
{
printf("Cannot open file bessel1.bin for reading\n");
printf("Aborting Now !!\a\n\n");
}
else
{
printf("File bessel1.bin successfully opened for reading\n\n");
file1Open = 1;
}
if (file1Open == 1)
{
oldHar = (float *) calloc(MaxFreqAllowed,sizeof(float));
printf("Reading in bessel1 function table\n\n");
for (besselOrder = BesselOrderMin; besselOrder <= BesselOrderMax; besselOrder++)
{
printf("Reading in Bessel1 function order %d\n",besselOrder);
for (besselIndex = BesselIndexMin; besselIndex <= BesselIndexMax; besselIndex = besselIndex + BesselIndexInc)
{
fread(&temp1Float,sizeof(float),1,inFile1);
bessel1[besselOrder+14][(int) (besselIndex/BesselIndexInc)] = temp1Float;
}
}
fclose(inFile1);
printf("\nReading in of bessel1 function table completed\n");
minFreq = 0;
maxFreq = 0;
freq1 = Freq1;
freq2 = Freq2;
minEvenN = MinEvenN;
maxEvenN = MaxEvenN;
minOddN = minEvenN + 1;
maxOddN = maxEvenN + 1;
/* Start of DOS Extender codes */
/*
xMax = XMax;
yMax = YMax;
GrInit();
GrSetColor(0,0,0,0);
*/
/* End of DOS Extender codes */
/* Start of IRIS codes for SGI */
prefposition(WindowXMin,WindowXMax,WindowYMin,WindowYMax);
winopen("Live Spectrum");
ortho2(0,MaxFreqToPlot,0,1.2);
color(BLACK);
clear();
printf("\n");
/* End of IRIS codes for SGI */
/* Start of IRIS codes for TEKTRONIX terminal */
/*
TEKinit();
*/
/* End of IRIS codes for TEKTRONIX terminal */
for (index1 = MinIndex1; index1 <= MaxIndex1; index1 = index1 + Index1Inc)
{
j0I1 = bessel1[0+14][(int) (index1/BesselIndexInc)];
for (index2 = MinIndex2; index2 <= MaxIndex2; index2 = index2 + Index2Inc)
{
calcHar = (float *) calloc(MaxFreqAllowed,sizeof(float));
j0I2 = bessel1[0+14][(int) (index2/BesselIndexInc)];
for (oddN = minOddN; oddN <= maxOddN; oddN = oddN + 2)
{
freqA = oddN * freq1;
ampA = j0I2 * bessel1[oddN+14][(int) (index1/BesselIndexInc)];
minFreq = FINDMIN(minFreq,freqA);
maxFreq = FINDMAX(maxFreq,freqA);
freqB = oddN * freq2;
ampB = j0I1 * bessel1[oddN+14][(int) (index2/BesselIndexInc)];
minFreq = FINDMIN(minFreq,freqB);
maxFreq = FINDMAX(maxFreq,freqB);
/* */
if (freqA < 0)
{
freqA = -freqA;
ampA = -ampA;
}
if (freqB < 0)
{
freqB = -freqB;
ampB = -ampB;
}
*(calcHar + freqA) = *(calcHar + freqA) + (float) ampA;
*(calcHar + freqB) = *(calcHar + freqB) + (float) ampB;
for (evenN = minEvenN; evenN <= maxEvenN; evenN = evenN + 2)
{
oddFreq1 = oddN * freq1;
oddFreq2 = oddN * freq2;
evenFreq1 = evenN * freq1;
evenFreq2 = evenN * freq2;
freqC = oddFreq1 + evenFreq2;
ampC = bessel1[oddN+14][(int) (index1/BesselIndexInc)] * bessel1[evenN+14][(int) (index2/BesselIndexInc)];
minFreq = FINDMIN(minFreq,freqC);
maxFreq = FINDMAX(maxFreq,freqC);
freqD = oddFreq1 - evenFreq2;
ampD = ampC;
minFreq = FINDMIN(minFreq,freqD);
maxFreq = FINDMAX(maxFreq,freqD);
freqE = oddFreq2 + evenFreq1;
ampE = bessel1[oddN+14][(int) (index2/BesselIndexInc)] * bessel1[evenN+14][(int) (index1/BesselIndexInc)];
minFreq = FINDMIN(minFreq,freqE);
maxFreq = FINDMAX(maxFreq,freqE);
freqF = oddFreq2 - evenFreq1;
ampF = ampE;
minFreq = FINDMIN(minFreq,freqF);
maxFreq = FINDMAX(maxFreq,freqF);
/* */
if (freqC < 0)
{
freqC = -freqC;
ampC = -ampC;
}
if (freqD < 0)
{
freqD = -freqD;
ampD = -ampD;
}
if (freqE < 0)
{
freqE = -freqE;
ampE = -ampE;
}
if (freqF < 0)
{
freqF = -freqF;
ampF = -ampF;
}
*(calcHar + freqC) = *(calcHar + freqC) + (float) ampC;
*(calcHar + freqD) = *(calcHar + freqD) + (float) ampD;
*(calcHar + freqE) = *(calcHar + freqE) + (float) ampE;
*(calcHar + freqF) = *(calcHar + freqF) + (float) ampF;
}
}
maxFreq = FINDMAX(maxFreq,-minFreq);
maxAmp = 0.0;
topFreq = 0;
for (freq = 0; freq <= MaxFreqToPlot; freq++)
{
temp1Float = *(calcHar + freq);
if (maxAmp < fabs((double) temp1Float))
{
maxAmp = fabs((double) temp1Float);
}
if (temp1Float != 0.0)
{
topFreq = freq;
}
}
/* Start of DOS Extender codes */
/*
printf("Freq1 = %ld Freq2 = %ld\n",freq1,freq2);
printf("Index1 = %.2f Index2 = %.2f\n",index1,index2);
outFile1 = fopen("dfmspec.dta","wt");
fprintf(outFile,"%ld\n%.2f\n%ld\n%.2f\n\n",freq1,index1,freq2,index2);
*/
/*
GrInit();
printf("Freq1 = %ld Freq2 = %ld\n",freq1,freq2);
printf("Index1 = %.2f Index2 = %.2f\n",index1,index2);
*/
/*
for (freq = 0; freq <= MaxFreqToPlot; freq++)
{
temp1Float = *(calcHar + freq);
temp2Float = *(oldHar + freq);
*(oldHar + freq) = temp1Float;
if (temp1Float != 0.0)
{
coordinate[1] = fabs((double) temp2Float)/oldMaxAmp;
x1 = (int) (((float) freq / (float) topFreq) * (float) xMax);
y1 = yMax - 1;
x2 = x1;
y2 = yMax - 1 - coordinate[1] * TopPixel;
GrLine(x1,y1,x2,y2,0);
if (temp1Float > 0.0)
{
dosColor = 14;
}
else
{
dosColor = 9;
}
coordinate[1] = fabs((double) temp1Float)/maxAmp;
x1 = (int) (((float) freq / (float) topFreq) * (float) xMax);
y1 = yMax - 1;
x2 = x1;
y2 = yMax - 1 - coordinate[1] * TopPixel;
GrLine(x1,y1,x2,y2,dosColor);
}
}
*/
/*
GrLine(0,yMax,xMax,yMax);
*/
/*
oldMaxAmp = maxAmp;
free(calcHar);
}
}
*/
/* End of DOS Extender codes */
/* Start of IRIS codes for SGI */
printf("Freq1 = %ld Index1 = %.2f Freq2 = %ld Index2 = %.2f\n",freq1,index1,freq2,index2);
for (freq = 0; freq <= MaxFreqToPlot; freq++)
{
temp1Float = *(calcHar + freq);
temp2Float = *(oldHar + freq);
*(oldHar + freq) = temp1Float;
if (temp1Float != 0.0)
{
bgnline();
color(BLACK);
coordinate[0] = (double) freq;
coordinate[1] = 0.0;
v2d(coordinate);
coordinate[1] = fabs((double) temp2Float)/oldMaxAmp;
v2d(coordinate);
endline();
bgnline();
if (temp1Float > 0.0)
{
color(YELLOW);
}
else
{
color(BLUE);
}
coordinate[0] = (double) freq;
coordinate[1] = 0.0;
v2d(coordinate);
coordinate[1] = fabs((double) temp1Float)/maxAmp;
v2d(coordinate);
endline();
}
}
oldMaxAmp = maxAmp;
free(calcHar);
}
}
/* End of IRIS codes SGI */
/* Start of IRIS codes for TEKTRONIX terminal */
/*
printf("%d %.2f %d %.2f",freq1,index1,freq2,index2);
*/
/*
for (freq = 0; freq <= MaxFreqToPlot; freq++)
{
temp1Float = *(calcHar + freq);
temp2Float = *(oldHar + freq);
*(oldHar + freq) = temp1Float;
if (temp1Float != 0.0)
{
*/
/*
fprintf(outFile1,"%.2f %ld %le\n",index2,freq,fabs(temp1Float));
*/
/*
TEKcolorIndex(0);
coordinate[1] = fabs((double) temp2Float)/oldMaxAmp;
x1 = (int) (((float) freq / (float) topFreq) * (float) TekWindowXMax);
y1 = TekWindowYMin + 1;
x2 = x1;
y2 = coordinate[1] * TopPixel + 1;
TEKdrawLine(x1,y1,x2,y2);
if (temp1Float > 0.0)
{
TEKcolorIndex(10);
}
else
{
TEKcolorIndex(9);
}
coordinate[1] = fabs((double) temp1Float)/maxAmp;
x1 = (int) (((float) freq / (float) topFreq) * (float) TekWindowXMax);
y1 = TekWindowYMin + 1;
x2 = x1;
y2 = coordinate[1] * TopPixel + 1;
TEKdrawLine(x1,y1,x2,y2);
}
}
oldMaxAmp = maxAmp;
free(calcHar);
}
}
TEKtoVT();
*/
/* End of IRIS codes for TEKTRONIX terminal */
}
ringbell();
sleep();
}