how to properly generate ideal sine test sample?

[radiolistener]

you're right, I missed the 1 difference since I was having a quick look, but for all intends and purposes that is HOW to generate a sine table. rounding, phase offset, scaling, etc. is the parameters you can play with if you want to optimize something

I already played with it a long time, and cannot find the way how version which starts from 0 3134 5792 is calculated. But this version has more clean sine than version calculated with simple round or without round.

Original sample has  -95.83 dBc for 3'rd harmonic. While round version has -87.91 dBc which is worse for almost 8 dB.

And now hamster_nz provided even more clean sine:

Code: [Select]

0  3134  5791  7566  8190  7566  5791  3134  0  -3134  -5791  -7566  -8190  -7566  -5791  -3134
it has third harmonic -105 dBc! 

The only problem is that it was generated with brute-force method. I'm looking the way to generate it with math expression.

Any idea why usual calculation with round produce result which is not the best?

[hamster_nz]

I'll leave the PC looking for better for a while, but this one seems to be quite good!

Code: [Select]

1598  4551  6811  8034  8034  6811  4551  1598 -1598 -4551 -6811 -8034 -8034 -6811 -4551 -1598

About 1/20th the  RMS error of the original. Be interested what your design makes of it.

(note that it is still scaled by 8181, but because of the phase the max value doesn't appear in the numbers...

[hamster_nz]

Any idea why usual calculation with round produce result which is not the best?

It is because of interplay between the rounding alters the level of the fundamental.

That and we are truly playing around it the weeds...

[radiolistener]

I'll leave the PC looking for better for a while, but this one seems to be quite good!

Code: [Select]

1598  4551  6811  8034  8034  6811  4551  1598 -1598 -4551 -6811 -8034 -8034 -6811 -4551 -1598

About 1/20th the  RMS error of the original. Be interested what your design makes of it.

(note that it is still scaled by 8181, but because of the phase the max value doesn't appear in the numbers...

Yes. All harmonics are bellow -103 dBc!   

[OM222O]

you're right, I missed the 1 difference since I was having a quick look, but for all intends and purposes that is HOW to generate a sine table. rounding, phase offset, scaling, etc. is the parameters you can play with if you want to optimize something

I already played with it a long time, and cannot find the way how version which starts from 0 3134 5792 is calculated. But this version has more clean sine than version calculated with simple round or without round.

Original sample has  -95.83 dBc for 3'rd harmonic. While round version has -87.91 dBc which is worse for almost 8 dB.

And now hamster_nz provided even more clean sine:

Code: [Select]

0  3134  5791  7566  8190  7566  5791  3134  0  -3134  -5791  -7566  -8190  -7566  -5791  -3134
it has third harmonic -105 dBc! 

The only problem is that it was generated with brute-force method. I'm looking the way to generate it with math expression.

Any idea why usual calculation with round produce result which is not the best?

To be honest I've always hated math and avoided it like the plague, but there's a much better method than brute force. For each offset (the variable) you can add up the decimal part of each sampling point (cost) and use something like stochastic gradient decent with momentum or some other optimization / approximation algorithm to minimize the cost. there are lots of examples of SGD with python online, so give those a shot. Brute force algorithm is bad because there are an infinite number of phase offsets you can try, but optimization algorithms converge fairly quickly. This minimizes quantization noise, but again, I'm not sure why you would want to do this.

[hamster_nz]

you're right, I missed the 1 difference since I was having a quick look, but for all intends and purposes that is HOW to generate a sine table. rounding, phase offset, scaling, etc. is the parameters you can play with if you want to optimize something

I already played with it a long time, and cannot find the way how version which starts from 0 3134 5792 is calculated. But this version has more clean sine than version calculated with simple round or without round.

Original sample has  -95.83 dBc for 3'rd harmonic. While round version has -87.91 dBc which is worse for almost 8 dB.

And now hamster_nz provided even more clean sine:

Code: [Select]

0  3134  5791  7566  8190  7566  5791  3134  0  -3134  -5791  -7566  -8190  -7566  -5791  -3134
it has third harmonic -105 dBc! 

The only problem is that it was generated with brute-force method. I'm looking the way to generate it with math expression.

Any idea why usual calculation with round produce result which is not the best?

To be honest I've always hated math and avoided it like the plague, but there's a much better method than brute force. For each offset (the variable) you can add up the decimal part of each sampling point (cost) and use something like stochastic gradient decent with momentum or some other optimization / approximation algorithm to minimize the cost. there are lots of examples of SGD with python online, so give those a shot. Brute force algorithm is bad because there are an infinite number of phase offsets you can try, but optimization algorithms converge fairly quickly. This minimizes quantization noise, but again,...

Talk is cheap - feel free to find the optimal answer . I know there is at least one better set of values out there for you to find, and I would be very interested in seeing your code.

... I'm not sure why you would want to do this.

Sometimes just playing with numbers is fun - and it is exploring a bit of the landscape where some of my own problems sometimes live.

It might come in useful one day, for exactly the same sort of reasons it is useful for Radiolistener. Then again it might not 

[Someone]

The only problem is that it was generated with brute-force method. I'm looking the way to generate it with math expression.

Any idea why usual calculation with round produce result which is not the best?

To be honest I've always hated math and avoided it like the plague, but there's a much better method than brute force. For each offset (the variable) you can add up the decimal part of each sampling point (cost) and use something like stochastic gradient decent with momentum or some other optimization / approximation algorithm to minimize the cost. there are lots of examples of SGD with python online, so give those a shot. Brute force algorithm is bad because there are an infinite number of phase offsets you can try, but optimization algorithms converge fairly quickly. This minimizes quantization noise, but again, I'm not sure why you would want to do this.

As already said:

If you want to minimize quantization error that's pretty easy to do mathematically.

Except it isnt! The cost function has many local minima, so there is no easy method.

If you'd like to propose a solution to this problem that has a convex solution we're all ears. But quantisation noise is still an unsolved optimisation problem, gradient descent fails hard.

[brucehoult]

The only problem is that it was generated with brute-force method. I'm looking the way to generate it with math expression.

That's mathematically impossible. The best you can do is search a state space systematically for good candidates, then do a 2D gradient descent on amplitude and phase from the best candidates to find their local optimums.

[OM222O]

The only problem is that it was generated with brute-force method. I'm looking the way to generate it with math expression.

Any idea why usual calculation with round produce result which is not the best?

To be honest I've always hated math and avoided it like the plague, but there's a much better method than brute force. For each offset (the variable) you can add up the decimal part of each sampling point (cost) and use something like stochastic gradient decent with momentum or some other optimization / approximation algorithm to minimize the cost. there are lots of examples of SGD with python online, so give those a shot. Brute force algorithm is bad because there are an infinite number of phase offsets you can try, but optimization algorithms converge fairly quickly. This minimizes quantization noise, but again, I'm not sure why you would want to do this.

As already said:

If you want to minimize quantization error that's pretty easy to do mathematically.

Except it isnt! The cost function has many local minima, so there is no easy method.

If you'd like to propose a solution to this problem that has a convex solution we're all ears. But quantisation noise is still an unsolved optimisation problem, gradient descent fails hard.

A very quick and dirty plot of the rms error with phase going from 0 to pi in 1000 steps shows a clear repeating pattern:


So realistically you only need to consider the first 100 steps or so, after that all of the answers are repeated and the same. I changed it to calculate with 10000 steps. 9358th value (phase offset of 2.9399024052293283 radians) had the lowest RMS error (0.5637751689365491) and the actual values were:

[ 1641 -1555 -4514 -6786 -8025 -8042 -6835 -4587 -1641  1555  4514  6786 8025  8042  6835  4587]

if you want to get an even better result, you can keep repeating this over a smaller and smaller range but I'll leave it here. The actual noise of the DAC and filter components will have a much higher impact in any real design, so this is all entirely pointless.

Edit: code for those who are curious

Code: [Select]

import numpy as np
from matplotlib import pyplot as plt

def create_table(gain, phase, points):
    idx = np.arange(points+1)
    sf = idx.max() / (2*np.pi)
    arr = idx[:-1]/sf
    arr = np.sin(arr + phase)*gain
    integer_arr = np.round(arr).astype('int32')
    return (arr, integer_arr)

def calculate_error(arr, integer_arr):
    diff = integer_arr - arr
    err = np.sqrt((diff**2).sum())
    return err

n = 10000
errors = []
increment = np.pi/n
for i in range(n):
    arr, integer_arr = create_table((2**13)-1, increment*i, 16)
    errors.append(calculate_error(arr, integer_arr))
errors = np.array(errors)
   
plt.figure(figsize=(16,9))
plt.plot(np.arange(len(errors)),errors)
plt.show()

arr, integer_arr = create_table((2**13)-1, increment*np.argmin(errors), 16)
print(np.argmin(errors))
print(increment*np.argmin(errors))
print(calculate_error(arr, integer_arr))
print(integer_arr)

If you want to run it online, remove the matplotlib import and the 3 plot lines (plt.figure, plt.plot and plt.show)

Edit 2: [ 1641  4587  6835  8042  8025  6786  4514  1555 -1641 -4587 -6835 -8042 -8025 -6786 -4514 -1555] is even better (100k steps over the first 0.4 radians)

[BrianHG]

... I'm not sure why you would want to do this.

Sometimes just playing with numbers is fun - and it is exploring a bit of the landscape where some of my own problems sometimes live.

It might come in useful one day, for exactly the same sort of reasons it is useful for Radiolistener. Then again it might not 

Here is a use: Using this sine table as a reference for multimode digital QAM modulator, every DB you can extract in purity especially with high quantization levels, the sturdier the SNR of your decoded data on the other side.

[SiliconWizard]

this results in:
[    0  3331  6087  7791  8147  7094  4815  1703 -1703 -4815 -7094 -8147 -7791 -6087 -3331     0]
you can easily change the number of samples , add phase offset or change the scaling factor etc.

2 sequential samples with the value 0 in the 16 phase loop?
Wouldn't that generate substantial distortion?

Yes, it's 16 samples while the period is computed over 15 samples. Looks like a typical off-by-one issue.

0 and 2 pi are the same point effectively. If this is an issue you can set the array to be of lenght N+1 and ignore the last value

Of course, but that means you only have 15 effective samples per period, and not 16.

[rstofer]

Edit: code for those who are curious

I was!  I am beginning to think there is some merit to Python, specifically in the libraries.  numpy and matplotlib are great examples.

[OM222O]

Edit: code for those who are curious

I was!  I am beginning to think there is some merit to Python, specifically in the libraries.  numpy and matplotlib are great examples.

For data analysis / scientific computing Python is the go-to language. It's extremely flexible and if performance is an issue for large computations, you can create C extensions for it using something like Cython. Numpy is actually written in C because it's used to crunch numbers.

I realized there is a small mistake in the code I posted. To calculate the RMS error, I didn't divide by N .The actual line should be:
 

Code: [Select]

err = np.sqrt((diff**2).sum()/len(arr))but this doesn't actually change the results.

[radiolistener]

Edit 2: [ 1641  4587  6835  8042  8025  6786  4514  1555 -1641 -4587 -6835 -8042 -8025 -6786 -4514 -1555] is even better (100k steps over the first 0.4 radians)

thanks, your version has -101.02 dBc for 3'rd harmonic.
And version generated with fixed code has -101.02 dBc for 3'rd harmonic.

By the way, hamster_nz provided me (in private message) with version which is even more clean: -109.13 dBc for 3'rd harmonic

[Someone]

Edit 2: [ 1641  4587  6835  8042  8025  6786  4514  1555 -1641 -4587 -6835 -8042 -8025 -6786 -4514 -1555] is even better (100k steps over the first 0.4 radians)

thanks, your version has -101.02 dBc for 3'rd and -105.96 dBc for 7'th harmonic.
And version generated with fixed code has -101.02 dBc for 3'rd and -105.97 dBc for 7'th harmonic.

By the way, hamster_nz provided me (in private message) with version which is even more clean: -109.13 dBc for 3'rd and -105.97 dBc for 7'th harmonic

Still not surpassing the earlier example:

[0 3107 5741 7501 8119 7501 5741 3107 0 -3107 -5741 -7501 -8119 -7501 -5741 -3107]

If you have a preference for maximum suppression in the 3rd harmonic it can be pushed at least -148dB (likely further).

[radiolistener]

If you have a preference for maximum suppression in the 3rd harmonic it can be pushed at least -148dB (likely further).

How? It will be useful for filter testing.

[radiolistener]

Still not surpassing the earlier example:

[0 3107 5741 7501 8119 7501 5741 3107 0 -3107 -5741 -7501 -8119 -7501 -5741 -3107]

wow, I didn't noticed this one. It is not full scale -0.08 dBFS, but its SFDR is 119.56 dB
How did you generated it?

PS: I was changed approach to find the harmonic index, previously I was used search local extremum around expected index, but such way is buggy, because it catch noise peaks near harmonic position. Now I just use N*index, where index is the first harmonic index. So, now it displays NaN dB for even harmonics (missing harmonic). The same value for even harmonics it displays for other sine examples provided it this topic. Odd harmonics levels with a new approach is almost the same.

[Someone]

Still not surpassing the earlier example:

[0 3107 5741 7501 8119 7501 5741 3107 0 -3107 -5741 -7501 -8119 -7501 -5741 -3107]

wow, I didn't noticed this one. It is not full scale -0.08 dBFS, but its SFDR is 119.56 dB
How did you generated it?

That one is relatively quick to find with a simple sweep of amplitude selecting the lowest RMS error, optimisation with phase or directly on the points finds the same answer.

[Berni]

If you are using this to test DAC then looping a single sine wave sample is the wrong way to do it.

This runs the converter trough the same cycle of codes, so it can give a false idea of the performance as it might hit codes that are particularly good or ones that are particularly bad.

The proper way requires generating a sinewave that is not divisible with the sample rate. In that case each wave of the sine wave passes trough different values, this exercises many more unique DAC codes. The spectral analysis can then be done over a window of 1000s of these sine waves to capture the average spectrum.

Generating the test data for this is quite a bit harder tho. You either have to hold a really long waveform in memory or calculate the sinewave on the fly. This involves things like CORDIC algorithms. I would generally just cheat here and use the Altera NCO IP Block that implements all of this (Including showing you a FFT of its performance while you are adjusting the design constraints)

[hamster_nz]

Still not surpassing the earlier example:

[0 3107 5741 7501 8119 7501 5741 3107 0 -3107 -5741 -7501 -8119 -7501 -5741 -3107]

wow, I didn't noticed this one. It is not full scale -0.08 dBFS, but its SFDR is 119.56 dB
How did you generated it?

That one is relatively quick to find with a simple sweep of amplitude selecting the lowest RMS error, optimisation with phase or directly on the points finds the same answer.

That is very impressive! I'm going to file that away,

it may be useful as a practical joke to play on somebody....  If I tested my FFT code on that vector and saw those results I would swear that it was broken!

[hamster_nz]

Here's a signal with low harmonics, and amplitude > +/-8191, but all samples within the range of +/- 8191.

Code: [Select]

1236,  4249,  6615,  7974,  8119,  7028,  4867,  1965, -1236, -4249, -6615, -7974, -8119, -7028, -4867, -1965

Might make an interesting way to look for odd math overflows in a DSP design.

[cbutlera]

If you are using this to test DAC then looping a single sine wave sample is the wrong way to do it.

This runs the converter trough the same cycle of codes, so it can give a false idea of the performance as it might hit codes that are particularly good or ones that are particularly bad.

The proper way requires generating a sinewave that is not divisible with the sample rate. In that case each wave of the sine wave passes trough different values, this exercises many more unique DAC codes. The spectral analysis can then be done over a window of 1000s of these sine waves to capture the average spectrum.

Generating the test data for this is quite a bit harder tho. You either have to hold a really long waveform in memory or calculate the sinewave on the fly. This involves things like CORDIC algorithms. I would generally just cheat here and use the Altera NCO IP Block that implements all of this (Including showing you a FFT of its performance while you are adjusting the design constraints)

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.

[radiolistener]

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?

[BrianHG]

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?

I think we covered this in an old MCU thread on how to generate a sine wave without any floating point or multiplication hardware.  After LUT and interpolation, we finally reached something similar to this, but with only 16 bit ints.

[Berni]

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"