Direct Digital Synthesis (DDS)

UPDATED

With all the recent work on software for the Si570 device, this got me thinking a bit about the process of direct digital synthesis.  I went back and read again a great tutorial on the topic by Analog Devices MT-085 Fundamentals of Direct Digital Synthesis and decided to experiment with a software DDS using the Arduino as a test platform.

Basically the architecture of a DDS system has a stable clock driving a look up of sinusoidal information where one or more cycles of a sine wave (or for that matter, any other arbitrary waveform) are fed in sequence to a digital to analog converter to produce the final output.  Each of these look ups can be stored in a read-only memory and the value fetched represents the corresponding digital amplitude of the signal at each clock tick.

For my experimentation, I will construct a sinusoidal waveform using 8 bit data representing the amplitude of the signal at each clock tick.  To construct the table of sinusoidal amplitude information, I will retrieve a value from 0-255 which will map to a signal where the sinusoidal waveform crosses zero at the middle of this range (e.g. 128) and 256 samples for one period of the waveform.

Each address in the look up table corresponds to a phase point on the sine wave from 0 to 360 degrees.  In reality only data for 90 degrees (one quadrant of a circle) would be required as the quadrant is indicated by the two most significant bits.  The look up table contains the corresponding digital amplitude information for one complete cycle of a sine wave used to drive a DAC (digital to analogue converter).

If the look up table contained an entry for each bit of a 32 bit index into the table, there would be 2^32 output values before the index overflowed.  A sine wave thus constructed would have a frequency equal to the input clock frequency divided by 2^32 as it would take that many data points to reconstruct the sine wave.  Let's call this the M=1 case (one increment per clock tick).

If we increment through the table a little faster, say M=2 then the table index will roll over twice as fast and the output frequency would be doubled.  For an n-bit index (in DDS systems this is called a phase accumulator) there are 2^n possible phase points (data points in our look up table).  Let's call M the amount that the phase accumulator (index) is incremented on each clock cycle.  If fc is the clock frequency, then the frequency of the output sine wave is determined by the following formulae, also known as the DDS tuning equation.

    fout = (M * fc) / (2^n)

The frequency resolution of a DDS system is equal to fc/2^n.  If n = 32, this represents a resolution greater than one part in four billion. (2^32 = 4,294,967,296).  This is rather impractical and in general unnecessary.  In most DDS systems not all of the bits of the phase accumulator are used in the look up table and only the upper (most significant) bits are used.  This can significantly reduce the size of the look up table without affecting the frequency resolution.  This truncation of the phase does however add a small amount of phase noise to the final output.

There is a lot more information in the reference I mention above (MT-085) and well worth a read by anyone interested in DDS technology.

The first task at hand is to construct a table of 256 sine values representing the amplitude of the sinusoid for one sine period.  I wrote a little C# application to generate the necessary data structure text:

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace GenSineData
{
  class Program
  {
    static void Main(string[] args)
    {
      Console.WriteLine("PROGMEM prog_uchar sine256[] = ");
      Console.WriteLine("{");
      Console.Write("  ");

      for (int i = 0; i < 256; i++)
      {
        double angle = 360.0 / 256.0 * Convert.ToDouble(i);
        double radians = angle * Math.PI / 180;
        double sin = Math.Sin(radians);
        
        // Map range from -1 to +1 to 0 to 255
        int amplitude = Convert.ToInt32(map(sin, -1, 1, 0, 255));
        Console.Write(amplitude.ToString() + (i != 255 ? "," : ""));
        if (i > 0 && (i+1) % 16 == 0)
        {
          Console.WriteLine();
          Console.Write("  ");
        }
      }
      Console.WriteLine("};");
    }

    // Map from one range to another
    static double map(double x, double in_min,
                      double in_max, double out_min, double out_max)
    {
      return (x - in_min) * (out_max - out_min) / (in_max - in_min) + out_min;
    }
  }
}

The output of this application can be pasted into my Arduino code producing a program memory-based lookup table.

PROGMEM prog_uchar sine256[] =
{
  128,131,134,137,140,143,146,149,152,155,158,162,165,167,170,173,
  176,179,182,185,188,190,193,196,198,201,203,206,208,211,213,215,
  218,220,222,224,226,228,230,232,234,235,237,238,240,241,243,244,
  245,246,248,249,250,250,251,252,253,253,254,254,254,255,255,255,
  255,255,255,255,254,254,254,253,253,252,251,250,250,249,248,246,
  245,244,243,241,240,238,237,235,234,232,230,228,226,224,222,220,
  218,215,213,211,208,206,203,201,198,196,193,190,188,185,182,179,
  176,173,170,167,165,162,158,155,152,149,146,143,140,137,134,131,
  128,124,121,118,115,112,109,106,103,100,97,93,90,88,85,82,
  79,76,73,70,67,65,62,59,57,54,52,49,47,44,42,40,
  37,35,33,31,29,27,25,23,21,20,18,17,15,14,12,11,
  10,9,7,6,5,5,4,3,2,2,1,1,1,0,0,0,
  0,0,0,0,1,1,1,2,2,3,4,5,5,6,7,9,
  10,11,12,14,15,17,18,20,21,23,25,27,29,31,33,35,
  37,40,42,44,47,49,52,54,57,59,62,65,67,70,73,76,
  79,82,85,88,90,93,97,100,103,106,109,112,115,118,121,124
  };

Importing this data into Excel and graphing it, we find a nice sinusoidal waveform with zero crossings at the value 128 and all amplitude values between 0 and 255 using only 256 data points and 8 bit amplitude values.

So far, so good...

To compute the necessary tuning word for the DDS for the desired output frequency:

  M = (2^32) * desired_frequency / clock_frequency

For my purposes, I intend to use a clock_frequency of approximately 32kHz by dividing down the Arduino master clock of 16MHz.  I will use a timer interrupt at 32kHz and a 32 bit phase accumulator.  The upper 8 bits will be used as the index into the above sine wave table used to drive the Arduino PWM DAC.  Of course a low pass filter will be required to remove high frequency components of the audio output.

Ok, so now we put together a little Arduino test application to pull all these concepts together.  I am using an ATMega2560 for my testing but any Arduino would work so long as the pin differences are taken into account between the platforms.  More on this later.

Here is my DDSTest sketch:

// Sine wave generator using DDS techniques
// Jeff Whitlatch - ko7m

#include "avr/pgmspace.h"

// Single period sine wave table.
// Amplitudes are 0-255 with 128 as zero crossing.  256 samples per period.
PROGMEM prog_uchar sine256[] =
{
  128,131,134,137,140,143,146,149,152,155,158,162,165,167,170,173,
  176,179,182,185,188,190,193,196,198,201,203,206,208,211,213,215,
  218,220,222,224,226,228,230,232,234,235,237,238,240,241,243,244,
  245,246,248,249,250,250,251,252,253,253,254,254,254,255,255,255,
  255,255,255,255,254,254,254,253,253,252,251,250,250,249,248,246,
  245,244,243,241,240,238,237,235,234,232,230,228,226,224,222,220,
  218,215,213,211,208,206,203,201,198,196,193,190,188,185,182,179,
  176,173,170,167,165,162,158,155,152,149,146,143,140,137,134,131,
  128,124,121,118,115,112,109,106,103,100,97,93,90,88,85,82,
  79,76,73,70,67,65,62,59,57,54,52,49,47,44,42,40,
  37,35,33,31,29,27,25,23,21,20,18,17,15,14,12,11,
  10,9,7,6,5,5,4,3,2,2,1,1,1,0,0,0,
  0,0,0,0,1,1,1,2,2,3,4,5,5,6,7,9,
  10,11,12,14,15,17,18,20,21,23,25,27,29,31,33,35,
  37,40,42,44,47,49,52,54,57,59,62,65,67,70,73,76,
  79,82,85,88,90,93,97,100,103,106,109,112,115,118,121,124
  };

// Useful macros for setting and resetting bits
#define cbi(sfr, bit) (_SFR_BYTE(sfr) &= ~_BV(bit))
#define sbi(sfr, bit) (_SFR_BYTE(sfr) |= _BV(bit))

// DDS frequency and reference clock
double dds_frequency = 0.0;
const double ref_frequency = (16000000/500);

// These must all be marked as volatile as they are used in an interrupt service routine
volatile byte sine_table_index;
volatile uint32_t phase_accumulator;
volatile uint32_t tuning_word;

void setup()
{
  // PWM output for timer2 is pin 10 on the ATMega2560
  // If you use an ATMega328 (such as the UNO) you need to make this pin 11
  // See spreadsheet here
  pinMode(10, OUTPUT);      // Timer 2 PWM output on mega256 is pin 10
  
  // Set up timer2 to a phase correct 32kHz clock
  timer2Setup();

  // disable interrupts to avoid timing distortion
  cbi (TIMSK0,TOIE0);    // Disable timer 0.  Breaks the delay() function
  sbi (TIMSK2,TOIE2);    // Enable timer 2.

  // Set up initial DDS frequency and calculate the timing word
  dds_frequency = 1000;
  tuning_word = pow(2,32) * dds_frequency / ref_frequency;
}

// Nothing to do here.  Everything is interrupt driven
void loop()
{
}

// Setup timer2 with prescaler = 1, PWM mode to phase correct PWM
// See the ATMega2560 datasheet for all the gory details
void timer2Setup()
{
  // Clock prescaler = 1
  sbi (TCCR2B, CS20);      // 001 = no prescaling
  cbi (TCCR2B, CS21);
  cbi (TCCR2B, CS22);

  // Phase Correct PWM
  cbi (TCCR2A, COM2A0);    // 10 = clear OC2A compare match
  sbi (TCCR2A, COM2A1);

  // Mode 1
  sbi (TCCR2A, WGM20);     // See table 20-8 in datasheet
  cbi (TCCR2A, WGM21);
  cbi (TCCR2B, WGM22);
}

// Timer 2 interrupt service routine (ISR) is used to generate
// the timebase reference clock for the DDS generator at 32kHz.
ISR(TIMER2_OVF_vect)
{
  // Update phase accumulator and extract the sine table index from it
  phase_accumulator += tuning_word;
  sine_table_index = phase_accumulator >> 24;  // Use upper 8 bits as index
  
  // Set current amplitude value for the sine wave being constructed.
  OCR2A = pgm_read_byte_near(sine256 + sine_table_index); 

}

Connecting a powered speaker to pin 10 allows me to hear the 1kHz tone.  The PWM signal on pin 10 looks like this, though when this picture was taken, a 1.5kHz tone was being generated.

So, far so good...  A low pass filter will be required or at least some level of integration to turn this into a sine wave.  Since we are sampling at 32kHz, Nyquist says that the highest frequency we can produce is 16kHz.  So by setting a simple RC low pass filter cutoff frequency near this frequency should be adequate.  I will choose 15kHz for simplicity.

  cutoff frequency = 1/ (2 * pi * R * C)

Using a 1000 ohm resistor for R:

  15000 = 1/(2000 * pi * C)
  C = 1/(30000000 * pi)
  C = 1.06103295e-8 farad or 0.0106103295 microfarads.  (.01 uF will do thanks...)

The following is the result using this simple integrator:

Fun!  Ok, so next installment, I hope to do something useful with all this...

More to come...