/*
 *   File: Amod.c      Created 17-SEP-2007  Richard B. Johnson
 *
 *   RichardBJohnson@comcast.net, rjohnson@analogic.com
 *
 *   Program to generate Amplitude Modulation waveform data.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <errno.h>
#include <math.h>

#define SAMPLES     10e6                // Number of samples
#define SECONDS     0.01                // Seconds of data
#define MOD_FREQ    1000.0              // Modulation frequency
#define CAR_FREQ    0.5e6               // Carrier frequency
#define CAR_LEVEL   1.0                 // Carrier level
#if !defined(PI)
# define PI 3.1415926535897932384626433832795029
#endif
#if !defined(TRUE)
# define TRUE 1
#endif
#if !defined(FALSE)
# define FALSE 0
#endif

#define FAIL -1
#define ERRORS(s) \
{                 \
    fprintf(stderr, "Error on line %d, file %s, (%s)\n", \
    __LINE__, __FILE__, strerror(errno));                \
    exit(EXIT_FAILURE);  \
}

#define No_CARRIER_SHIFT

//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//
//   Generate a sine-wave with respect to the input parameters
//
static void make_wave(double *samp, double freq, double sec, double mag)
{
    double twopi, rate,  turns, theta;
    size_t i, lim;
    twopi      = PI * 2.0;
    rate       = SAMPLES;
    lim        = (size_t)(sec * SAMPLES);
    for(i= 0;  i < lim ; i++)
    {
        turns = ((double)i / rate) * twopi;
        theta = turns * freq;
        *samp++ = sin(theta) * mag;
    }
}
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//
//   Sum every element in an array with a constant
//
static void sum_array_const(double *array, double val, size_t len)
{
    size_t i;
    for(i=0; i< len; i++)
        *array++ += val;
}
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//
//   Multiply the values in two source arrays and put the result in
//   the destination array
//
static void mult_array(double *dest, double *src1, double *src2, size_t len)
{
    size_t i;
    for(i=0; i< len; i++)
        *dest++ = *src1++ * *src2++;
}
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//
//  Make every element in the array positive. Don't call fabs()!
//
static void abs_array(double *array, size_t len)
{
    size_t i;
    for(i=0; i< len; i++)
    {
        if(*array < 0.0)
            *array *= -1.0;
        array++;
    }
}
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//
//   Filter the array elements using a simple IIR filter.
//
static void flt_array(double *array, size_t len)
{
    double var, *ptr;
    size_t i, j;
    for(j = 0; j < 16; j++)
    {
        ptr = array;
        var = *ptr;
        for(i=0; i< len; i++)
        {
            *ptr += var;
            *ptr /= 2.0;
            var = *ptr++;
        }
    }
}
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//
//  Clip any negative levels from the array.
//
static void clip_array(double *array, size_t len)
{
    size_t i;
    for(i=0; i< len; i++)
    {
        if(*array < 0.0)
            *array = 0.0;
        array++;
    }
}
#ifdef CARRIER_SHIFT
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//
//  Return the average value of all the elements in an array.
//
static double ave_array(double *array, size_t len)
{
    size_t i;
    double ave=0.0;
    for(i=0; i< len; i++)
        ave += *array++;
    ave /= (double) len;
    return ave;
}
#endif
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//
//  Define some work to be done.
//
static struct {
    const char *fname;
    double carrier;
    double pcent;
    double freq;
    int    clip;
    int    det;
    } work[] = {
    { "100pc.txt"        , CAR_LEVEL, 1.00, MOD_FREQ, FALSE, FALSE },
    { "100pc-det.txt"    , CAR_LEVEL, 1.00, MOD_FREQ, FALSE, TRUE  },
    { "100dsb.txt"       ,         0, 1.00, MOD_FREQ, FALSE, FALSE },
    { "125pc.txt"        , CAR_LEVEL, 1.25, MOD_FREQ, FALSE, FALSE },
    { "125pc-det.txt"    , CAR_LEVEL, 1.25, MOD_FREQ, FALSE, TRUE  },
    { "125pclip.txt"     , CAR_LEVEL, 1.25, MOD_FREQ, TRUE,  FALSE },
    { "125pclip-det.txt" , CAR_LEVEL, 1.25, MOD_FREQ, TRUE,  TRUE  },
    { "150pc.txt"        , CAR_LEVEL, 1.50, MOD_FREQ, FALSE, FALSE },
    { "150pc-det.txt"    , CAR_LEVEL, 1.50, MOD_FREQ, FALSE, TRUE  },
    { "150pclip.txt"     , CAR_LEVEL, 1.50, MOD_FREQ, TRUE,  FALSE },
    { "150pclip-det.txt" , CAR_LEVEL, 1.50, MOD_FREQ, TRUE,  TRUE  },
    { "200pc-det.txt"    , CAR_LEVEL, 2.00, MOD_FREQ, FALSE, TRUE  },
    { "200pclip.txt"     , CAR_LEVEL, 2.00, MOD_FREQ, TRUE,  FALSE },
    { "200pclip-det.txt" , CAR_LEVEL, 2.00, MOD_FREQ, TRUE,  TRUE  },
    };
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
 
#define ArraySize(a) (sizeof(a) / sizeof(a[0]))

int main(int args, char *argv[])
{
    double *carrier, *modulation, *signal;
    FILE *file;
    size_t i, j, len, lim;
    len = (size_t)(SAMPLES * SECONDS) * sizeof(double);
    if((carrier = malloc(len)) == NULL)
        ERRORS("malloc");
    if((modulation = malloc(len)) == NULL)
        ERRORS("malloc");
    if((signal = malloc(len)) == NULL)
        ERRORS("malloc");
    lim = (size_t)(SAMPLES * SECONDS);
    make_wave(carrier, CAR_FREQ, SECONDS, CAR_LEVEL);

    for(j = 0; j < ArraySize(work); j++)
    { 
        make_wave(modulation, work[j].freq, SECONDS, work[j].pcent);
        sum_array_const(modulation, work[j].carrier, lim); 
        if(work[j].clip == TRUE)
            clip_array(modulation, lim);
#ifdef CARRIER_SHIFT
        printf("Carrier shift = %2.2f%%\n",
            (ave_array(modulation, lim) - work[j].carrier) * 100.0);
#endif
        mult_array(signal, modulation, carrier, lim);
        if(work[j].det == TRUE)
        {
            abs_array(signal, lim);     // Detect the signal
            flt_array(signal, lim);     // Filter RF out
        }
        (void)remove(work[j].fname);
        printf("Creating file %s\n", work[j].fname);
        if((file = fopen(work[j].fname, "w")) == NULL)
            ERRORS("fopen");
        for(i = 0;  i < lim ; i++)
            fprintf(file, "%u\t%f\n", i, signal[i]);
        (void)fclose(file);
    }
    free(carrier);
    free(modulation);
    free(signal);
    return 0;
}