Low-overhead high quality sine wave generation

Resonant filters are often used for this purpose,
but with a remark that they need some way of amplitude
control, because the numeric errors can make it
rise without control (gain > 1) or decay to 0 (gain < 1).
I don't know about you, but for me it was a great
surprise to discover that for certain "magic" values
of amplitude the error integrated over the entire
period is *exactly* 0, which means that the filter
begins the next period in exactly the same state it
started, so no AGC is needed. For example, for
omega=2*PI/1024 and 24-bit arithmetic the stable
amplitude value is 2163013 (with two_cos = 0x7fff62).

I thought it might be an interesting result to share.

	Best regards, Piotr

I have used the following resonant filter oscillator structure for a long time. 

Integrator1 -> mult by sqrt(2-2cos(theta)) -> integrator2 -> mult by sqrt(2-2cos(theta)) ->  back to integrator1

Where integrator1 is delayed by 1 sample and integrator2 has no delay. You start it with an impulse and the amplitude is stable forever, even when you use finite wordlengths. 
However, if you change the frequency on-the-fly, the amplitude will shift. 
I forget the name of this structure but I think it was published in the 60's

Bob

radams2000@gmail.com wrote:

> Integrator1 -> mult by sqrt(2-2cos(theta)) -> integrator2 -> mult by sqrt(2-2cos(theta)) ->  back to integrator1

Interesting, I'll need to simulate it to evaluate its performance.
Especially that you didn't mention any magic constants to ensure 
stability. Mine is the much simpler, classic structure:

x[n] = two_cos*x[n-1] - x[n-2]

with x[-1] and x[-2] selected according to the desired initial phase. In 
my application one full sine cycle should happen after 1024 phase
increments in order to divide the grid frequency into that number of
intervals.

	Best regards, Piotr

On Sunday, September 24, 2017 at 4:41:09 AM UTC-4, Piotr Wyderski wrote:
> radams2000@gmail.com wrote:
> 
> > Integrator1 -> mult by sqrt(2-2cos(theta)) -> integrator2 -> mult by sqrt(2-2cos(theta)) ->  back to integrator1
> 
> Interesting, I'll need to simulate it to evaluate its performance.
> Especially that you didn't mention any magic constants to ensure 
> stability. Mine is the much simpler, classic structure:
> 
> x[n] = two_cos*x[n-1] - x[n-2]
> 
> with x[-1] and x[-2] selected according to the desired initial phase. In 
> my application one full sine cycle should happen after 1024 phase
> increments in order to divide the grid frequency into that number of
> intervals.

i used to think this method could conceivably either die off or go unstable because of roundoff error.  but the constant that determines if the sinusoid is decreasing or increasing is what multiplies x[n-2].  James McCartney straightened me out about this maybe two decades ago.

if you need a quadrature pair of oscillators, you can do it pretty simply but there *is* that concern of the sinusoids growing slowly (unstable) to dying off slowly.  there is a pretty simple fix to that to make the amplitude constant.

r b-j

Agreed. Even though I've used this trick in the past, when processing power was really low, I haven't used it recently. Just use a standard interpolated sine lookup table and you'll sleep better at night, plus you can get quadrature by just offsetting the pointers modulo N.  

Bob

Sat, 23 Sep 2017 11:07:56 +0200 : Piotr Wyderski:

> Resonant filters are often used for this purpose, but with a remark that
> they need some way of amplitude control, because the numeric errors can
> make it rise without control (gain > 1) or decay to 0 (gain < 1). I
> don't know about you, but for me it was a great surprise to discover
> that for certain "magic" values of amplitude the error integrated over
> the entire period is *exactly* 0, which means that the filter begins the
> next period in exactly the same state it started, so no AGC is needed.
> For example, for omega=2*PI/1024 and 24-bit arithmetic the stable
> amplitude value is 2163013 (with two_cos = 0x7fff62).
> 
> I thought it might be an interesting result to share.
> 
> 	Best regards, Piotr

You may consider alternatives:
- Table lookup and interpolation
- Approximation methods (Taylor power series, Chebychev polynomial)

BR, Alois

On Tuesday, September 26, 2017 at 8:50:55 AM UTC-4, radam...@gmail.com wrote:
> Agreed. Even though I've used this trick in the past, when processing power was really low, I haven't used it recently. Just use a standard interpolated sine lookup table and you'll sleep better at night, plus you can get quadrature by just offsetting the pointers modulo N.  
> 

it is the simplest and easiest method in the 99% of the cases where you have the memory for the LUT.

for a particular synthesizer for a particular synthesizer company in a particular town in Massachusetts i live in, i once did it by generating a simple sawtooth, turning that into a triangle, and running that through a 5th-order, odd-symmetry polynomial.

no memory for LUT, 2 states (only used one), 16 instructions (of a very quirky ASIC DSP) and good for 70 dB sine to harmonic ratio.  the other Bob was surprized i could do a 5th order with those resources.

- the 3rd Bob.

robert bristow-johnson wrote on 9/26/2017 1:43 PM:
> On Tuesday, September 26, 2017 at 8:50:55 AM UTC-4, radam...@gmail.com wrote:
>> Agreed. Even though I've used this trick in the past, when processing power was really low, I haven't used it recently. Just use a standard interpolated sine lookup table and you'll sleep better at night, plus you can get quadrature by just offsetting the pointers modulo N.
>>
>
> it is the simplest and easiest method in the 99% of the cases where you have the memory for the LUT.

The LUT is sized by the resolution required.  It works pretty well be using 
half the required phase bits for the address on the LUT and the remaining 
phase bits for the interpolation.  Unless you need a very high resolution 
phase word, the LUT isn't so large for adequate values of "so large".


> for a particular synthesizer for a particular synthesizer company in a particular town in Massachusetts i live in, i once did it by generating a simple sawtooth, turning that into a triangle, and running that through a 5th-order, odd-symmetry polynomial.

Am I correct in thinking that is at least five multiplies?  Or is it 10? 
Depending on your hardware and your sample rate that's a lot of computation.


> no memory for LUT, 2 states (only used one), 16 instructions (of a very quirky ASIC DSP) and good for 70 dB sine to harmonic ratio.  the other Bob was surprized i could do a 5th order with those resources.

I'm more of an FPGA guy and the resources for a 32 bit result are usually 
available with the LUT, interpolation technique.

-- 

Rick C

Viewed the eclipse at Wintercrest Farms,
on the centerline of totality since 1998

On Tuesday, September 26, 2017 at 6:20:17 PM UTC-4, rickman wrote:
> robert bristow-johnson wrote on 9/26/2017 1:43 PM:
> > On Tuesday, September 26, 2017 at 8:50:55 AM UTC-4, radam...@gmail.com wrote:
> >> Agreed. Even though I've used this trick in the past, when processing power was really low, I haven't used it recently. Just use a standard interpolated sine lookup table and you'll sleep better at night, plus you can get quadrature by just offsetting the pointers modulo N.
> >>
> >
> > it is the simplest and easiest method in the 99% of the cases where you have the memory for the LUT.
> 
> The LUT is sized by the resolution required.  It works pretty well be using 
> half the required phase bits for the address on the LUT and the remaining 
> phase bits for the interpolation.

well, that's how we normally do it.

>  Unless you need a very high resolution 
> phase word, the LUT isn't so large for adequate values of "so large".

i would think so.  this ASIC can do 64 voices and has a sampler built in (so it can do a sine-wave from the sampler), but it also has these little DSP blocks with 8 instructions, two inputs, one state, one output (which can also be a state) per block.  the instructions can simultaneously multiply, add, and do some other instruction like limit or rectify. but it didn't have a simple LUT.

 
> > for a particular synthesizer for a particular synthesizer company in a particular town in Massachusetts i live in, i once did it by generating a simple sawtooth, turning that into a triangle, and running that through a 5th-order, odd-symmetry polynomial.
> 
> Am I correct in thinking that is at least five multiplies?  Or is it 10? 
> Depending on your hardware and your sample rate that's a lot of computation.

you start with "x".

then you have x^2  (1 Mult)

and you have x^4  (2 Mults)

c_2 x^4   (3 Mults)

c_1 x^2   (4 Mults)

c_0 + c_1 x^2 + c_2 x^4   (addition)

x (c_0 + c_1 x^2 + c_2 x^4) (5 Mults)


> 
> Viewed the eclipse at Wintercrest Farms,
> on the centerline of totality since 1998

i was on the Cerulean-Hopkinsville Road in Kentucky.  about 50 feet from the point of greatest eclipse.

r b-j

robert bristow-johnson wrote on 9/26/2017 10:35 PM:
> On Tuesday, September 26, 2017 at 6:20:17 PM UTC-4, rickman wrote:
>> robert bristow-johnson wrote on 9/26/2017 1:43 PM:
>>> On Tuesday, September 26, 2017 at 8:50:55 AM UTC-4, radam...@gmail.com wrote:
>>>> Agreed. Even though I've used this trick in the past, when processing power was really low, I haven't used it recently. Just use a standard interpolated sine lookup table and you'll sleep better at night, plus you can get quadrature by just offsetting the pointers modulo N.
>>>>
>>>
>>> it is the simplest and easiest method in the 99% of the cases where you have the memory for the LUT.
>>
>> The LUT is sized by the resolution required.  It works pretty well be using
>> half the required phase bits for the address on the LUT and the remaining
>> phase bits for the interpolation.
>
> well, that's how we normally do it.
>
>>  Unless you need a very high resolution
>> phase word, the LUT isn't so large for adequate values of "so large".
>
> i would think so.  this ASIC can do 64 voices and has a sampler built in (so it can do a sine-wave from the sampler), but it also has these little DSP blocks with 8 instructions, two inputs, one state, one output (which can also be a state) per block.  the instructions can simultaneously multiply, add, and do some other instruction like limit or rectify. but it didn't have a simple LUT.

Specialized devices often are very good at the task they were designed for, 
but not so good at others.


>>> for a particular synthesizer for a particular synthesizer company in a particular town in Massachusetts i live in, i once did it by generating a simple sawtooth, turning that into a triangle, and running that through a 5th-order, odd-symmetry polynomial.
>>
>> Am I correct in thinking that is at least five multiplies?  Or is it 10?
>> Depending on your hardware and your sample rate that's a lot of computation.
>
> you start with "x".
>
> then you have x^2  (1 Mult)
>
> and you have x^4  (2 Mults)
>
> c_2 x^4   (3 Mults)
>
> c_1 x^2   (4 Mults)
>
> c_0 + c_1 x^2 + c_2 x^4   (addition)
>
> x (c_0 + c_1 x^2 + c_2 x^4) (5 Mults)

Maybe I'm missing something.  Maybe it's the odd symmetry.  Are you saying 
you don't need to calculate terms for x^2 and x^4 times appropriate 
coefficients?  What is the final expression you needed to calculate?  The 
last line shown would produce

c_0 x + c_1 x^3 + c_2 x^5

Isn't that right?  Or did you leave out the combining step that might produce

c_0 + c_0 x + c_1 x^2 + c_1 x^3 + c_2 x^4 + c_2 x^5


>> Viewed the eclipse at Wintercrest Farms,
>> on the centerline of totality since 1998
>
> i was on the Cerulean-Hopkinsville Road in Kentucky.  about 50 feet from the point of greatest eclipse.

The buildup was pretty cool and totality was awesome although short but the 
denouement was rather vapid.  Maybe we didn't have enough eclipse dancers.

-- 

Rick C

Viewed the eclipse at Wintercrest Farms,
on the centerline of totality since 1998