Archive for March, 2015

Fun With Direct Digital Synthesis

You know that good feeling you get when a project comes together? Thats what I felt when I saw this sine wave on my oscilloscope.  It’s a Pulse Width Modulation (PWM) signal generated by an Arduino using a Direct Digital Synthesis (DDS) technique and a lookup table of amplitude values. It was the result of searching for a way to improve the tone quality of a morse code keyer project and I had gotten half wave, then quarter sine wave symmetry working in a demonstration sketch.

Quarter Wave Symmetry

Quarter Wave Symmetry
1000 Hz 256 Step Table

Digital to Analog Conversion is best done with optimized precision circuits. The relatively simple processor in an Arduino does not have a DAC device but it can do Pulse Width Modulation which can, in most cases, do the job adequately.

Pulse Width Modulation is widely documented. Basically it is a technique for generating an analog voltage from a digital signal. The device generates a high speed pulse train with varying pulse widths.  This is fed through a filter circuit which averages the energy in the pulses.  Fat pulses average to a higher voltage, skinny pulses average to a low voltage. By changing the pulse width over time, digital software can create a good analog representation of an arbitrary periodic signal. Other interesting applications are possible.

An Arduino generates PWM by setting up a timer. This is a counter with a circuit to trigger an interrupt at a specified count. The counter is incremented by the processor clock or, through a prescaler, a submultiple of the processor clock.  PWM is achieved by the Interrupt Processing Routine (ISR) manipulating the specified count at which the timer fires. DDS applications synthesize an arbitrary waveform by using a lookup table of values that the ISR accesses sequentially to change the terminal count.

Common non-DDS applications might be dimming an LED , feeding an RC servo, or controlling the speed of a motor, and the Arduino environment has a convenient analogWrite() function that handles all the details for you. Hybrid cars use PWM to control the traction motors.

I found many articles on generating DDS sound using the Arduino platform. I picked one from KO7M as a base to learn more about the method. Read Jeff’s weblog entry with source code here. His sketch sets up timer 2 for PWM with a 32 KHz base frequency.  Timer 2 is normally used by the Arduino for the tone() function and since my aim is to supersede the tone() function with something better, his example was a perfect starting point. In this article, I will only describe the changes I made in Jeff’s Interrupt Service Routine (ISR) and in the tables. Details of setting up the timers are covered in the source code and in many Internet sources. Here is a drawing of the Arduino connections:

Schematic for DDS Demo

Schematic for DDS Demo

Only five external parts needed including a button switch that stops and starts the sweep.

My test bed started out simple but it grew some (gruesome?).

Test Bed

Test Bed

Part 0: Analyzing and Parting Out the KO7M Sketch

After confirming that the sketch actually worked, I added code to sweep the frequency range from 50 to 2000 HZ so I could hear what it would sound like, and observe fidelity on my oscilloscope. Later I also added lines to generate a pulse at the start of the sine wave on a different pin. Using this pulse as an external sync source allowed the scope to display a steady waveform. I removed the sync code in the examples in this post for clarity.

This is the ISR from the original sketch (with small changes by me). TuningWord is the only parameter passed from the main part of the program. It is adjusted to set the particular PWM output frequency desired:

ISR(TIMER2_OVF_vect) {
  byte sine_table_index;
  static uint32_t phase_accumulator;
  
  // 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); 
}

Note Jeff is storing his sine table in flash memory. If you’re not familiar with that technique, read this tutorial. It is a good method for reducing the RAM footprint of an Arduino sketch. You can find my version of Jeff’s sketch here.  This is the 1000 Hz output from the original sketch with floating point math:

Original Floating Point

Original Floating Point
1000 Hz 256 Table

Part 1: Integer Conversion and Multiple Sine Tables

At this time, all floating point math in the example was converted to integer. To my surprise, it still works. There is more jitter on the waveforms due to integer truncation, but for keyer sidetone it is acceptable and the sketch sized dropped by about half.

I wanted to experiment with multiple lookup tables to get an idea of the minimum memory footprint possible while still producing acceptable sound. Libre Office Calc easily produced data with a 0-255 range for tables of length 256, 128, 64, 32 and 16 samples. See my spreadsheet here. Calc could also graph the tables to show how rough the waveform might be.  It was a simple matter to paste a table column into VI and join the values into C style  statements. These lines of data were inserted into the sketch and wrapped with #ifdef and #endif statements. At the top of the sketch is this code:

// Table sizes. Uncomment one and only one
//#define Table256
//#define Table128
#define Table64
//#define Table32
//#define Table16

Which allows me to choose a table size to evaluate at compile time. The wrapped tables look like:

#ifdef Table16
// 16 samples per period
#define Bits 12               // # of bits to shift 
PROGMEM prog_uchar sineTable[16] = {
  128,177,218,245,255,245,218,177,
  127,79,38,11,1,11,38,79
};
#endif // Table16

Bits is a constant definition that controls how far the ISR has to right shift the phase accumulator variable so the number that remains does not exceed the sine table size.

The interrupt handler changed little and looks like this:

ISR(TIMER2_OVF_vect) {
byte sineTableindex;
static unsigned int phaseAccumulator;

  // Update phase accumulator and extract the sine table index from it
  phaseAccumulator += tuningWord;
se
  // Right shift because we're only using the most significant bits of the INT
  sineTableindex = phaseAccumulator >> Bits;

  // Look up current amplitude value for the sine wave being constructed.
  OCR2A = pgm_read_byte_near(sineTable + sineTableindex);
}

Code for this full table sketch is here. These are waveforms from the Integer version:

Integer Version

Integer Version
Three Table Sizes

Part 2: Using Half Wave Symmetry of Sine Waves to Reduce Table Size

During my Google research I found articles pointing out that the size of the sine table could be reduced by a factor of two or four if you exploited the symmetry of a sine wave. Several pages described how to do this in FPGA hardware but I could not find a good example of symmetry use in C code and nothing at all for an Arduino. I tackled symmetry in two stages.

For half wave symmetry it is only necessary to duplicate the first half of the wave for the second half, inverting the second half as it is generated. First you need to know which half is being processed during the interrupt. The Most Signifigant Bit (MSB) of the phase accumulator variable has that information. It will be zero in the first half of the sine wave, one in the second half. I added a variable to record this bit and use a mask of 0x8000 to extract the bit. Then a new constant MSBMask was created inside the table definitions to turn off the high order MSB during table lookup. The tables now look like this example:

#ifdef Table16
// 16 samples per period
#define Bits 12
#define MSBMask 0x07
PROGMEM prog_uchar sineTable[8] = {
 128,177,218,245,255,245,218,177
};
#endif // Table16

Note the table sizes are half the original. This 16 stepe per period table is now only 8 bytes long.  The modified ISR with half wave symmetry now looks like:

ISR(TIMER2_OVF_vect) {
byte sineTableindex;                     // calculated index into sine table     
static unsigned int phaseAccumulator;    // running total of phase offset
unsigned int whichHalf;                  // for sine symmetry switching

  // Update phase accumulator and extract the sine table index from it
  phaseAccumulator += tuningWord;
  whichHalf = phaseAccumulator & 0x8000; // record which half of wave

  // Right shift because we're only using the most significant bits
  // of the INT. Leave only enough bits for one half table range
  sineTableindex = (phaseAccumulator >> Bits) & MSBMask;

  // Set current amplitude value for the sine wave being constructed.
  OCR2A = pgm_read_byte_near(sineTable + sineTableindex);

  // invert the second half of sine wave
  if(whichHalf) OCR2A = 255 - OCR2A;
}

The only additions to the ISR were the whichHalf variable, MSBMask which is defined with the sine table to remove the MSB, and the if statement at the end which does the actual inversion.

Code for the half wave symmetry version is here. These photos are from the half wave version:

Half Wave Symmetry Version

Half Wave Symmetry Version
Three Table Sizes

Part 3: Using Quarter Wave Sine Symmetry to Reduce Table Size

Exploiting quarter wave symmetry is more complicated. It is a left/right symmetry while half wave symmetry is up/down. First you have to know which quarter of the wave is being generated. The second most significant bit in the phase accumulator variable has this information. It will be 0 if in the first or third quarter, 1 if second or fourth. A new state variable and a mask of 0x4000 records this bit.  Since we are using less of the lookup tables, the MSBMask is smaller, it now needs to discard two MSB bits. The tables now look like:

#ifdef Table16
// 16 samples per period
#define Bits 12
#define MSBMask 0x03
PROGMEM prog_uchar sineTable[5] = {
  128, 177, 218, 245, 255
};
#endif // Table16

Note that the sixteen sample table has been reduced to five (not four) bytes. The 256 level table is only 65 bytes long! The odd byte at the end was necessary to work around a glitch – the table read from flash would look up the same value twice at the beginning and end of the even quadrant. Adding one to the index and one extra table byte solves that issue. The ISR now looks like this:

ISR(TIMER2_OVF_vect) {
byte sineTableindex;                     // calculated index into sine table     
static unsigned int phaseAccumulator;    // running total of phase offset
unsigned int whichHalf;                  // for sine symmetry
unsigned int whichQtr;                   // for sine symmetry
static boolean pulsed;                   // for unique scope sync

  // Update phase accumulator and extract the sine table index from it
  phaseAccumulator += tuningWord;
  whichHalf = phaseAccumulator & 0x8000; // record which half of wave
  whichQtr =  phaseAccumulator & 0x4000; // record which quarter of wave

  // Right shift because we're only using the most significant bits
  // of the INT. Leave only enough bits for one quarter table range
  sineTableindex = (phaseAccumulator >> Bits) & MSBMask;

  // Look up current amplitude value for the sine wave being constructed.
  if(whichQtr) {                         // in second or fourth quarter
  // +1 works around a glitch where same data at the quarter 
  // transitions was looked up twice in a row
    OCR2A = pgm_read_byte_near(sineTable + (MSBMask - sineTableindex + 1));
  } 
  else {                                 // in first or third quarter
    OCR2A = pgm_read_byte_near(sineTable + sineTableindex);
  }
  // invert the second half of sine wave
  if(whichHalf) OCR2A = 255 - OCR2A;
}

This code does quarter wave symmetry first, then wraps half wave symmetry around that for the final output.  The changes to add quarter wave are the whichQtr variable and it’s masking line, and the if(whichQtr) statement which inverts the table lookup left to right. MSBMask in in the byte read statement because it is the number of steps that exactly form a quarter of the wave form.

Code for the quarter wave symmetry sketch is here. Here are scope traces of the quarter wave symmetry version:

Quarter Wave Symmetry Version

Quarter Wave Symmetry Version
Three Table Sizes

Conclusion:

These changes to the original Interrupt Service Routine allowed the table sizes to be reduced by a factor of four with hopefully, not too much additional overhead. Any ISR overhead added is executed 32000 every second so you have to be conservative.  CPU time in the ISR is traded for a smaller program.

All three of these sketches can be useful for arbitrary waveform synthesis by constructing an appropriate table.  The base sketch would work for any periodic waveform, including those that have no symmetry either up/down or left/right. A sawtooth wave for example.  Half wave symmetry would be good for representing an underdamped or overdamped square wave.  Quarter wave symmetry works for sine waves, as shown here, or a triangle wave.

The oscilloscope traces are interesting by themselves. These sketches and a scope could be used to demonstrate microcontroller techniques at your local science fair. My scope is an 80’s vintage Tektronix clone. On the screen you can see the generated sine waves and by zooming in, see the artifacts of PWM happening at a 32 KHz rate. Use the coarse 16 step table and the staircase model of DDS is clear. You can view the raw PWM by simply disconnecting the 0.1 Ufd integrating capacitor from the scope input. Interesting though, disconnecting the filter capacitor does not affect the audio much because your ears can’t respond to 32 Khz sound.

DDS quantization is most visible at low frequencies, using the smaller tables. Here you can clearly see the 16 steps and 32Khz fuzz on the trace.

50 Hz 16 Step Waveforms

50 Hz 16 Step Waveforms

Footnote:

The 64 step sine synthesizer with quarter wave symmetry was successfully added to my keyer project.  I see no degradation of the morse timing at 50 WPM.  A brief video is here.  Also I recorded MP3s of the tone before and after.  The tone produced by the stock Ardino tone() function is here.  The tone from the 64 step (16 values) DDS function is here.

Advertisements