//
//
//   File: how2.cpp      Created 12-MAY-2007        Richard B. Johnson
//                                                  RichardBJohnson@comcast.net
//
//   The data, text, and source-code within these files are presumed to
//   be accurate and free of encumbrances. Use at your own risk.
//
//
//  This specifically shows how to convert a time-domain array of data
//  to a frequency-domain (spectrum analyzer) display.
//

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <fcntl.h>
#include <time.h>
#include <sys/stat.h>
#include <io.h>
#include <sys/timeb.h>
#include <atlstr.h>
#include <process.h>
#include <direct.h>

#include "IScope_Common.h"
#include "MeasureLibrary.h"
#include "altmath.h"
#include "how2.h"

#define NR_FILTER 0x1000

extern MSL msl;
extern SIG sig;
HOW how;
double HOW::srate = SAMPLE_RATE;

static struct HOW::file_dat fd[] = {
	{ "CLEAN-00-NOWIN",  &HOW::acquire_clean_data,  &SIG::MakeRectWIN,         0.0 },
    { "CLEAN+45-NOWIN",  &HOW::acquire_clean_data,  &SIG::MakeRectWIN,        45.0 },
    { "CLEAN+90-NOWIN",  &HOW::acquire_clean_data,  &SIG::MakeRectWIN,        90.0 },
    { "CLEAN-90-NOWIN",  &HOW::acquire_clean_data,  &SIG::MakeRectWIN,       -90.0 },
    { "CLEAN-45-NOWIN",  &HOW::acquire_clean_data,  &SIG::MakeRectWIN,       -45.0 },
    { "DIRTY-00-NOWIN",  &HOW::acquire_dirty_data,  &SIG::MakeRectWIN,         0.0 },
    { "DIRTY+45-NOWIN",  &HOW::acquire_dirty_data,  &SIG::MakeRectWIN,        45.0 },
    { "DIRTY+90-NOWIN",  &HOW::acquire_dirty_data,  &SIG::MakeRectWIN,        90.0 },
    { "DIRTY-90-NOWIN",  &HOW::acquire_dirty_data,  &SIG::MakeRectWIN,       -90.0 },
    { "DIRTY-45-NOWIN",  &HOW::acquire_dirty_data,  &SIG::MakeRectWIN,       -45.0 },
    { "NOISY-00-NOWIN",  &HOW::acquire_noisy_data,  &SIG::MakeRectWIN,         0.0 },
    { "NOISY+45-NOWIN",  &HOW::acquire_noisy_data,  &SIG::MakeRectWIN,        45.0 },
    { "NOISY+90-NOWIN",  &HOW::acquire_noisy_data,  &SIG::MakeRectWIN,        90.0 },
    { "NOISY-90-NOWIN",  &HOW::acquire_noisy_data,  &SIG::MakeRectWIN,       -90.0 },
    { "NOISY-45-NOWIN",  &HOW::acquire_noisy_data,  &SIG::MakeRectWIN,       -45.0 },
    { "DISTO-00-NOWIN",  &HOW::acquire_dist_data,   &SIG::MakeRectWIN,         0.0 },
    { "DISTO+45-NOWIN",  &HOW::acquire_dist_data,   &SIG::MakeRectWIN,        45.0 },
    { "DISTO+90-NOWIN",  &HOW::acquire_dist_data,   &SIG::MakeRectWIN,        90.0 },
    { "DISTO-90-NOWIN",  &HOW::acquire_dist_data,   &SIG::MakeRectWIN,       -90.0 },
    { "DISTO-45-NOWIN",  &HOW::acquire_dist_data,   &SIG::MakeRectWIN,       -45.0 },
    { "CLEAN-00-GAUSS",  &HOW::acquire_clean_data,  &SIG::MakeWeedonGauss,     0.0 },
    { "CLEAN+45-GAUSS",  &HOW::acquire_clean_data,  &SIG::MakeWeedonGauss,    45.0 },
    { "CLEAN+90-GAUSS",  &HOW::acquire_clean_data,  &SIG::MakeWeedonGauss,    90.0 },
    { "CLEAN-90-GAUSS",  &HOW::acquire_clean_data,  &SIG::MakeWeedonGauss,   -90.0 },
    { "CLEAN-45-GAUSS",  &HOW::acquire_clean_data,  &SIG::MakeWeedonGauss,   -45.0 },
    { "DIRTY-00-GAUSS",  &HOW::acquire_dirty_data,  &SIG::MakeWeedonGauss,     0.0 },
    { "DIRTY+45-GAUSS",  &HOW::acquire_dirty_data,  &SIG::MakeWeedonGauss,    45.0 },
    { "DIRTY+90-GAUSS",  &HOW::acquire_dirty_data,  &SIG::MakeWeedonGauss,    90.0 },
    { "DIRTY-90-GAUSS",  &HOW::acquire_dirty_data,  &SIG::MakeWeedonGauss,   -90.0 },
    { "DIRTY-45-GAUSS",  &HOW::acquire_dirty_data,  &SIG::MakeWeedonGauss,   -45.0 },
    { "NOISY-00-GAUSS",  &HOW::acquire_noisy_data,  &SIG::MakeWeedonGauss,     0.0 },
    { "NOISY+45-GAUSS",  &HOW::acquire_noisy_data,  &SIG::MakeWeedonGauss,    45.0 },
    { "NOISY+90-GAUSS",  &HOW::acquire_noisy_data,  &SIG::MakeWeedonGauss,    90.0 },
    { "NOISY-90-GAUSS",  &HOW::acquire_noisy_data,  &SIG::MakeWeedonGauss,   -90.0 },
    { "NOISY-45-GAUSS",  &HOW::acquire_noisy_data,  &SIG::MakeWeedonGauss,   -45.0 },
    { "DISTO-00-GAUSS",  &HOW::acquire_dist_data,   &SIG::MakeWeedonGauss,     0.0 },
    { "DISTO+45-GAUSS",  &HOW::acquire_dist_data,   &SIG::MakeWeedonGauss,    45.0 },
    { "DISTO+90-GAUSS",  &HOW::acquire_dist_data,   &SIG::MakeWeedonGauss,    90.0 },
    { "DISTO-90-GAUSS",  &HOW::acquire_dist_data,   &SIG::MakeWeedonGauss,   -90.0 },
    { "DISTO-45-GAUSS",  &HOW::acquire_dist_data,   &SIG::MakeWeedonGauss,   -45.0 },
    { "CLEAN-00-LANCZ",  &HOW::acquire_clean_data,  &SIG::MakeLanczos,         0.0 },
    { "CLEAN-00-BKHAR",  &HOW::acquire_clean_data,  &SIG::MakeBlackmanHarris,  0.0 },
	{ "CLEAN-00-BKMAN",  &HOW::acquire_clean_data,  &SIG::MakeBlackman,        0.0 },
    { "CLEAN+45-BKHAR",  &HOW::acquire_clean_data,  &SIG::MakeBlackmanHarris, 45.0 },
    { "CLEAN+90-BKHAR",  &HOW::acquire_clean_data,  &SIG::MakeBlackmanHarris, 90.0 },
    { "CLEAN-90-BKHAR",  &HOW::acquire_clean_data,  &SIG::MakeBlackmanHarris,-90.0 },
    { "CLEAN-45-BKHAR",  &HOW::acquire_clean_data,  &SIG::MakeBlackmanHarris,-45.0 },
    { "DIRTY-00-BKHAR",  &HOW::acquire_dirty_data,  &SIG::MakeBlackmanHarris,  0.0 },
    { "DIRTY+45-BKHAR",  &HOW::acquire_dirty_data,  &SIG::MakeBlackmanHarris, 45.0 },
    { "DIRTY+90-BKHAR",  &HOW::acquire_dirty_data,  &SIG::MakeBlackmanHarris, 90.0 },
    { "DIRTY-90-BKHAR",  &HOW::acquire_dirty_data,  &SIG::MakeBlackmanHarris,-90.0 },
    { "DIRTY-45-BKHAR",  &HOW::acquire_dirty_data,  &SIG::MakeBlackmanHarris,-45.0 },
    { "NOISY-00-BKHAR",  &HOW::acquire_noisy_data,  &SIG::MakeBlackmanHarris,  0.0 },
    { "NOISY+45-BKHAR",  &HOW::acquire_noisy_data,  &SIG::MakeBlackmanHarris, 45.0 },
    { "NOISY+90-BKHAR",  &HOW::acquire_noisy_data,  &SIG::MakeBlackmanHarris, 90.0 },
    { "NOISY-90-BKHAR",  &HOW::acquire_noisy_data,  &SIG::MakeBlackmanHarris,-90.0 },
    { "NOISY-45-BKHAR",  &HOW::acquire_noisy_data,  &SIG::MakeBlackmanHarris,-45.0 },
    { "DISTO-00-BKHAR",  &HOW::acquire_dist_data,   &SIG::MakeBlackmanHarris,  0.0 },
	{ "DIST0-00-BKMAN",  &HOW::acquire_dist_data,   &SIG::MakeBlackman,        0.0 },
    { "DISTO+45-BKHAR",  &HOW::acquire_dist_data,   &SIG::MakeBlackmanHarris, 45.0 },
    { "DISTO+90-BKHAR",  &HOW::acquire_dist_data,   &SIG::MakeBlackmanHarris, 90.0 },
    { "DISTO-90-BKHAR",  &HOW::acquire_dist_data,   &SIG::MakeBlackmanHarris,-90.0 },
    { "DISTO-45-BKHAR",  &HOW::acquire_dist_data,   &SIG::MakeBlackmanHarris,-45.0 },
    { "SQUAR-00-NOWIN",  &HOW::acquire_square_data, &SIG::MakeRectWIN,         0.0 },
    { "SQUAR-00-GAUSS",  &HOW::acquire_square_data, &SIG::MakeWeedonGauss,     0.0 },
    { "GLITCH-00-BKHAR", &HOW::acquire_glitch_data, &SIG::MakeBlackmanHarris,  0.0 }, 
    { "SQUAR-00-BKHAR",  &HOW::acquire_square_data, &SIG::MakeBlackmanHarris,  0.0} };


static const char Data[]="Data\\";

//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//  This gets some noisy floating-point data for tests.
//
//
void HOW::get_noisy_fpdat(double *buf, int len)
{
	short *sp;
	sp = new short [len];
    acquire_noisy_data(sp, 0.0, len);      // Acquire "previous" data
	sig.CvtShort2Double(buf, sp, len);     // Convert the data, put into "previous"
	delete[] sp;
}

#define ArraySize(a) (sizeof(a) / sizeof(a[0]))
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//
//  This tells you the actual frequency of the 'carrier' component.
//  So that we are using the same waveform, I will use a square wave
//  since it will work with our internal MeasureFrequency() procedure.
//
void HOW::show_frequency(double *buf, int len)
{
    int i, xing;
    int sample_rate;
    double f, dif;
    sample_rate = (len / NODES);
//
//  This square-wave alternates between -1 and +1 with zero time
//  in between. Not a practical waveform, but guranteed to have
//  NODES (currently 32) cycles in the waveform.
//
    for(i = 0; i< len; i++)       /* Make a square wave */
        buf[i] = ((i / (sample_rate / 2)) & 1) ? 1.0 : -1.0;
//
//  The frequency is the number of times per second a signal
//  crosses the  same axis in the same direction. So, let's
//  verify that the signal has 32 such crossings.
//
    f = 0.0;
    for(i = xing= 0; i< len; i++)       // Verify a square wave 
    {
        if(buf[i] > f)
        {
            f = buf[i];
            xing++;
        }
        else if(buf[i] < f)
            f = buf[i];
    }
    if(xing != NODES)
    {
        fprintf(stderr, "The square wave had only %d cycles in line %d\n",
                         xing, __LINE__);
        exit(EXIT_FAILURE);
    }
    printf("The square wave is verified to have %d cycles\n", xing);
//
//  Now measure and display the period. It must be 1/32 = 0.03125
//  or something is wrong.
//
    f = msl.MeasurePeriod(buf, len, (double)len, 1, 50.0);
    dif = f - 0.03125;
    if((dif = ((dif < 0) ? -dif:dif))> DYN)
    {
        fprintf(stderr, "Error from line %d, should have been %f, was %f\n",
               __LINE__, 0.03125, f);
        exit(EXIT_FAILURE);
    }
    printf("Period from MeasurePeriod() = %.5f\n", f);
//
//  Now get the frequency from the MeasureFrequency code.
//  It should be 32.0 or something is wrong.
//
    f = msl.MeasureFrequency(buf, len, (double)len, 9, 50.0);
    dif = f - 32.0;
    if((dif = ((dif < 0) ? -dif:dif))> DYN)
    {
        fprintf(stderr, "Error from line %d, should have been %f, was %f\n",
               __LINE__, 32.0, f);
        exit(EXIT_FAILURE);
    }
    printf("Frequency from MeasureFrequency = %.1f\n", f);
#if 0
    sig.MulConst(buf, 1e-4, len);
    f = msl.MeasureFrequency(buf, len, (double)len, 9, 50.0);
    dif = f - 32.0;
    if((dif = ((dif < 0) ? -dif:dif))> DYN)
    {
        fprintf(stderr, "Error from line %d, should have been %f, was %f\n",
               __LINE__, 32.0, f);
        exit(EXIT_FAILURE);
    }
    printf("Frequency from MeasureFrequency = %.1f\n", f);
#endif

//
//  Now Get the frequency by doing a FFT. This must return the same
//  value or something is wrong. It will never return the __exact__
//  same value because it is floating-point, but it should be very close.
//
    f = sig.do_FFT(buf, buf, SIG::WINDOW_TYPE_RCOS, len);
    dif = f - 32.0;
    if((dif = ((dif < 0) ? -dif:dif))> 1e-6)
    {
        fprintf(stderr, "Error from line %d, should have been %f, was %f\n",
               __LINE__, 32.0, f);
        exit(EXIT_FAILURE);
    }
    printf("Frequency from doFFT() = %.1f\n", f);
    printf("This is %.1f cycles in %d samples.\n", f, len);
    printf("With a sample rate of %.1f Hz, the frequency is %.1f Hz\n",
        srate, f / (srate / (double) len) );
}
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//
//   This makes a ring-down on each edge.
//
void HOW::make_ring_down(short *buf, int len)
{
    double theta;
    double sample_rate = (double) len;
    double val;
    int i;
//
//  Generate the fundamental component (5 x 2 * PI).
//
    for(i=0; i< len; i++)
    {
        theta  = ((double)i / sample_rate) * 10.0 * PI;
        val =  cos(theta) * 400.0 * (double) (len - i) / (double) len;
        buf[i] += (short) (val + PAD);
    }
}
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//
//  This pretends to acquire square-waves with ring-down glitches.
//
void HOW::acquire_glitch_data(short *buf, double phase, int len)
{
    int i;
    int sample_rate = len / NODES;
    phase = phase;
    for(i=0; i< len; i++)
        buf[i] = ((i / (sample_rate / 2)) & 1) ? 10000 : -10000;

    for(i=0; i< len; i++)
        if((i / (sample_rate / 2) & 1) && ((len - i) > (sample_rate/2)))
             make_ring_down(&buf[i], sample_rate / 4);
}
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//
//  This pretends to acquire noisy data. It generates data as short
//  (16-bit) integers.
//
void HOW::acquire_noisy_data(short *indat, double phase, int len)
{
    double theta;
    double sample_rate = (double) (len / NODES);
    int i;
	unsigned int munge;
	short rnd;
    phase = (phase / 360.0) * 2.0 * PI;
//
//  Generate the fundamental component (2 x PI).
//
    for(i=0; i< len; i++)
    {
        theta  = ((double)i / sample_rate) * 2.0 * PI;
        theta += phase;
        indat[i] = (short) ((cos(theta) * 16384.0) + PAD);
    }
//
//  Add some noise. This is not real good noise, but it'll do.
//

	while(rand() & 1)
		;          // Remove correlation

    for(i=0; i< len; i++)
	{
		munge = (unsigned int) rand();
		rnd = (short) munge & 0xffff;
		munge >>= 0x10;
		rnd |= (short) munge;
		rnd /= 2;

        indat[i] += rnd;
	}
}
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//
//  This pretends to acquire clean data. It generates data as short
//  (16-bit) integers. A virtual oscilloscopw would substitute its
//  data input routine here.
//
void HOW::acquire_clean_data(short *indat, double phase, int len)
{
    double theta;
    double sample_rate = (double) (len / NODES);
    int i;
    phase = (phase / 360.0) * 2.0 * PI;
/*
 *  Generate the fundamental component (2 x PI).
 */
    for(i=0; i< len; i++)
    {
        theta  = ((double)i / sample_rate) * 2.0 * PI;
        theta += phase;
        indat[i] = (short) ((cos(theta) * 32767.5) + PAD);
    }
}
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//
//  This pretends to acquire harmonic-riddled data. It generates data
//  as short (16-bit) integers. We are going to have 4 harmonics
//  that are 10 percent of the fundamental, i.e., -20 dB. They can
//  all sit on top of each other, in sync, so that leaves 100% - 40% =
//  60% available for the fundamental. The maximum fundamental signal
//  in the 16-bit world is 32767. Therefore, we will use 0.6 x 32767 =
//  19,660 for the fundamental.
//

void HOW::acquire_dirty_data(short *indat, double phase, int len)
{
    double theta;
    double sample_rate = (double) (len / NODES);
    int i;
    phase = (phase / 360.0) * 2.0 * PI;
/*
 *  Generate the fundamental component (2 x PI).
 */
    for(i=0; i< len; i++)
    {
        theta = ((double)i / sample_rate) * 2.0 * PI;
        indat[i] = (short) ((cos(theta) * 19660.0) + PAD);
    }
/*
 *  Generate second harmonic (4 x PI), 20 dB (10%) down.
 */
    for(i=0; i< len; i++)
    {
        theta  = ((double)i / sample_rate) * 4.0 * PI;
        theta += phase;
        indat[i] += (short) ((cos(theta) * 1966.0) + PAD);
    }
/*
 *  Generate third harmonic (6 x PI), 20 dB (10%) down.
 */
    for(i=0; i< len; i++)
    {
        theta  = ((double)i / sample_rate) * 6.0 * PI;
        theta += phase;
        indat[i] += (short) ((cos(theta) * 1966.0) + PAD);
    }
/*
 *  Generate fourth harmonic (8 x PI), 20 dB (10%) down.
 */
    for(i=0; i< len; i++)
    {
        theta  = ((double)i / sample_rate) * 8.0 * PI;
        theta += phase;
        indat[i] += (short) ((cos(theta) * 1966.0) + PAD);
    }
}
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//
//   Generate a signal with a second harmonic. The second harmonic
//   will be generated slightly out-of-phase. The phase-difference
//   is the input parameter, phase.  The maximum level is 32767.
//   We will have a second harmonic that is 20 dB down, i.e., 10%.
//   Therefore the fundamental will be 32767 - 3277 = 29490.
//
void HOW::acquire_dist_data(short *indat, double phase, int len)
{
    double theta;
    double sample_rate = (double) (len / NODES);
    int i;
    phase = (phase / 360.0) * 2.0 * PI;
//
//  Generate the fundamental component (2 x PI).
//
    for(i=0; i< len; i++)
    {
        theta = ((double)i / sample_rate) * 2.0 * PI;
        indat[i] = (short) ((cos(theta) * 29490.0) + PAD);
    }
//
//  Generate second harmonic (4 x PI), 20 dB (10%) down.
//
    for(i=0; i< len; i++)
    {
        theta  = ((double)i / sample_rate) * 4.0 * PI;
        theta += phase;
        indat[i] += (short) ((cos(theta) * 2949.0) + PAD);
    }
}
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//
//  Make a square-wave. This should have all odd-order harmonics.
//
void HOW::acquire_square_data(short *buf, double phase, int len)
{
    int i;
    int sample_rate = len / NODES;
    phase = phase;
    for(i = 0; i< len; i++)       /* Make a square wave */
        buf[i] = ((i / (sample_rate / 2)) & 1) ? 32767 : -32767;
}
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

int main(int argc, char *argv[])

{
	return how.showit();
}

int HOW::showit()
{
    char fbuf[0x100];
	char *cmds[0x04];
    int i, j, m, p2, len;
    short  *sp;
    double *win;
    double *indat;
    double *outdat;
    double *a, *b, *c;
    double f, fund, second, third, fourth, fifth;
    double db_const;
    FILE *file;
    len = ARRAY_SIZE;
	errno_t err;

	(void)_mkdir("Data");

//
//  First get some memory for our buffers.
//
    if((sp = new short [len]) == NULL)
    {
        fprintf(stderr, "Error from line %d, memory allocation\n", __LINE__);
        exit(EXIT_FAILURE);
    }
    if((indat = new double [len]) == NULL)
    {
        fprintf(stderr, "Error from line %d, memory allocation\n", __LINE__);
        exit(EXIT_FAILURE);
    }
    if((outdat = new double [len]) == NULL)
    {
        fprintf(stderr, "Error from line %d, memory allocation\n", __LINE__);
        exit(EXIT_FAILURE);
    }
    if((win    = new double [len]) == NULL)
    {
        fprintf(stderr, "Error from line %d, memory allocation\n", __LINE__);
        exit(EXIT_FAILURE);
    }
//
//  Next, before we acquire any data, we want to create a window with
//  which we will multiply the data.
//
//  The windows available are:
//     Hamming           Good for looking close to the carrier.
//     Hanning           Good for looking at spurs well away from carrier.
//     Kaiser            Good for looking at a broad spectrum.
//     Analogic          Good for looking close to carrier.
//     Weedon/Gauss      Good for looking at harmonics.
//     Raised cosine     Good for looking near strong spurs.
//     Blackman/Harris   Good for making Jeff shut up.
//
//  I will generate a Weedon/Gauss window in this demo.
//
//
    sig.MakeWeedonGauss(win, len);                  // Generate kernel
//
//  Since I will log all the data, as a natural log, I can predefine
//  a constant that will directly convert the result to 20 log10(A/B)
//  after the subtraction (instead of division) of A and B. This speeds
//  things up when converting to decibels, the normal spectrum display.
//
    db_const = 20.0 / log(10.0);                // Const for logN
//
//  Now I will acquire and process data in a loop as it would be
//  done on a continuous display. If I ever change the length, len,
//  I need to re-make a kernel with the new length. Otherwise I
//  use the same kernel over and over again.
//
    for(i=0; i< 4; i++)
    {
        acquire_dirty_data(sp, 90.0, len);       // Acquire the data
        sig.CvtShort2Double(indat, sp, len);     // Convert the data
        sig.MultArray(indat, win, len);          // Window the data
        f = sig.do_FFT(outdat, indat, SIG::WINDOW_TYPE_NONE, len);
//
//   Note that "do_FFT()" returns a fractional bin value that points
//   to the fundamental in the array. Why? You need to multiply that
//   value by a real number to point to a harmonic. Suppose you needed
//   the 10th harmonic of a signal with an actual carrier 1/2 way
//   between the two bins, 10 and 11. If you used integers, you would
//   look at 10 * 10 = 100 for the 10th harmonic  when, in fact you
//   need to look at 10.5 * 100 = 105 for the 10th harmonic.
//
//   The do_FFT() procedure, before it returns, calls a function
//   that looks at the data and returns its best estimate of where
//   the signal peak would have been. It does this by using the
//   amplitude of the two signals either side of the actual integer
//   bin.
//
//   In every case, perform the floating-point operations first, then
//   cast to an integer to index the array, i.e., i = (int)(f * 10.0);
//
//   Note that you only use 1/2 of the data returned from the FFT. The
//   second half is an alias.
//
//   There is a potential problem with very clean data. The result
//   of a FFT could be zero. We can't (or shouldn't) take the log
//   of zero so, we will sum some small amount of offset into the
//   data. We call this 'ZERO' in this code.
//
        sig.AddConst(outdat, ZERO, len/2);       // No minus infinity
        sig.LogN(outdat, len/2);                 // Log display data
        fund = outdat[(int)f];                   // Fundamental
        second = outdat[(int)(f * 2.0)];         // Second harmonic
        third  = outdat[(int)(f * 3.0)];         // Third harmonic
        fourth = outdat[(int)(f * 4.0)];         // Fourth harmonic
        fifth  = outdat[(int)(f * 5.0)];         // Fifth harmonic
        printf("Dirty data:\n");
        printf("Fund. frequency = %4.4f (fractional samples)\n", f);
        printf("Second harmonic = %4.4f dB\n", (second-fund) * db_const);
        printf("Third harmonic  = %4.4f dB\n", (third-fund)  * db_const);
        printf("Fourth harmonic = %4.4f dB\n", (fourth-fund) * db_const);
        printf("Fifth harmonic  = %4.4f dB\n", (fifth-fund)  * db_const);
        acquire_clean_data(sp, 90.0, len);      // Acquire the data
        sig.CvtShort2Double(indat, sp, len);    // Convert the data
        sig.MultArray(indat, win, len);         // Window the data
        f = sig.do_FFT(outdat, indat, SIG::WINDOW_TYPE_NONE, len);
//
//   Note that you only use 1/2 of the data returned from the FFT. The
//   second half is an alias.
//
        sig.AddConst(outdat, ZERO, len/2);       // No minus infinity
        sig.LogN(outdat, len/2);                 // Log display data
        fund = outdat[(int)f];                   // Fundamental
        second = outdat[(int)(f * 2.0)];         // Second harmonic
        third  = outdat[(int)(f * 3.0)];         // Third harmonic
        fourth = outdat[(int)(f * 4.0)];         // Fourth harmonic
        fifth  = outdat[(int)(f * 5.0)];         // Fifth harmonic
        printf("Clean data:\n");
        printf("Fund. frequency = %4.4f (fractional samples)\n", f);
        printf("Second harmonic = %4.4f dB\n", (second-fund) * db_const);
        printf("Third harmonic  = %4.4f dB\n", (third-fund)  * db_const);
        printf("Fourth harmonic = %4.4f dB\n", (fourth-fund) * db_const);
        printf("Fifth harmonic  = %4.4f dB\n", (fifth-fund)  * db_const);
//
//  If the screen doesn't have scaling software, we can do the following
//  to normalize everything for direct display on the screen.
//
        sig.SubConst(outdat, fund, len/2);       // Subtract Fo magnitude
        sig.MulConst(outdat, db_const, len/2);   // Scale to decibels
        printf("Clean data (pre windowed (1)) :\n");
        printf("Fundamental dB  = %4.4f dB\n", outdat[(int)f]);
        printf("Second harmonic = %4.4f dB\n", outdat[(int)(f * 2.0)]);
        printf("Third harmonic  = %4.4f dB\n", outdat[(int)(f * 3.0)]);
        printf("Fourth harmonic = %4.4f dB\n", outdat[(int)(f * 4.0)]);
        printf("Fifth harmonic  = %4.4f dB\n", outdat[(int)(f * 5.0)]);
    }
//
//   For folks who want to do all the dirty-work themselves, this executes
//   the FFT routine directly without the help of do_FFT(). The FFT does
//   and 'in place' FFT and requires that the result be extracted from
//   both the real and imaginary buffers as the sqrt(a^2 + b^2). Function
//   get_mag() does that for you. If speed is a problem, allocate the
//   vector buffers a, and b, only once in the program. However, you
//   __must__ zero those buffers, regardless.
//
    if((a = new double [len]) == NULL)
    {
        fprintf(stderr, "Error from line %d, memory allocation\n", __LINE__);
        exit(EXIT_FAILURE);
    }
    if((b = new double [len]) == NULL)
    {
        fprintf(stderr, "Error from line %d, memory allocation\n", __LINE__);
        exit(EXIT_FAILURE);
    }
    fprintf(stdout,"Writing %d files....", ArraySize(fd) * 3);
    fflush(stdout);

    for(i = 0; i < ArraySize(fd); i++)
    {
        for(m=0, p2=1; m< 32; m++, p2 <<= 1)     // m = log2(n)
            if(p2 >= len)
                break;
        memset(b, 0x00, p2 * sizeof(double));      // Zero imaginary array
        p2 >>=1;                                   // Divide by two
        (how.*fd[i].routine)(sp, fd[i].phase, len); // Get the data
        (sig.*fd[i].window)(win, len);             // Create the window
        sig.CvtShort2Double(a, sp, len);         // Convert the data
		sprintf_s(fbuf, sizeof(fbuf), "%s%s%s", Data, fd[i].fname, ".RAW");
        if((err = fopen_s(&file, fbuf, "w")) == 0)
        {
            for(j=0; j< len; j++)                // Write the whole array
                fprintf(file, "%d\t%f\n", j, a[j]);
            fclose(file);
        }
        sig.MultArray(a, win, len);              // Window the data
        sig.parse_fft(a, b, m, 0);               // Do actual FFT
        sig.get_mag(outdat, a, b, p2);           // Extract magnitude
        sig.AddConst(outdat, 0.1, p2);           // Prot against infinity
        sig.LogN(outdat, p2);                    // Log display data
        f = sig.IndexToFirstPeak(outdat, p2);    // Fundamental is here
        fund = outdat[(int)f];                   // Get fundamental mag
        sig.SubConst(outdat, fund, p2);          // Subtract Fo magnitude
		sig.MulConst(outdat, db_const, p2);
		sprintf_s(fbuf, sizeof(fbuf), "%s%s%s", Data, fd[i].fname, ".DAT");
        if((err = fopen_s(&file, fbuf, "w")) == 0)
        {
			fprintf(file, "#@graph.title  Spectrum\n");
            fprintf(file, "#@graph.xaxis.title Frequency\n");
            fprintf(file, "#@graph.yaxis.title Amplitude decibels\n");
            for(j=0; j< (int)(f * 10.0); j++)    // Write to 10th harmonic
                fprintf(file, "%d\t%f\n", j, outdat[j]);
            fclose(file);
        }
        sig.get_phase(outdat, a, b, p2);         // Get FFT phase
        sprintf_s(fbuf, sizeof(fbuf), "%s%s%s", Data, fd[i].fname, ".PHA");
        if((err = fopen_s(&file, fbuf, "w")) == 0)
        {
			fprintf(file, "#@graph.title  Phase\n");
            fprintf(file, "#@graph.xaxis.title Frequency\n");
            fprintf(file, "#@graph.yaxis.title Angle\n");
            for(j=0; j< (int)(f * 10.0); j++)    // Write to 10th harmonic
                fprintf(file, "%d\t%f\n", j, outdat[j]);
            fclose(file);
        }
    }
    fprintf(stdout, "Done!\n");
	show_frequency(outdat, len);

//
//  Clean up.
//

    delete[] a;
    delete[] b;
    delete[] sp;
    delete[] win;
    delete[] indat;
    delete[] outdat;

//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//
//  This shows how to perform the RSS filtering on noisy data
//

    a  = new double [len];                 // New data
	b  = new double [len];                 // Previous data
	c  = new double [len];                 // Filtered data
	
	cmds[0]= "Presto.exe";
	cmds[1]= "-xy";
	cmds[2]= fbuf;
	cmds[3]= NULL;
 
    get_noisy_fpdat(c, len);               // * Get noisy data
	sprintf_s(fbuf, sizeof(fbuf),"%s%s", Data, "Noisy.dat");
    fprintf(stdout, "Creating a noisy data file, %s\n", fbuf);
    (void) remove(fbuf);
	if((err = fopen_s(&file, fbuf, "w")) == 0)
    {
        fprintf(file, "#@graph.title Generated data\n");
        fprintf(file, "#@graph.xaxis.title Time\n");
        fprintf(file, "#@graph.yaxis.title Amplitude\n");
  
        for(j=0; j < len; j++)
                fprintf(file, "%d\t%f\n", j, c[j]);
        fclose(file);
		_spawnvp(_P_DETACH, cmds[0], cmds);
    }
    fprintf(stdout, "Filtering the data\n");
    for(i = 0; i < NR_FILTER; i++)
	{
 	    memcpy(b, c, len * sizeof(double));  // Copy previous
		get_noisy_fpdat(a, len);             // Get more data
		msl.FiltBySum(c, a, b, 2.0, len);    // Filter it
	}
	sprintf_s(fbuf, sizeof(fbuf), "%s%s", Data, "Filter.dat");
    fprintf(stdout, "Creating a filtered data file, %s\n", fbuf);
    (void) remove(fbuf);
	if((err = fopen_s(&file, fbuf, "w")) == 0)
    {
		fprintf(file, "#@graph.title Filtered data\n");
        fprintf(file, "#@graph.xaxis.title Time\n");
        fprintf(file, "#@graph.yaxis.title Amplitude\n");

        for(j=0; j < len; j++)
                fprintf(file, "%d\t%f\n", j, c[j]);
        fclose(file);
		_spawnvp(_P_DETACH, cmds[0], cmds);
    }
    fprintf(stdout, "Filtering the data another way\n");
    for(i = 0; i < NR_FILTER; i++)
	{
 	    memcpy(b, c, len * sizeof(double));  // Copy previous
		get_noisy_fpdat(a, len);             // Get more data
		sig.AddArray(c, a, b, len);          // Sum the data
	}
	sig.DivConst(c, (double)NR_FILTER, len);
    sprintf_s(fbuf, sizeof(fbuf), "%s%s", Data, "OtherFilter.dat");
    fprintf(stdout, "Creating a filtered data file, %s\n", fbuf);
    (void) remove(fbuf);
	if((err = fopen_s(&file, fbuf, "w")) == 0)
    {
		fprintf(file, "#@graph.title Filtered data\n");
        fprintf(file, "#@graph.xaxis.title Time\n");
        fprintf(file, "#@graph.yaxis.title Amplitude\n");

        for(j=0; j < len; j++)
                fprintf(file, "%d\t%f\n", j, c[j]);
        fclose(file);
		_spawnvp(_P_DETACH, cmds[0], cmds);
    }
    delete[] a;
    delete[] b;
    delete[] c;
    return 0;
}
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=    E N D    =-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

