Listing for AFMSPRT.C

/* AFMSPRT.C */

/* This program attempts finds the Fourier spectrum for a given

set of parameters of the AFM algorithm.

*/

#include <math.h>

#include <stdio.h>

#include <malloc.h>

#include <gl.h>

#include </usr/people/postgd/gansl/dsp/tek4014.h>

void main()

{

#define BesselOrderMin -14

#define BesselOrderMax 14

#define BesselIndexMin 0.0

#define BesselIndexMax 15.0

#define BesselIndexInc 0.01

/* 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

#define FINDMIN(a,b) ((a < b) ? a:b)

#define FINDMAX(a,b) ((a > b) ? a:b)

#define MinN -12

#define MaxN 12

#define MaxFreqAllowed 50000

#define MaxFreqToPlot 10000

#define FreqCarrier 440

#define FreqModulator 440

#define MinIndex 0.7

#define MaxIndex 2.5

#define IndexInc 0.05

#define MinRatio 0.6

#define MaxRatio 1.4

#define RatioInc 0.05

double jni,temp0,temp1,temp2;

double index,ratio;

double ampA,maxAmp,oldMaxAmp=1.0;

long freqA,freq,minFreq,maxFreq;

long freqCarrier,freqModulator,topFreq;

int runningN,minN,maxN;

FILE *inFile1,*outFile1;

int file1Open=0;

float *calcHar,*oldHar;

float temp1Float,temp2Float;

double tempDouble;

double coordinate[4];

int x1,y1,x2,y2,xMax,yMax;

float bessel1[29][1501];

int besselOrder;

float besselIndex;

outFile1 = fopen("AFMtest.dat","w");

fprintf(outFile1,"FreqCarrier = %ld Hz FreqModulator = %ld Hz\n",FreqCarrier,FreqModulator);

fprintf(outFile1,"Index ratio Freq Amplitude\n");

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;

freqCarrier = FreqCarrier;

freqModulator = FreqModulator;

minN = MinN;

maxN = MaxN;

/* 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 (index = MinIndex; index <= MaxIndex; index = index + IndexInc)

{

for (ratio = MinRatio; ratio <= MaxRatio; ratio = ratio + RatioInc)

{

calcHar = (float *) calloc(MaxFreqAllowed,sizeof(float));

for (runningN = minN; runningN <= maxN; runningN = runningN + 1)

{

freqA = freqCarrier + (runningN * freqModulator);

jni = bessel1[runningN+14][(int) (index/BesselIndexInc)];

ampA = jni * pow(ratio,(double) (runningN));

minFreq = FINDMIN(minFreq,freqA);

maxFreq = FINDMAX(maxFreq,freq);

if (freqA < 0)

{

freqA = -freqA;

ampA = -ampA;

}

*(calcHar + freqA) = *(calcHar + freqA) + (float) ampA;

}

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;

}

}

/*

printf("FreqCarrier = %ld Index = %.2f ratio = %.2f\n",freqCarrier,index,ratio);

*/

/* Start of IRIS codes for Console Graphics */

/*

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 for Console Graphics */

/* Start of IRIS codes for TEKTRONIX terminal */

/*

printf("%d %.2f %.2f",freqCarrier,index,ratio);

*/

for (freq = 0; freq <= MaxFreqToPlot; freq++)

{

temp1Float = *(calcHar + freq);

temp2Float = *(oldHar + freq);

*(oldHar + freq) = temp1Float;

if (temp1Float != 0.0)

{

fprintf(outFile1,"%.2f %.2f %ld %le\n",index,ratio,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);

}

}

/* End of IRIS codes for TEKTRONIX terminal */

}

else

{

printf("Bessel1.bin not successfully read into memory : file1Open = %d\n",file1Open);

printf("Program Aborted !\n");

}

}