Arduino source code: 64-bit fixed-point DDS calculations

DDSCalcTest.ino

// Fixed point exercise for 60 kHz crystal tester

#include <avr/pgmspace.h>

char Buffer[10+1+10+1];               // string buffer for long long conversions

#define GIGA 1000000000LL
#define MEGA 1000000LL
#define KILO 1000LL

struct ll_fx {
  uint32_t low;
  uint32_t high;
};

union ll_u {
  uint64_t fx_64;
  struct ll_fx fx_32;
};

union ll_u CtPerHz;                   // will be 2^32 / 125 MHz
union ll_u HzPerCt;                   // will be 125 MHz / 2^32

union ll_u One;                       // 1.0 as fixed point
union ll_u Tenth;                     // 0.1 as fixed point

union ll_u TenthHzCt;                 // 0.1 Hz in counts

union ll_u Oscillator;                // nominal oscillator frequency
union ll_u OscOffset;                 // oscillator calibration offset

union ll_u TestFreq,TestCount;        // useful variables
union ll_u TempFX;


//-----------
// Round scaled fixed point to specific number of decimal places: 0 through 8
// You should display the value with only Decimals characters beyond the point
// Must calculate rounding value as separate variable to avoid mystery error

uint64_t RoundFixedPt(union ll_u TheNumber,unsigned Decimals) {
union ll_u Rnd;

//  printf(" round before: %08lx %08lx\n",TheNumber.fx_32.high,TheNumber.fx_32.low);

  Rnd.fx_64 = (One.fx_64 / 2) / (pow(10LL,Decimals));
//  printf("         incr: %08lx %08lx\n",Rnd.fx_32.high,Rnd.fx_32.low);

  TheNumber.fx_64 = TheNumber.fx_64 + Rnd.fx_64;
//  printf("        after: %08lx %08lx\n",TheNumber.fx_32.high,TheNumber.fx_32.low);

  return TheNumber.fx_64;
}


//-----------
// Multiply two unsigned scaled fixed point numbers without overflowing a 64 bit value
// The product of the two integer parts mut be < 2^32

uint64_t MultiplyFixedPt(union ll_u Mcand, union ll_u Mplier) {
union ll_u Result;

  Result.fx_64  = ((uint64_t)Mcand.fx_32.high * (uint64_t)Mplier.fx_32.high) << 32; // integer parts (clear fract)

  Result.fx_64 += ((uint64_t)Mcand.fx_32.low * (uint64_t)Mplier.fx_32.low) >> 32;   // fraction parts (always < 1)

  Result.fx_64 += (uint64_t)Mcand.fx_32.high * (uint64_t)Mplier.fx_32.low;          // cross products

  Result.fx_64 += (uint64_t)Mcand.fx_32.low * (uint64_t)Mplier.fx_32.high;

  return Result.fx_64;
}


//-----------
// Long long print-to-buffer helpers
//  Assumes little-Endian layout

void PrintHexLL(char *pBuffer,union ll_u FixedPt) {
  sprintf(pBuffer,"%08lx %08lx",FixedPt.fx_32.high,FixedPt.fx_32.low);
}

// converts all 9 decimal digits of fraction, which should suffice

void PrintFractionLL(char *pBuffer,union ll_u FixedPt) {
  union ll_u Fraction;

  Fraction.fx_64 = FixedPt.fx_32.low;             // copy 32 fraction bits, high order = 0
  Fraction.fx_64 *= GIGA;                         // times 10^9 for conversion
  Fraction.fx_64 >>= 32;                          // align integer part in low long
  sprintf(pBuffer,"%09lu",Fraction.fx_32.low);    // convert low long to decimal
}

void PrintIntegerLL(char *pBuffer,union ll_u FixedPt) {
  sprintf(pBuffer,"%lu",FixedPt.fx_32.high);
}

void PrintFixedPt(char *pBuffer,union ll_u FixedPt) {
  PrintIntegerLL(pBuffer,FixedPt);  // do the integer part
  pBuffer += strlen(pBuffer);       // aim pointer beyond integer
  *pBuffer++ = '.';                 // drop in the decimal point, tick pointer
  PrintFractionLL(pBuffer,FixedPt);
}

void PrintFixedPtRounded(char *pBuffer,union ll_u FixedPt,unsigned Decimals) {
char *pDecPt;
//char *pBase;

//  pBase = pBuffer;

  FixedPt.fx_64 = RoundFixedPt(FixedPt,Decimals);

  PrintIntegerLL(pBuffer,FixedPt);  // do the integer part

//  printf(" Buffer int: [%s]\n",pBase);

  pBuffer += strlen(pBuffer);       // aim pointer beyond integer


  pDecPt = pBuffer;                 // save the point location
  *pBuffer++ = '.';                 // drop in the decimal point, tick pointer

  PrintFractionLL(pBuffer,FixedPt);

//  printf(" Buffer all: [%s]\n",pBase);

  if (Decimals == 0)
    *pDecPt = 0;                    // 0 places means discard the decimal point
  else
    *(pDecPt + Decimals + 1) = 0;   // truncate string to leave . and Decimals chars

//  printf(" Buffer end: [%s]\n",pBase);


}


//-- Helper routine for printf()

int s_putc(char c, FILE *t) {
  Serial.write(c);
}

//-----------

void setup ()
{
  Serial.begin (115200);
  fdevopen(&s_putc,0);                            // set up serial output for printf()

  Serial.println (F("DDS calculation exercise"));
  Serial.println (F("Ed Nisley - KE4ZNU - May 2017\n"));

// set up useful constants

  TempFX.fx_64 = -1;
  PrintFixedPt(Buffer,TempFX);
  printf("Max fixed point: %s\n",Buffer);

  One.fx_32.high = 1;                      // Set up 1.0, a very useful constant
  PrintFixedPt(Buffer,One);
  printf("\n1.0: %s\n",Buffer);

  Tenth.fx_64 = One.fx_64 / 10;             // Likewise, 0.1
  PrintFixedPt(Buffer,Tenth);
  printf("\n0.1: %s\n",Buffer);
  PrintFixedPtRounded(Buffer,Tenth,9);                    // show rounded value
  printf("0.1 to 9 dec: %s\n",Buffer);

  TestFreq.fx_64 = RoundFixedPt(Tenth,3);                 // show full string after rounding
  PrintFixedPt(Buffer,TestFreq);
  printf("0.1 to 3 dec: %s (full string)\n",Buffer);

  PrintFixedPtRounded(Buffer,Tenth,3);                    // show truncated string with rounded value
  printf("0.1 to 3 dec: %s (truncated string)\n",Buffer);

  CtPerHz.fx_64 = -1;                       // Set up 2^32 - 1, which is close enough
  CtPerHz.fx_64 /= 125 * MEGA;              // divide by nominal oscillator
  PrintFixedPt(Buffer,CtPerHz);
  printf("\nCt/Hz = %s\n",Buffer);

  printf("Rounding: \n");
  for (int d = 9; d >= 0; d--) {
    PrintFixedPtRounded(Buffer,CtPerHz,d);
    printf("     %d: %s\n",d,Buffer);
  }

  HzPerCt.fx_64 = 125 * MEGA;               // 125 MHz / 2^32, without actually shifting!
  PrintFixedPt(Buffer,HzPerCt);
  printf("\nHz/Ct: %s\n",Buffer);

  TenthHzCt.fx_64 = MultiplyFixedPt(Tenth,CtPerHz);     // 0.1 Hz as delta-phase count
  PrintFixedPt(Buffer,TenthHzCt);
  printf("\n0.1 Hz as ct: %s\n",Buffer);

  printf("Rounding: \n");
  for (int d = 9; d >= 0; d--) {
    PrintFixedPtRounded(Buffer,TenthHzCt,d);
    printf("           %d: %s\n",d,Buffer);
  }

// Try out various DDS computations

  TestFreq.fx_64 = One.fx_64 * (60 * KILO);               // set 60 kHz
  PrintFixedPt(Buffer,TestFreq);
  printf("\nTest frequency: %s\n",Buffer);
  PrintFixedPtRounded(Buffer,TestFreq,1);
  printf("         round: %s\n",Buffer);

  TestCount.fx_64 = MultiplyFixedPt(TestFreq,CtPerHz);    // convert to counts
  PrintFixedPt(Buffer,TestCount);
  printf("Delta phase ct: %s\n",Buffer);
  PrintFixedPtRounded(Buffer,TestCount,0);
  printf("  round to int: %s\n",Buffer);

  TestFreq.fx_64 += Tenth.fx_64;                          // set 60000.1 kHz
  PrintFixedPt(Buffer,TestFreq);
  printf("\nTest frequency: %s\n",Buffer);
  PrintFixedPtRounded(Buffer,TestFreq,1);
  printf("         round: %s\n",Buffer);

  TestCount.fx_64 = MultiplyFixedPt(TestFreq,CtPerHz);    // convert to counts
  PrintFixedPt(Buffer,TestCount);
  printf("Delta phase ct: %s\n",Buffer);
  PrintFixedPtRounded(Buffer,TestCount,0);
  printf("  round to int: %s\n",Buffer);

  TestFreq.fx_64 -= Tenth.fx_64 * 2;                      // set 59999.9 kHz
  PrintFixedPt(Buffer,TestFreq);
  printf("\nTest frequency: %s\n",Buffer);
  PrintFixedPtRounded(Buffer,TestFreq,1);
  printf("         round: %s\n",Buffer);

  TestCount.fx_64 = MultiplyFixedPt(TestFreq,CtPerHz);    // convert to counts
  PrintFixedPt(Buffer,TestCount);
  printf("Delta phase ct: %s\n",Buffer);
  PrintFixedPtRounded(Buffer,TestCount,0);
  printf("  round to int: %s\n",Buffer);

  TestFreq.fx_64 = (599999LL * One.fx_64) / 10;           // set 59999.9 kHz differently
  PrintFixedPt(Buffer,TestFreq);
  printf("\nTest frequency: %s\n",Buffer);
  PrintFixedPtRounded(Buffer,TestFreq,1);
  printf("         round: %s\n",Buffer);

  TestCount.fx_64 = MultiplyFixedPt(TestFreq,CtPerHz);    // convert to counts
  PrintFixedPt(Buffer,TestCount);
  printf("Delta phase ct: %s\n",Buffer);
  PrintFixedPtRounded(Buffer,TestCount,0);
  printf("  round to int: %s\n",Buffer);

  TempFX.fx_64 = RoundFixedPt(TestCount,0);               // compute frequency from integer count
  TestFreq.fx_64 = MultiplyFixedPt(TempFX,HzPerCt);
  PrintFixedPt(Buffer,TestFreq);
  printf("Int ct -> freq: %s\n",Buffer);
  PrintFixedPtRounded(Buffer,TestFreq,1);
  printf("         round: %s\n",Buffer);

  TestFreq.fx_64 = One.fx_64 * (10 * MEGA);               // set 10 MHz
  PrintFixedPt(Buffer,TestFreq);
  printf("\nTest frequency: %s\n",Buffer);
  PrintFixedPtRounded(Buffer,TestFreq,1);
  printf("         round: %s\n",Buffer);

  TestCount.fx_64 = MultiplyFixedPt(TestFreq,CtPerHz);    // convert to counts
  PrintFixedPt(Buffer,TestCount);
  printf("Delta phase ct: %s\n",Buffer);
  PrintFixedPtRounded(Buffer,TestCount,0);
  printf("  round to int: %s\n",Buffer);

  TempFX.fx_64 = RoundFixedPt(TestCount,0);               // compute frequency from integer count
  TestFreq.fx_64 = MultiplyFixedPt(TempFX,HzPerCt);
  PrintFixedPt(Buffer,TestFreq);
  printf("Int ct -> freq: %s\n",Buffer);
  PrintFixedPtRounded(Buffer,TestFreq,1);
  printf("         round: %s\n",Buffer);

  TestFreq.fx_64 = One.fx_64 * (10 * MEGA);               // set 10 MHz + 0.1 Hz
  TestFreq.fx_64 += Tenth.fx_64;
  PrintFixedPt(Buffer,TestFreq);
  printf("\nTest frequency: %s\n",Buffer);
  PrintFixedPtRounded(Buffer,TestFreq,1);
  printf("         round: %s\n",Buffer);

  TestCount.fx_64 = MultiplyFixedPt(TestFreq,CtPerHz);    // convert to counts
  PrintFixedPt(Buffer,TestCount);
  printf("Delta phase ct: %s\n",Buffer);
  PrintFixedPtRounded(Buffer,TestCount,0);
  printf("  round to int: %s\n",Buffer);

  TempFX.fx_64 = RoundFixedPt(TestCount,0);               // compute frequency from integer count
  TestFreq.fx_64 = MultiplyFixedPt(TempFX,HzPerCt);
  PrintFixedPt(Buffer,TestFreq);
  printf("Int ct -> freq: %s\n",Buffer);
  PrintFixedPtRounded(Buffer,TestFreq,1);
  printf("         round: %s\n",Buffer);

  TestFreq.fx_64 = One.fx_64 * (10 * MEGA);               // set 10 MHz - 0.1 Hz
  TestFreq.fx_64 -= Tenth.fx_64;
  PrintFixedPt(Buffer,TestFreq);
  printf("\nTest frequency: %s\n",Buffer);
  PrintFixedPtRounded(Buffer,TestFreq,1);
  printf("         round: %s\n",Buffer);

  TestCount.fx_64 = MultiplyFixedPt(TestFreq,CtPerHz);    // convert to counts
  PrintFixedPt(Buffer,TestCount);
  printf("Delta phase ct: %s\n",Buffer);
  PrintFixedPtRounded(Buffer,TestCount,0);
  printf("  round to int: %s\n",Buffer);

  TempFX.fx_64 = RoundFixedPt(TestCount,0);               // compute frequency from integer count
  TestFreq.fx_64 = MultiplyFixedPt(TempFX,HzPerCt);
  PrintFixedPt(Buffer,TestFreq);
  printf("Int ct -> freq: %s\n",Buffer);
  PrintFixedPtRounded(Buffer,TestFreq,1);
  printf("         round: %s\n",Buffer);

  Oscillator.fx_64 = One.fx_64 * (125 * MEGA);
  Serial.println("Oscillator tune: CtPerHz");
  PrintFixedPtRounded(Buffer,Oscillator,2);
  printf(" Oscillator: %s\n",Buffer);

  for (int i=-10; i<=10; i++) {
    OscOffset.fx_64 = i * One.fx_64;
    CtPerHz.fx_64 = 1LL << 63;
    CtPerHz.fx_64 /= (Oscillator.fx_64 + OscOffset.fx_64) >> 33;
    PrintFixedPt(Buffer,CtPerHz);
    printf(" %+3d -> %s\n",i,Buffer);
  }

}

//-----------

void loop () {

}

DDSCalcTest.txt

DDS calculation exercise
Ed Nisley - KE4ZNU - May 2017

Max fixed point: 4294967295.999999999

1.0: 1.000000000

0.1: 0.099999999
0.1 to 9 dec: 0.100000000
0.1 to 3 dec: 0.100499999 (full string)
0.1 to 3 dec: 0.100 (truncated string)

Ct/Hz = 34.359738367
Rounding:
     9: 34.359738368
     8: 34.35973837
     7: 34.3597384
     6: 34.359738
     5: 34.35974
     4: 34.3597
     3: 34.360
     2: 34.36
     1: 34.4
     0: 34

Hz/Ct: 0.029103830

0.1 Hz as ct: 3.435973831
Rounding:
           9: 3.435973832
           8: 3.43597383
           7: 3.4359738
           6: 3.435974
           5: 3.43597
           4: 3.4360
           3: 3.436
           2: 3.44
           1: 3.4
           0: 3

Test frequency: 60000.000000000
         round: 60000.0
Delta phase ct: 2061584.302070550
  round to int: 2061584

Test frequency: 60000.099999999
         round: 60000.1
Delta phase ct: 2061587.738044382
  round to int: 2061588

Test frequency: 59999.900000000
         round: 59999.9
Delta phase ct: 2061580.866096718
  round to int: 2061581

Test frequency: 59999.899999999
         round: 59999.9
Delta phase ct: 2061580.866096710
  round to int: 2061581
Int ct -> freq: 59999.914551639
         round: 59999.9

Test frequency: 10000000.000000000
         round: 10000000.0
Delta phase ct: 343597383.678425103
  round to int: 343597384
Int ct -> freq: 10000000.014506079
         round: 10000000.0

Test frequency: 10000000.099999999
         round: 10000000.1
Delta phase ct: 343597387.114398935
  round to int: 343597387
Int ct -> freq: 10000000.114506079
         round: 10000000.1

Test frequency: 9999999.900000000
         round: 9999999.9
Delta phase ct: 343597380.242451271
  round to int: 343597380
Int ct -> freq: 9999999.914506079
         round: 9999999.9
Oscillator tune: CtPerHz
 Oscillator: 125000000.00
 -10 -> 34.359741116
  -9 -> 34.359741116
  -8 -> 34.359740566
  -7 -> 34.359740566
  -6 -> 34.359740017
  -5 -> 34.359740017
  -4 -> 34.359739467
  -3 -> 34.359739467
  -2 -> 34.359738917
  -1 -> 34.359738917
  +0 -> 34.359738367
  +1 -> 34.359738367
  +2 -> 34.359737818
  +3 -> 34.359737818
  +4 -> 34.359737268
  +5 -> 34.359737268
  +6 -> 34.359736718
  +7 -> 34.359736718
  +8 -> 34.359736168
  +9 -> 34.359736168
 +10 -> 34.359735619

More details on my blog at http://wp.me/poZKh-6MX