DECLARE SUB FFTCalc (xreal!(), yimag!(), numdat%)

DECLARE FUNCTION FFTMagnitude (xr!(), yi!(), ns%, i%)
DECLARE FUNCTION FFTPhase (xr!(), yi!(), ns%, i%)

'****************************************************************************
'Module SPECTRUM.BAS
'These subroutines do spectral analysis on data samples.  They assume
'synchronous data, so a complete number of cycles is present in the data
'and there are no end effects.  No windowing is done.
'
'Zachary Wolf
'9/18/95
'****************************************************************************

SUB spcalcfft (ns%, samp!(), fftm!(), fftp!())

'****************************************************************************
'This subroutine calculates the FFT magnitude and phase of a set of samples.
'
'Input:
' ns%, the number of samples, power of 2
' samp!(0 to ns%-1), the samples
'
'Output:
' fftm!(0 to ns%/2), the magnitude of the FFT coefficients
' fftp!(0 to ns%/2), the phase of the FFT coeficients in degrees, -180 to 180
'
'Zachary Wolf
'9/12/95
'****************************************************************************

'Initialize
DIM xr!(0 TO ns% - 1)
DIM yi!(0 TO ns% - 1)
FOR i% = 0 TO ns% - 1
 xr!(i%) = samp!(i%)
 yi!(i%) = 0!
NEXT i%

'Get the frequency spectrum using an FFT
CALL FFTCalc(xr!(), yi!(), ns%)

'Calculate the magnitude and phase of each term in the spectrum
FOR i% = 0 TO ns% / 2
 fftm!(i%) = FFTMagnitude(xr!(), yi!(), ns%, i%)
 fftp!(i%) = FFTPhase(xr!(), yi!(), ns%, i%)
NEXT i%

'Convert from radians to degrees
'Convert from (0,360) to (-180,180)
FOR i% = 0 TO ns% / 2
 fftp!(i%) = fftp!(i%) * 180! / 3.141592654#
 IF fftp!(i%) > 180! THEN fftp!(i%) = -(360! - fftp!(i%))
NEXT i%

END SUB

SUB spcalcfit (ncy%, ns%, samp!(), a!, b!, c!, phi!, res!)

'****************************************************************************
'This subroutine fits ns% data points to a sine wave having n cycles in the
'ns points.
'The fit is to y = a!*sin(n*2pi*i/ns) + b!*cos(n*2pi*i/ns)
'                = c!*cos(n*2pi*i/ns + phi!)
'The rms of the residuals is also returned.
'The conversion from a and b to c and phi goes as follows:
'y = c!*cos(n*2pi*i/ns - nphi!)
'nphi! = -phi!
'y = c!*cos(n*2pi*i/ns)*cos(nphi!) + c!*sin(n*2pi*i/ns)*sin(nphi!)
'a! = c!*sin(nphi!), b! = c!*cos(nphi!)
'c!^2 = a!^2 + b!^2, c! = sqr(a!^2 + b!^2)
'a!/b! = tan(nphi!), nphi! = atan(a!/b!), consider different quadrants
'
'Input:
' ncy%, the number of cycles in the fitted sine wave
' ns%, the number of data points
' samp!(0 to ns%-1), the data points
'
'Output:
' a!, the coeffecient of the sin term from the fit
' b!,            "           cos        "
' c!, the fitted amplitude
' phi!, the fitted phase in degrees
' res!, the rms of the residuals
'
'Zachary Wolf
'8/19/94, 5/4/95
'****************************************************************************

'Parameters
pi! = 3.1415926#

'Compute various sums
Sss! = 0!
Scc! = 0!
Ssc! = 0!
Sys! = 0!
Syc! = 0!
FOR i% = 0 TO ns% - 1
 Sss! = Sss! + SIN(ncy% * 2 * pi! * i% / ns%) * SIN(ncy% * 2 * pi! * i% / ns%)
 Scc! = Scc! + COS(ncy% * 2 * pi! * i% / ns%) * COS(ncy% * 2 * pi! * i% / ns%)
 Ssc! = Ssc! + SIN(ncy% * 2 * pi! * i% / ns%) * COS(ncy% * 2 * pi! * i% / ns%)
 Sys! = Sys! + samp!(i%) * SIN(ncy% * 2 * pi! * i% / ns%)
 Syc! = Syc! + samp!(i%) * COS(ncy% * 2 * pi! * i% / ns%)
NEXT i%

'Compute the determinant
d! = Sss! * Scc! - Ssc! * Ssc!

'Compute the sin and cos coefficients
a! = (Sys! * Scc! - Ssc! * Syc!) / d!
b! = (Sss! * Syc! - Sys! * Ssc!) / d!

'Compute the amplitude coefficient
c! = SQR(a! ^ 2 + b! ^ 2)

'Compute the phase coefficient, nphi! = -phi!
IF a! > 0 AND b! > 0 THEN nphi! = ATN(a! / b!)
IF a! > 0 AND b! < 0 THEN nphi! = pi! - ATN(-a! / b!)
IF a! < 0 AND b! > 0 THEN nphi! = -ATN(-a! / b!)
IF a! < 0 AND b! < 0 THEN nphi! = -pi! + ATN(a! / b!)
IF a! = 0 AND b! >= 0 THEN nphi! = 0!
IF a! = 0 AND b! < 0 THEN nphi! = pi!
IF b! = 0 AND a! > 0 THEN nphi! = pi! / 2
IF b! = 0 AND a! < 0 THEN nphi! = -pi! / 2
phi! = -nphi!
phi! = phi! * 180! / pi!        'convert to degrees

'Compute the rms of the residuals
res! = 0!
FOR i% = 0 TO ns% - 1
 res! = res! + (samp!(i%) - a! * SIN(ncy% * 2 * pi! * i% / ns%) - b! * COS(ncy% * 2 * pi! * i% / ns%)) ^ 2
NEXT i%
res! = SQR(res! / ns%)

END SUB

SUB spcalchar (ns%, fftm!(), fftp!(), nfun%, nh%, harm!(), harp!())

'****************************************************************************
'This subroutine returns the harmonics of a given input fundamental frequency.
'It assumes nfun% complete cycles of the fundamental frequency in the
'original data record.  The FFT coefficients of this frequency and
'its harmonics is returned.
'
'Input:
' ns%, the number of samples of the original data
' fftm!(0 to ns%/2), fft magnitudes
' fftp!(0 to ns%/2), fft phases in degrees
' nfun%, the number of cycles of the fundamental frequency in the data record
' nh%, the number of harmonics of the fundamental frequency to return
'
'Output:
' harm!(1 to nh%), magnitude at each harmonic
' harp!(1 to nh%), phase in degrees at each harmonic
'
'Zachary Wolf
'9/15/95
'****************************************************************************

'Return the requested values from the input fft values
FOR i% = 1 TO nh%
 IF i% * nfun% < ns% / 2 THEN
  harm!(i%) = fftm!(i% * nfun%)
  harp!(i%) = fftp!(i% * nfun%)
 ELSE
  harm!(i%) = 0!
  harp!(i%) = 0!
 END IF
NEXT i%

END SUB

SUB spgensample (ns%, nfun%, nmain%, i%, samp!)

'****************************************************************************
'This subroutine generates data samples for testing.
'
'Input:
' ns%, the number of samples in a complete record
' nfun%, the number of fundamental oscillation cycles in the data record
' nmain%, the number of the main harmonic
' i%, the sample number to generate
'
'Output:
' samp!, the data sample
'
'Zachary Wolf
'6/17/96
'****************************************************************************

'Define parameters
pi! = 3.141592654#
phi! = pi! / 4!

'Generate the sample
samp! = 1! * COS(nfun% * nmain% * 2! * pi! * i% / ns% + phi!)
samp! = samp! + .005
samp! = samp! + .001 * COS(1! * nfun% * 2! * pi! * i% / ns% + phi!)
samp! = samp! + .002 * COS(2! * nfun% * 2! * pi! * i% / ns% + phi!)
samp! = samp! + .003 * COS(3! * nfun% * 2! * pi! * i% / ns% + phi!)
samp! = samp! + .004 * COS(4! * nfun% * 2! * pi! * i% / ns% + phi!)
samp! = samp! + .005 * COS(5! * nfun% * 2! * pi! * i% / ns% + phi!)
samp! = samp! + .006 * COS(6! * nfun% * 2! * pi! * i% / ns% + phi!)
samp! = samp! + .007 * COS(7! * nfun% * 2! * pi! * i% / ns% + phi!)
samp! = samp! + .008 * COS(8! * nfun% * 2! * pi! * i% / ns% + phi!)
samp! = samp! + .009 * COS(9! * nfun% * 2! * pi! * i% / ns% + phi!)
samp! = samp! + .01 * COS(10! * nfun% * 2! * pi! * i% / ns% + phi!)
samp! = samp! + .011 * COS(11! * nfun% * 2! * pi! * i% / ns% + phi!)
samp! = samp! + .012 * COS(12! * nfun% * 2! * pi! * i% / ns% + phi!)
samp! = samp! + .013 * COS(13! * nfun% * 2! * pi! * i% / ns% + phi!)
samp! = samp! + .014 * COS(14! * nfun% * 2! * pi! * i% / ns% + phi!)
samp! = samp! + .015 * COS(15! * nfun% * 2! * pi! * i% / ns% + phi!)

END SUB

SUB spgensamples (ns%, nfun%, nmain%, samp!())

'****************************************************************************
'This subroutine generates data samples for testing.
'
'Input:
' ns%, the number of samples to generate, power of 2
' nfun%, the number of fundamental oscillation cycles in the data record, power of 2
' nmain%, the number of the main harmonic
'
'Output:
' samp!(0 to ns%-1), the data samples
'
'Zachary Wolf
'9/18/95
'****************************************************************************

'Define parameters
pi! = 3.141592654#
phi! = pi! / 4!

'Message
PRINT "SPGENSAMPLES: generating test data"

'Generate the samples
FOR i% = 0 TO ns% - 1
 samp!(i%) = 1! * COS(nfun% * nmain% * 2! * pi! * i% / ns% + phi!)
 samp!(i%) = samp!(i%) + .005
 samp!(i%) = samp!(i%) + .001 * COS(1! * nfun% * 2! * pi! * i% / ns% + phi!)
 samp!(i%) = samp!(i%) + .002 * COS(2! * nfun% * 2! * pi! * i% / ns% + phi!)
 samp!(i%) = samp!(i%) + .003 * COS(3! * nfun% * 2! * pi! * i% / ns% + phi!)
 samp!(i%) = samp!(i%) + .004 * COS(4! * nfun% * 2! * pi! * i% / ns% + phi!)
 samp!(i%) = samp!(i%) + .005 * COS(5! * nfun% * 2! * pi! * i% / ns% + phi!)
 samp!(i%) = samp!(i%) + .006 * COS(6! * nfun% * 2! * pi! * i% / ns% + phi!)
 samp!(i%) = samp!(i%) + .007 * COS(7! * nfun% * 2! * pi! * i% / ns% + phi!)
 samp!(i%) = samp!(i%) + .008 * COS(8! * nfun% * 2! * pi! * i% / ns% + phi!)
 samp!(i%) = samp!(i%) + .009 * COS(9! * nfun% * 2! * pi! * i% / ns% + phi!)
 samp!(i%) = samp!(i%) + .01 * COS(10! * nfun% * 2! * pi! * i% / ns% + phi!)
 samp!(i%) = samp!(i%) + .011 * COS(11! * nfun% * 2! * pi! * i% / ns% + phi!)
 samp!(i%) = samp!(i%) + .012 * COS(12! * nfun% * 2! * pi! * i% / ns% + phi!)
 samp!(i%) = samp!(i%) + .013 * COS(13! * nfun% * 2! * pi! * i% / ns% + phi!)
 samp!(i%) = samp!(i%) + .014 * COS(14! * nfun% * 2! * pi! * i% / ns% + phi!)
 samp!(i%) = samp!(i%) + .015 * COS(15! * nfun% * 2! * pi! * i% / ns% + phi!)
NEXT i%
GOTO done

'Another coil voltage generator
coil1:
'v = N vth Br L
'  = N R omega G R sin(2 theta) L
'  = N GL R^2 2 pi f sin(2 theta)
gl! = 1!
f0! = 2!
FOR i% = 0 TO ns% - 1
 f! = f0! * (1! + .01 * SIN(2! * pi! * i% / ns%))
 samp!(i%) = coilnturnsP% * gl! * P! ^ 2 * 2 * pi! * f! * SIN(2! * 2! * pi! * i% / ns%)
NEXT i%
GOTO done

'Another coil voltage generator
coil3:
'v = N vth Br L
'  = N R omega 1/2 S R^2 sin(3 theta) L
'  = N SL 1/2 R^3 2 pi f sin(3 theta)
sl! = 1!
f0! = 2!
FOR i% = 0 TO ns% - 1
 f! = f0!
' f! = f0! * (1! + .01 * SIN(2! * pi! * i% / ns%))
 samp!(i%) = coilnturnsP% * sl! * .5 * coilradiusmP! ^ 3 * 2 * pi! * f! * SIN(3! * 2! * pi! * i% / ns%)
NEXT i%
GOTO done

done:

END SUB

SUB spplotfft (ns%, fftm!(), fftp!())

'****************************************************************************
'This subroutine gets voltage samples from the coil and puts them in
'a file for plotting.
'
'Zachary Wolf
'5/4/95
'****************************************************************************

'Initialize a plot file
filenum% = FREEFILE
OPEN "fft.plt" FOR OUTPUT AS filenum%

'Write the samples to the plot file
FOR k% = 0 TO ns% / 2
 PRINT #filenum%, USING "####"; k%;
 PRINT #filenum%, USING " ###.########"; fftm!(k%);
 PRINT #filenum%, USING " ####.#####"; fftp!(k%)
NEXT k%

'Close the log file
CLOSE filenum%

END SUB

SUB spplothar (nh%, harm!(), harp!())

'****************************************************************************
'This subroutine prints the magnitudes and phases of the harmonics
'of a given frequency in the input data record.
'
'Input:
' nh%, the number of harmonics
' harm%(1 to nh%), the magnitude of each harmonic
' harp%(1 to nh%), the phase of each harmonic
'
'Zachary Wolf
' 9/15/95
'****************************************************************************

'Print the first 16 harmonics
IF nh% > 16 THEN
 nh16% = 16
ELSE
 nh16% = nh%
END IF

'Initialize a print file
filenum% = FREEFILE
OPEN "har.plt" FOR OUTPUT AS filenum%

'Write the samples to the print file
PRINT #filenum%, "i, mag(i), phase(i)"
FOR k% = 1 TO nh16%
 PRINT #filenum%, USING "####"; k%;
 PRINT #filenum%, USING "  ##.########"; harm!(k%);
 PRINT #filenum%, USING "  ####.###"; harp!(k%)
NEXT k%

'Close the log file
CLOSE filenum%

END SUB

SUB spplotsamp (ns%, samp!())

'****************************************************************************
'This subroutine gets voltage samples from the coil and puts them in
'a file for plotting.
'
'Zachary Wolf
'5/4/95
'****************************************************************************

'Initialize a plot file
filenum% = FREEFILE
OPEN "samp.plt" FOR OUTPUT AS filenum%

'Write the samples to the plot file
FOR k% = 0 TO ns% - 1
 PRINT #filenum%, USING "####"; k%;
 PRINT #filenum%, USING " ###.########"; samp!(k%)
NEXT k%

'Close the log file
CLOSE filenum%

END SUB

SUB spprntfit (ncy%, a!, b!, c!, phi!, res!)

'****************************************************************************
'This subroutine logs the parameters of the fit to the samples.
'
'Input:
' ncy%, the number of cycles in the fitted sine wave
' a!, the coeffecient of the sin term from the fit
' b!,            "           cos        "
' c!, the fitted amplitude
' phi!, the fitted phase in degrees
' res!, the rms of the residuals
'
'Zachary Wolf
'1/5/95, 5/4/95
'****************************************************************************

'Open the log file
filenum% = FREEFILE
OPEN "fit.dat" FOR OUTPUT AS filenum%

'Write the fit parameters
PRINT #filenum%,
PRINT #filenum%, TIME$; "  Data Fit: "
PRINT #filenum%, "Yfit(i) = ";
PRINT #filenum%, USING "##.#######"; c!;
PRINT #filenum%, USING " * cos(##*2pi*i/ns + "; ncy%;
PRINT #filenum%, USING "####.##"; phi!;
PRINT #filenum%, " deg)"

PRINT #filenum%, "        = ";
PRINT #filenum%, USING "##.#######"; a!;
PRINT #filenum%, USING " * sin(##*2pi*i/ns) + "; ncy%;
PRINT #filenum%, USING "##.#######"; b!;
PRINT #filenum%, USING " * cos(##*2pi*i/ns)"; ncy%

PRINT #filenum%, "RMS of the residuals = ";
PRINT #filenum%, USING "##.#######"; res!

'Close the file
CLOSE filenum%

END SUB