DDS.c

#include "Header/main.h"

//The look-up-table for the sine values / phase accumulator

static double lut[LUT_SIZE];

//Phaseaccumulator

static double phase;

void DDS_init()

{

  //Fill LUT

  Uint8 aa = 0;

  for(aa = 0; aa < LUT_SIZE; aa++)

  {

    //Calculate sine samples, Int16  -32.768 to 32.767 => amplify by 32767 (2^15-1)

    lut[aa] = (double) sin(PI * aa / LUT_SIZE);

  }

  //Preset variable for the phaseaccumulator

  phase = 0;

}

/// <summary>

/// Direct Digital Synthesis (DDS) using Lanczos-Splines for interpolation.

/// It is desinged for a frequency range from 100Hz-10kHz

/// This method accumulates phase jumps non-quantized and returns corresponding

/// sine values from a look-up-table, in which only half a wave is stored.

/// Values between two LUT-indexes are interpolated by a Lanczos-Spline. The spline starts one index before

/// the arbitrary value, and ends two indexes after.

/// </summary>

/// <param name="A">Amplitude</param>

/// <param name="f">Frequency</param>

/// <returns>Array with calculated sine values</returns>

void DDS(Int16 *target, double A, Int16 f)

{

  //Binary tuning word (Normally multiplied by 2 but since we are using 2 Buffers ping/pong

  //we have to output only 25 times the ping, and 25 times the pong buffer

  double jumpSize = (double)f * LUT_SIZE / F_CLK;

  //Current position in the LUT

  Int16 index = 0;

  //Current iteration step of the accumulation

  Int32 iterator = 0;

  //Count variable to loop through the output buffer

  Uint16 aa = 0;

  for(aa = 0; aa < OUT_SIZE; aa++)

  {

    // Lanczos a=3 (x0' = floor(x0))

    // Î(x0) = L(x0-x0'-1)*lut(x0'-1) + L(x0-x0')*lut(x0') + L(x0-x0'+1)*lut(x0'+1) + L(x0-x0'+2)*lut(x0'+2)

    //Result

    double result = 0;

    Uint8 bb = 0;

    //Loop through the sum (Sigma)

    for(bb = 0; bb < 4; bb++)

    {

      //Set iterator

      iterator = (Int32)phase + bb - 1;

      //Set index for the LUT

      index = iterator % (2 * LUT_SIZE);

      //Get the sine value from the LUT, corresponding to the index.

      //Actually the LUT is a ringbuffer, if the index is negativ =>

      //Get the value from the end. Since we have only stored half

      //a sine wave, we also have to invert the value

      //Also calculate a step of the accumulation: L(x0-x0'-bb-1) * lut(x0'-bb-1)

      if(index < 0)

      {

        //If index is below the size of the lut, set to end and invert

        result += Lanczos(phase - iterator) * -lut[LUT_SIZE + index];

      }

      else if(index < LUT_SIZE)

      {

        //First half

        result += Lanczos(phase - iterator) * lut[index];

      }

      else

      {

        //If index is beyond the size of the lut => second half, set to start but invert value

        result += Lanczos(phase - iterator) * -lut[index - LUT_SIZE];

      }

    }

    //Amplify the result by amplitude A

    *(target + aa) = A * result * 32767;

    //Increment phase

    phase += jumpSize;

  }

  //Decrement phase;

  phase = (Int32)phase % OUT_SIZE;

}