ADS131M08: How are the metrology-fir-table.c generated in ADCEnergyMetrologyLibrary SW_1.03?

/*******************************************************************************
 *  generate-phase-shift-table.c -
 *
 *  Copyright (C) 2011 Texas Instruments Incorporated - http://www.ti.com/ 
 * 
 *  Redistribution and use in source and binary forms, with or without 
 *  modification, are permitted provided that the following conditions 
 *  are met:
 *
 *    Redistributions of source code must retain the above copyright 
 *    notice, this list of conditions and the following disclaimer.
 *
 *    Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the 
 *    documentation and/or other materials provided with the   
 *    distribution.
 *
 *    Neither the name of Texas Instruments Incorporated nor the names of
 *    its contributors may be used to endorse or promote products derived
 *    from this software without specific prior written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 
 *  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 
 *  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 *  A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 
 *  OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 
 *  SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
 *  LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 *  DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 *  THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
 *  (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
 *  OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 ******************************************************************************/

/* Generate a table of single tap FIR coeffs for phase shifting signals by
   +-1/2 a sample. */

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

//gcc -o generate-phase-shift-table.exe generate-phase-shift-table.c
// ./generate-phase-shift-table


int betas[5000];
int As[5000];


#define SAMPLE_RATE           2,929.69
#define FREQUENCY             50
#define STEPS                 256      
#define X_PARAMETER           1

//#define  M_PI                 3.14159265      


void generate_fir_table(float sampling_rate, float frequency, int steps)
{
    double beta;
    double omega;
    double phi;
    double D;
    double A;
    double a;
    int i;

    omega = 2.0*M_PI*frequency/sampling_rate;
    phi = omega/2.0;
    for (i = 0;  i < steps;  i++, phi -= omega/(double) steps)
    {
        /* phi is now in radians */
        D = phi/omega;
        a = (1.0 - 2.0*D)*(1.0 - 2.0*D)*cos(omega)*cos(omega) + 4.0*D*(1.0 - D);
        beta = -((1.0 - 2.0*D)*cos(omega) - sqrt(a))/(2.0*(1.0 - D));
        A = 1.0/sqrt((cos(omega) + beta)*(cos(omega) + beta) + sin(omega)*sin(omega));
        /* Scale things */
        beta *= 32768.0;
        /* Never allow beta to quite reach 32768 */
        if (beta >= 32768.0)
            beta -= 1.0;
        A *= 16384.0;
        betas[i] = beta;
        As[i] = A;
    }
}

void fine_tune_filter_table_gains(float sampling_rate, float frequency, int steps)
{
    double omega;
    double rms_in;
    double rms_out;
    double last_amp;
    double amp;
    int i;
    int j;

    /* Optimise the filter taps, through simulation of the filter */
    omega = 2.0*M_PI*frequency/sampling_rate;
    //printf("omega = %f %f\n", omega, omega*180.0/M_PI);
    for (i = 0;  i < steps;  i++)
    {
        do
        {
            rms_in = 0;
            rms_out = 0;
            last_amp = 0.0;
            for (j = 1;  j <= 10000;  j++)
            {
                amp = 1500*sin(j*omega);
                rms_in += amp*amp;
                amp += (last_amp*betas[i])/32768.0;
                amp *= As[i]/16384.0;
                rms_out += amp*amp;
                last_amp = amp;
            }
            rms_in = sqrt(rms_in/10000.0);
            rms_out = sqrt(rms_out/10000.0);
            //printf ("%d %d %f %f\n", betas[i], As[i], rms_in, rms_out);
            if (rms_out < rms_in)
                As[i]++;
        }
        while (rms_out < rms_in);
    }
#if 0
    /* See how well the tune went */
    for (i = 0;  i < steps;  i++)
    {
        rms_in = 0;
        rms_out = 0;
        last_amp = 0.0;
        for (j = 1;  j <= 10000;  j++)
        {
            amp = 1500.0*sin(j*omega);
            rms_in += amp*amp;
            amp += (last_amp*betas[i])/32768.0;
            amp *= As[i]/16384.0;
            //rms_out += amp*amp;
            rms_out += amp*1500.0*sin(j*omega + (i - steps/2)*omega/(double) steps);
            last_amp = amp;
        }
        //rms_in = sqrt(rms_in/10000.0);
        //rms_out = sqrt(rms_out/10000.0);
        printf ("%d %d %f %f %f\n", betas[i], As[i], rms_in, rms_out, rms_out/rms_in);
    }
#endif
}

void output_fir_table(float sampling_rate, float frequency, int steps, const char *name)
{
    double omega;
    double phi;
    int i;

    omega = 2.0*M_PI*frequency/sampling_rate;
    printf("const int16_t %s[%d][2] =\n", name, steps);
    printf("{\n");
    phi = omega/2.0;
    for (i = 0;  i < steps - 1;  i++, phi -= omega/(double) steps)
        printf("    { %6d, %6d }, /* %3d: %7.4f deg */\n", betas[i], As[i], i, phi*180.0/M_PI);
    printf("    { %6d, %6d }  /* %3d: %7.4f deg */\n", betas[i], As[i], i, phi*180.0/M_PI);
    printf("};\n");
}

int main( void )
{
    float sampling_rate;
    float frequency;
    int steps;
    int fine_tune_gain;

   
    sampling_rate = SAMPLE_RATE;
    frequency = FREQUENCY;
    steps = STEPS;
    
    fine_tune_gain = (X_PARAMETER == 'x');

    generate_fir_table(sampling_rate, frequency, steps);
    if (fine_tune_gain)
        fine_tune_filter_table_gains(sampling_rate, frequency, steps);

    printf("/* This file was generated automatically by generate-phase-shift-table */\n");
    printf("#include <inttypes.h>\n");
    printf("\n");
    printf("#include <emeter-toolkit.h>\n");
    printf("\n");
    printf("#include \"emeter-metrology.h\"\n");
    printf("#include \"metrology-structs.h\"\n");
    printf("\n");
    printf("#if !defined(__HAS_SD_ADC__)  ||  defined(REACTIVE_POWER_BY_QUADRATURE_SUPPORT)\n");
    printf("\n");
    printf("/* This FIR table is used to create accurate 90 degree phase shifts */\n");
    printf("\n");
    printf("/* FIR coefficients to correct specific phase shift angles at %.2fHz\n", frequency);
    printf("   and sampling %.2f times per second. There are %d steps\n", sampling_rate, steps);
    printf("   in this table. */\n");
    printf("\n");
    printf("#if SAMPLES_PER_10_SECONDS != %d\n", (int) (sampling_rate*10.0f));
    printf("#error This FIR table is for %.2f samples per second.\n", sampling_rate);
    printf("#endif\n");
    printf("#if MAINS_NOMINAL_FREQUENCY != %d\n", (int) frequency);
    printf("#error This FIR table is for %.2fHz.\n", frequency);
    printf("#endif\n");
    printf("#if PHASE_SHIFT_FIR_TABLE_ELEMENTS != %d\n", steps);
    printf("#error This FIR table is for %d steps per sample.\n", steps);
    printf("#endif\n");
    printf("\n");

    output_fir_table(sampling_rate, frequency, steps, "phase_shift_fir_coeffs");

    printf("\n");
    printf("#endif\n");
    return  0;
}