how to properly generate ideal sine test sample?

[brucehoult]

Here is a simple algorithm that will generate a sine-wave with 12.433075... samples per cycle.  The amplitude of the sine-wave will vary a little over time, but will never differ by more than 1ppm from its initial value, no matter how long it runs.  So more than good enough to test a 14 bit DAC, and exercise almost every code.

Code: [Select]

     x = 2079038281;
     a = 0;
     for (;;) {
          x += (a >> 1);
          a -= (x >> 1);
     }

'x' and 'a' are 32 bit signed integers. The DAC samples are the upper 14 bits of either 'x' or 'a' and have an amplitude of 8191. To avoid a 1/2 LSB DC offset in the output, add 0x20000 to the value in 'x' or 'a' before extracting the upper 14 bits.

this is a real magic 

How it works?

Nice trick there. The results are impressive, Cheap to implement in a FPGA too.

By the looks of it this is two integrators looped into each other. This sort of setup oscillates infinitely as the values move from one integrator to the other and back again (much like a ideal LC circuit with 0 damping). Similar to how a NOT gate ring oscillator works except this is in "analog values"

It's a rotation matrix, sneakily using the already-updated value of x for the -sin(theta) bit. Or two shearing matrix concatenated.

[radiolistener]

I see different initial value doesn't affects frequency, but affects amplitude in a non-obvious way.
How to select initial value for that algorithm?

Is there any way to change frequency but keep low distortions?

[cbutlera]

Here is a simple algorithm that will generate a sine-wave with 12.433075... samples per cycle.  The amplitude of the sine-wave will vary a little over time, but will never differ by more than 1ppm from its initial value, no matter how long it runs.  So more than good enough to test a 14 bit DAC, and exercise almost every code.

Code: [Select]

     x = 2079038281;
     a = 0;
     for (;;) {
          x += (a >> 1);
          a -= (x >> 1);
     }

'x' and 'a' are 32 bit signed integers. The DAC samples are the upper 14 bits of either 'x' or 'a' and have an amplitude of 8191. To avoid a 1/2 LSB DC offset in the output, add 0x20000 to the value in 'x' or 'a' before extracting the upper 14 bits.

this is a real magic 

How it works?

This is an algorithm that I discovered while playing around with numbers back in the 1980s.  I was writing DSP code for voice-band modems at the time, on quite limited hardware, so tricks like this looked interesting.

There is nothing special about the initial value for x.  It just sets the amplitude, although some choices can result in a little more variation in the amplitude with time than others.  There is also nothing special about the right shift.  It is just a 32 bit multiply by 0.5.

An example of the more general algorithm that can generate any arbitrary frequency with no significant long term amplitude drift is:

Code: [Select]

     p = 2 * sin(pi / n)
     x = 1;
     a = -p / 2;
     for (;;) {
          x += a * p;
          a -= x * p;
     }

Here x has an amplitude of 1, and n is the number of samples per cycle of the generated sine-wave.  The accumulated amplitude error over time is typically less than 1:10¹² after 10¹² samples, with all variables as IEEE 80-bit extended precision floating point.

The algorithm actually works with any fixed point or floating point number system.

Note that x and a are not orthogonal if regarded as simultaneous samples.  However, they are orthogonal as time interleaved samples, where the in-phase and quadrature DACs are alternately updated.

The key to the amplitude stability of this algorithm is that the frequency is set by a single constant, so the fact that this constant is of finite precision has no impact.  The only errors are those due to finite arithmetic, and those errors have no intrinsic bias to either increase or decrease the amplitude.  The accumulated amplitude error has the appearance of a random walk, so the area of the annulus containing the orthogonalised equivalent to x and a (let's call these x and y) only increases with the square root of the number of steps.  Therefore over time, the probability of the algorithm landing on an x, y pair that it has visited before increases.  Once this happens there will obviously be no further increase in the peak to peak variation of the amplitude.

Here is the maths (If I haven't fouled up the MathJax).

Assuming an amplitude of 1, a starting phase of zero and an angular increment of \$\theta\$, at iteration n we have:
$$
x_n=sin(n \theta) \\
y_n=cos(n \theta - \theta / 2)
$$
Using the standard expansion of \$cos(A+B)\$, \$y_n\$ can be written as
$$
y_n=cos(n \theta)cos(\theta/2 ) + sin(n \theta)sin(\theta/2)
$$
Similarly, \$y_{n+1}\$ can be written as:
$$
y_{n+1}=cos(n \theta + \theta / 2) \\
y_{n+1}=cos(n \theta)cos(\theta/2 ) - sin(n \theta)sin(\theta/2)
$$
So \$y_{n+1}\$ can therefore be calculated from \$y_n\$, \$x_n\$ and the constant \$2sin(\theta/2) \$:
$$
y_{n+1}=y_n - 2 sin(\theta/2)x_n
$$
Using the standard expansion of \$sin(A+B) \$, \$x_n \$ can be written as:
$$
x_n=sin(n \theta)=sin(n \theta + \theta / 2 - \theta / 2) \\
x_n=sin(n \theta + \theta / 2)cos(\theta / 2) - cos(n \theta + \theta / 2)sin(\theta / 2)
$$
Similarly, \$x_{n+1} \$ can be written as:
$$
x_{n+1}=sin(n \theta + \theta / 2 + \theta / 2) \\
x_{n+1}=sin(n \theta + \theta / 2)cos(\theta / 2) + cos(n \theta + \theta / 2)sin(\theta / 2)
$$
So \$x_{n+1} \$ can therefore be calculated from \$x_n \$, \$y_{n+1} \$ and the constant \$ 2sin(\theta/2) \$:
$$
x_{n+1}=x_n + 2sin(\theta/2)y_{n+1}
$$

[cbutlera]

I see different initial value doesn't affects frequency, but affects amplitude in a non-obvious way.
How to select initial value for that algorithm?

With reference to the more general algorithm.

Code: [Select]

     p = 2 * sin(pi / n)
     x = 1;
     a = -p / 2;
     for (;;) {
          x += a * p;
          a -= x * p;
     }

If a is initialised to 0 the initial value for x is the required amplitude divided by cos(pi/n).  If a is initialised as above, then the initial value for x is just the amplitude.

That set's the overall amplitude, but the precise initialisation value can then be adjusted a little by trial and error, to minimise the amplitude noise.

The precise value also changes the cycle time of the algorithm before it repeats the same sequence of a/x pairs.  In the case of the first example I posted, with a starting x of 2079038281, the cycle time is a little over 500 million samples.

[emece67]

.

[cbutlera]

AFAIK this algorithm is called the state-variable method. There was a recent thread here where it appeared https://www.eevblog.com/forum/programming/compact-sin(x)-which-gets-to-within-0-1-x-0-to-2xpi/msg3899078/#msg3899078

IMHO, its good behaviour here is due to its period (expressed as samples) being quite high. I'm referring here not to the period of the sampled sine, but the period of the sequence (how many samples it takes to repeat itself).

Thanks for finding that mention of the algorithm.  Do you know of a reference that investigates its stability?

The algorithm itself is almost trivial.  The interesting thing that I discovered back in 1980s was its long-term stability.  At the time, I assumed that what I had discovered was well known.  I'm still sure that I couldn't possibly have been the first to discover its long-term stability, but it would be nice to have a reference that dates from before the mid 1980s, so that I could know for sure.

After reading this thread, I tried to search for a reference to it myself, and all that I could find was a possible link to a 1996 paper by John Edwards https://www.embedded.com/whats-your-sine-finding-the-right-algorithm-for-digital-frequency-synthesis-on-a-dsp/

In the mid to late 1980s, John Edwards was a student on summer placement at the company where I was developing DSP code for voice-band modems.  He was assisting me with DSP code and algorithm modelling, so I would almost certainly have discussed what I had found with him.

EDIT:
Since writing the above, I have just been in touch with John Edwards https://www.numerix-dsp.com/.  It seems that the solution to the problem of error accumulation with recursive rotation that he proposed, and that was mentioned in the above link, is not the algorithm that I have been talking about.

[emece67]

.

[emece67]

.

[SiliconWizard]

Here is a simple algorithm that will generate a sine-wave with 12.433075... samples per cycle.  The amplitude of the sine-wave will vary a little over time, but will never differ by more than 1ppm from its initial value, no matter how long it runs.  So more than good enough to test a 14 bit DAC, and exercise almost every code.

Code: [Select]

     x = 2079038281;
     a = 0;
     for (;;) {
          x += (a >> 1);
          a -= (x >> 1);
     }

'x' and 'a' are 32 bit signed integers. The DAC samples are the upper 14 bits of either 'x' or 'a' and have an amplitude of 8191. To avoid a 1/2 LSB DC offset in the output, add 0x20000 to the value in 'x' or 'a' before extracting the upper 14 bits.

this is a real magic 

How it works?

Nice trick there. The results are impressive, Cheap to implement in a FPGA too.

By the looks of it this is two integrators looped into each other. This sort of setup oscillates infinitely as the values move from one integrator to the other and back again (much like a ideal LC circuit with 0 damping). Similar to how a NOT gate ring oscillator works except this is in "analog values"

It's a rotation matrix, sneakily using the already-updated value of x for the -sin(theta) bit. Or two shearing matrix concatenated.

Which in its given form, I also consider a particular case of CORDIC.

[Someone]

this is simply to test the idea that higher period sequences do behave better.

Longer sequences can have lower noise in any given bin, but its just processing gain:
https://ccrma.stanford.edu/~jos/sasp/Processing_Gain.html
That is not a guarantee the noise will be lower just by having a longer sequence.

The OP was trying to work out how to find the idealized rounding, there is not a direct solution. Rounding to integers usually introduces correlated noise, which does not disappear through processing gain.

[Someone]

It's a rotation matrix, sneakily using the already-updated value of x for the -sin(theta) bit. Or two shearing matrix concatenated.

Which in its given form, I also consider a particular case of CORDIC.

Expanding out the equation I didnt get the classical Toeplitz symmetry of a rotation matrix:
cos(a) -sin(a)
sin(a)  cos(a)

but instead, expanding the given code:

the more general algorithm.

Code: [Select]

     p = 2 * sin(pi / n)
     x = 1;
     a = -p / 2;
     for (;;) {
          x += a * p;
          a -= x * p;
     }


gives something decidedly not symmetric:
  1   p
-p  1-p^2

* sorry, cant be bothered to work out mathjax formatting!

[radiolistener]

despite the fact it has no harmonics, the noise (at noise floor) is very synthetic and contains a lot of a small spurs.
Is it possible to add some dithering in order to get "soft" and flat noise floor?

[Berni]

despite the fact it has no harmonics, the noise (at noise floor) is very synthetic and contains a lot of a small spurs.
Is it possible to add some dithering in order to get "soft" and flat noise floor?

14bit only has has a dynamic range of 86dB so the noise floor being down there in the -125dB range means it is performing many times better than you wanted in the original post.

[radiolistener]

14bit only has has a dynamic range of 86dB so the noise floor being down there in the -125dB range means it is performing many times better than you wanted in the original post.

Actually it has -111 dBFS noise floor at FFT size 8192.
So, there is process gain about 10*log10(8192/2)=36.12 dB and normalized SNR = 111-36 = 75 dBFS.

The picture above shows RMS noise for a small piece of bandwidth = 2.7 kHz, while ADC sample rate is 96000 kHz, so it shows noise level with process gain 10*log10(48000/2.7) = 42.5 dB. And normalized SNR is about 126-42.5 = 83.5 dBFS.

But I want to get rid of these spurs, to not confuse it with images from signal processing... its not an issue if the noise will be a little higher to hide these spurs.

[Berni]

You can sum in some random white noise to mask it using a larger flat noise floor.

Still this is a good result on actual hardware. Getting any significantly better noise floors out of it likely involves narrow bandpass filters to select just the frequency you want.