Quadrature oscillator question

Thanks Al.  I read the article and it was very well-written.  Unfortunately,
with the precision, accuracy, and distortion I'm after, it looks like all
the described methods require more memory than my application can afford.

"Al Clark" <dsp@danvillesignal.com> wrote in message
news:Xns93F8536A77C48aclarkdanvillesignal@66.133.130.30...
> I wrote a tutorial article called "A Low Distortion, High Accuracy Signal
> Generator Using the ADSP-218x or ADSP-219x Families". It was published in
> the Analog Devices publication: DSP Connection last year. This article is
> available on our web site.
>
> Al Clark
> Danville Signal Processing, Inc.
> --------------------------------------------------------------------
> Purveyors of Fine DSP Hardware and other Cool Stuff
> Available at http://www.danvillesignal.com

"Glen Herrmannsfeldt" <gah@ugcs.caltech.edu> wrote in message
news:h9v9b.464987$YN5.314538@sccrnsc01...
>
> "Jon Harris" <jon_harrisTIGER@hotmail.com> wrote in message
> > > > 1. Have you seen this form of oscillator before?  Does it go by any
> > > > particular name?
> > >
> > > Well, it comes from d(sin(t))=cos(t) dt, and d(cos(t))= -sin(t) dt.
> > >
> > > sin(t+dt)=sin(t)+d(sin(t)) =sin(t)+cos(t) dt
> > >
> > > cos(t+dt)=cos(t)+d(cos(t))=cos(t)-sin(t) dt
> >
> > I might add that these can be derived from the standard angle sum
formulas
> > if you use the approximations cos(dt) = 1 and sin(dt) = dt (valid only
for
> > small dt and exact in the limit as dt->0).
>
> If you follow the definition of derivative that is exactly the way it will
> come out.

Right!

> > Minus signs not withstanding, your sin formula exactly matches what I'm
> > doing.  However, my cos formula uses the newly calculated or present
value
> > for sin whereas the one you list above is based on the past value for
> sin.3
> > More on this important point below...
>
> (snip)
>
> > I noticed if I use your formulas, then the amplitude does grow
> exponentially
> > and the higher the frequency the faster the exponential growth.
However,
> if
> > I use the formula for cos that uses the current value of sin instead of
> the
> > past value, it is much better behaved.  There is still the frequency
error
> > problem, and the amplitude is higher than expected for higher
frequencies
> > (though at least it is stable rather than exponentially increasing).
This
> > must be a key approximation that makes this technique work.
>
> Are you sure it doesn't do anything else in between, other than assign sin
> to sin_last, and cos to cos_last.   For one, if you wanted that you could
> assign directly to sin and cos, using sin and cos, without any sin_last
and
> cos_last.  Also, I think it still grows exponentially, somewhat depending
on
> f.

That's all it does.  I'm not sure I follow you--you definitely need to store
past values for sine and cos for this to work.

Take a look at Olli's post.  I think he cracked this nut and was able to put
it into more standard z-domain notation.  Here's what the algorithm looks
like re-written in more standard form (thanks Olli!):

  sin(t) = sin(t-1)-(freq*cos(t-1))                  (1)
  cos(t) = cos(t-1)+(freq*sin(t))                    (2)

> Which variables are integer, and which are floating point?

They are all floating point.

> > OK, it sounds like there may not be any simple way to linearize this.  I
> > suppose I could upsample this even further to stay further in the
"almost
> > linear" range, but that becomes computationally burdensome.
>
> I think the more usual one would be a phase accumulator and look-up table,
> but it depends somewhat on the allowed values and required precision.

Yes, for many applications that would be effective.  Unfortunately, I'm
after essentially arbitrary frequency resolution (or at least sub-Hz) and
distortion of -100dB or better, so this makes for a pretty big look-up table
and/or complex interpolation schemes.

Olli, you are the man!  I think you cracked the nut.  Applying your
suggested equation for frequency works great.  Also, you showed me how to
modify the algorithm to generate just a single tone without the quadrature,
so that may be useful as well.

One other question: do you know why the amplitude grows near the Nyquist
rate?  The amplitude is essentially constant at frequencies up to about 1/4
Nyquist.  Then it grows fairly rapidly, becoming about 30dB larger right
near Nyquist.  I'm doing everything floating-point with 40-bit precision.

Thanks for the insights,
-Jon

"Olli Niemitalo" <oniemita@mail.student.oulu.fi> wrote in message
news:Pine.GSO.4.58.0309161318280.25373@paju.oulu.fi...
> On Mon, 15 Sep 2003, Jon Harris wrote:
>
> > I have been using a quadrature oscillator that I _think_ I got from
> > Chamberlain's "Musical Applications of Microprocessors" book.  However,
I
> > don't have that reference available right now and am trying to figure
out
> > some problems with the oscillator.
> >
> > The oscillator is defined as follows:
> >
> > sin=sin_last-(freq*cos_last)
> > cos=cos_last+(freq*sin)
> > where sin_last and cos_last are the values of sin and cos from the
previous
> > sample period.
>
> It is easier to analyze the program if you write it as:
>
>   sin(t) = sin(t-1)-(freq*cos(t-1))                  (1)
>   cos(t) = cos(t-1)+(freq*sin(t))                    (2)
>
> Substituting sin(t) in (2) with its definition from (1) we get:
>
>   cos(t) = (1-freq^2)*cos(t-1) + freq*sin(t-1)       (3)
>
> We would like to have cos() defined in terms of cos() only. So we want
> to get rid of that sin(t-1). From (2) we can solve:
>
>   sin(t) = (cos(t)-cos(t-1)) / freq
>
> If we substitute t with t-1, we get the equivalent:
>
>   sin(t-1) = (cos(t-1)-cos(t-2)) / freq              (4)
>
> By substituting sin(t-1) in (3) with this definition, we have cos()
> defined only in terms of cos():
>
>   cos(t) = (2-freq^2)*cos(t-1) - cos(t-2)
>
> That's a recurrence relation from which we can go on to derive the
> transfer function and to find the poles. The transfer function is:
>
>                     ...
>   H(z) = --------------------------                  (5)
>          1 - (2-freq^2)*z^-1 + z^-2
>
> The transfer function of a typical oscillator has poles at e^(f*I)
> and e^(-f*I), where f is the oscillator frequency:
>
>                         ...
>   H(z) = ----------------------------------
>          (z^-1 - e^(f*I) * (z^-1 - e^(-f*I)
>
>                    ...
>        = -----------------------                     (6)
>          1 - 2*COS(f)*z^-1 + z^2
>
> That is delightfully close to what we have in (5)! If they are indeed the
> same transfer function, we can write:
>
>     2-freq^2 = 2*COS(f)
>  => freq = sqrt(2-2*COS(f))                          (7)
>
> If you want to give the frequency in Hz:
>
>   freq = Frequency * (2*pi / SampleRate)
>
> Please see if this works or if I missed a sign or something somewhere...
> :-)
>
> -olli

Hello Jon,

Try here:

http://personal.atl.bellsouth.net/p/h/physics/oscpaper.pdf

Clay




"Jon Harris" <jon_harrisTIGER@hotmail.com> wrote in message news:<3f66072b_2@newsfeed.slurp.net>...
> I couldn't find the article with a Google search, nor could I find Clay's
> home page.  Any pointers?
> 
> "Ian Buckner" <Ian_Buckner@agilent.com> wrote in message
> news:1063707218.615245@cswreg.cos.agilent.com...
> >
> >
> > Hi Jon:
> >
> > Have a look at Clay Turner's "Recursive Discrete Time Sinusoidal
> > Oscillators",
> > IEEE Signal Processing Mag, May 2003, p103-111. It is also available
> > on the
> > web, probably in Clay's home page.
> >
> > Regards
> >     Ian

"Jon Harris" <jon_harrisTIGER@hotmail.com> wrote in message
news:3f660a93$1_4@newsfeed.slurp.net...
> "Glen Herrmannsfeldt" <gah@ugcs.caltech.edu> wrote in message
> news:h9v9b.464987$YN5.314538@sccrnsc01...
> >
> > "Jon Harris" <jon_harrisTIGER@hotmail.com> wrote in message
> > > > > 1. Have you seen this form of oscillator before?  Does it go by
any
> > > > > particular name?
> > > >
> > > > Well, it comes from d(sin(t))=cos(t) dt, and d(cos(t))= -sin(t) dt.
> > > >
> > > > sin(t+dt)=sin(t)+d(sin(t)) =sin(t)+cos(t) dt
> > > >
> > > > cos(t+dt)=cos(t)+d(cos(t))=cos(t)-sin(t) dt

(snip)

> > > I noticed if I use your formulas, then the amplitude does
> > >grow exponentially and the higher the frequency the faster the
> > >exponential growth.
> > > However, if I use the formula for cos that uses the current
> > > value of sin instead of the past value, it is much better behaved.
> > > There is still the frequency error problem, and the amplitude
> > > is higher than expected for higher frequencies
> > > (though at least it is stable rather than exponentially increasing).
> > >This must be a key approximation that makes this technique work.
> >
> > Are you sure it doesn't do anything else in between, other than assign
sin
> > to sin_last, and cos to cos_last.   For one, if you wanted that you
could
> > assign directly to sin and cos, using sin and cos, without any sin_last
> and
> > cos_last.  Also, I think it still grows exponentially,
> > somewhat depending on f.
>
> That's all it does.  I'm not sure I follow you--you definitely need to
store
> past values for sine and cos for this to work.

What I meant was that you can just write:

 sin=sin-(freq*cos)
 cos=cos+(freq*sin)

instead of using sin_last and cos_last, and then assigning them later.

> Take a look at Olli's post.  I think he cracked this nut and was able to
put
> it into more standard z-domain notation.  Here's what the algorithm looks
> like re-written in more standard form (thanks Olli!):
>
>   sin(t) = sin(t-1)-(freq*cos(t-1))                  (1)
>   cos(t) = cos(t-1)+(freq*sin(t))                    (2)
>
> > Which variables are integer, and which are floating point?
>
> They are all floating point.

Oh, I was expecting mostly fixed point.  That changes it a little.

If you use the one where cos_last and sin_last are used.  (If you consider
the limit dt-->0, it doesn't matter either way.)

 sin=sin_last-(freq*cos_last)
 cos=cos_last+(freq*sin_last)

and look at the sin**2+cos**2, which should stay constant,

sin**2=sin_last**2-2*sin_last*(freq*cos_last)+(freq*cos_last)**2
cos**2=cos_last**2+2*cos_last*(freq*sin_last)+(freq*sin_last)**2

so sin**2+cos**2=(1+freq**2)(sin_last**2+cos_last**2)

If instead you use the original:

> sin=sin_last-(freq*cos_last)
> cos=cos_last+(freq*sin)

sin**2=sin_last**2-2*sin_last*(freq*cos_last)+(freq*cos_last)**2
cos**2=cos_last**2+2*cos_last*(freq*sin)+(freq*sin)**2

=cos_last**2+2*cos_last*freq*(sin_last-freq*cos_last)+(freq*sin_last-freq**2
*cos_last)**2

then  sin**2+cos**2=(sin_last**2+cos_last**2)
+(sin_last**2-cos_last**2)*freq**2
        -2*freq**3*sin_last*cos_last+freq**4*cos_last**2

on the average the =(sin_last**2-cos_last**2) term will be zero, though it
will sometimes be positive and sometimes negative.  The final two terms I
don't believe will cancel, especially as freq gets larger.

I had thought you were doing it in fixed point, I don't know how much
difference it makes.

(snip)
> > I think the more usual one would be a phase accumulator and look-up
table,
> > but it depends somewhat on the allowed values and required precision.
>
> Yes, for many applications that would be effective.  Unfortunately, I'm
> after essentially arbitrary frequency resolution (or at least sub-Hz) and
> distortion of -100dB or better, so this makes for a pretty big look-up
table
> and/or complex interpolation schemes.

In 32 bit fixed point arithmetic it should be pretty much sub-Hz, even 16
bit will be close.

Consider a fixed point number with the binary point on the left, so it is
only the fractional part, representing the phase.  That is, 0 to 1
represents 0 to 360 degrees or 0 to 2pi radians.  For each sample you add
the appropriate delta phase, and compute the sin of the new value.  At
44.1kHz the smallest phase increment is 1/65536, for 0.67Hz, and frequencies
increasing in 0.67Hz increments.

In 32 bits at 44.1kHz the smallest frequency and frequency increment are
about 1e-5 Hz.

Not knowing how much space you have available, or processor power, I think a
medium sized lookup table in linear interpolation gets to 100dB pretty
easily.

I still have a data book somewhere from the days of 256x8 ROMs, for a ROM
set to do sin() lookups.  One ROM does the high bits, two more to do linear
interpolation, and an adder to finish off the result.  64kx16 is a little
large, though, I agree.

A 32 bit accumulator allows fine frequency resolution, while only the high
16 bits are used for the output.

-- glen

On Tue, 16 Sep 2003, Jon Harris wrote:

> I have been using a quadrature oscillator that I _think_ I got from
> Chamberlain's "Musical Applications of Microprocessors" book.  However,
> I don't have that reference available right now and am trying to figure
> out some problems with the oscillator.

The oscillator is defined as:

  sin(t) = sin(t-1)-(freq*cos(t-1))                         (1)
  cos(t) = cos(t-1)+(freq*sin(t))                           (2)

On Tue, 16 Sep 2003, Jon Harris wrote:

> Olli, you are the man!  I think you cracked the nut.  Applying your
> suggested equation for frequency works great.

Ah great! :-) You know what, it simplifies some more for f = 0..Pi:

  freq = sqrt(2-2*COS(f))
       = 2*SIN(f/2)

> One other question: do you know why the amplitude grows near the Nyquist
> rate?  The amplitude is essentially constant at frequencies up to about
> 1/4 Nyquist.  Then it grows fairly rapidly, becoming about 30dB larger
> right near Nyquist.

Let's find out! We concluded that the cos() part was fine:

  cos(t) = COS(t*f)                                         (3)

The question then is, is sin() at the right phase to cos()? We would
expect that:

  sin(t) = SIN(t*f)

Recall that from (2) we can solve sin() in terms of cos():

  sin(t) = (cos(t)-cos(t-1)) / freq                    (4)

We can promptly see that sin() is of the same frequency as cos(), but
perhaps the phase is off. If we use the trigonometric definitions of
sin() and cos() and the latest definition of freq, we can easily confirm
that something is wrong by assigning some random values, say t=1 and f=Pi:

  SIN(t*f) = (COS(t*f)-COS(t*f-f)) / (SIN(f/2)*2)
   SIN(Pi) = (COS(Pi)-COS(0)) / (SIN(Pi/2)*2)
         0 = (-1 - 1) / (1*2)
         0 = -1

Hence, initializing the algorithm by cos(0)=1 and sin(0)=0 will fail as
sin() is not in quadrature phase with cos(). Plugging into (4) the
trigonometric definitions of cos() and freq hints that setting sin(0) =
SIN(f/2) gives the correct amplitude (the phase will still be off,
though). If this is the case, the error of setting sin(0)=0 is, as you had
found out, worse for higher frequencies.

-olli

A look up table with linear interpolation give a distortion of

3.35 - 40 * log10(N) [dB] where N is the total length of the table.

-100 dB distortion levels will need a full table of >260 values (ie 512).=


With a 16 bit DSP use for the phase increment and the interpolation and=20
the output double precision.

Ren=E9

other stuff snipped.

Hello Olli,

There is a much simpler way to see that this is not a quadrature oscillator.
Just write it in matrix form and see if the upper left and lower right
matrix elements aren't the same. Since for this example, they are not, then
this is not a quadrature osc.

Jon had determined the iteration as



[c' ]     [ 1-w^2        w  ] [ c ]
[   ]  =  [                 ] [   ]
[s' ]     [   -w         1  ] [ s ]

In the limit of small w, then the upper left and lower right elements become
the same so the relative phase between c and s tends to become 90 degrees,
but for non small w, it is not quadrature.

Clay







>
>   SIN(t*f) = (COS(t*f)-COS(t*f-f)) / (SIN(f/2)*2)
>    SIN(Pi) = (COS(Pi)-COS(0)) / (SIN(Pi/2)*2)
>          0 = (-1 - 1) / (1*2)
>          0 = -1
>
> Hence, initializing the algorithm by cos(0)=1 and sin(0)=0 will fail as
> sin() is not in quadrature phase with cos(). Plugging into (4) the
> trigonometric definitions of cos() and freq hints that setting sin(0) =
> SIN(f/2) gives the correct amplitude (the phase will still be off,
> though). If this is the case, the error of setting sin(0)=0 is, as you had
> found out, worse for higher frequencies.
>
> -olli

"Olli Niemitalo" <oniemita@mail.student.oulu.fi> wrote in message
news:Pine.GSO.4.58.0309170924100.27078@paju.oulu.fi...
>
> On Tue, 16 Sep 2003, Jon Harris wrote:
>
> > I have been using a quadrature oscillator that I _think_ I got from
> > Chamberlain's "Musical Applications of Microprocessors" book.  However,
> > I don't have that reference available right now and am trying to figure
> > out some problems with the oscillator.
>
> The oscillator is defined as:
>
>   sin(t) = sin(t-1)-(freq*cos(t-1))                         (1)
>   cos(t) = cos(t-1)+(freq*sin(t))                           (2)
>
> On Tue, 16 Sep 2003, Jon Harris wrote:
>
> > Olli, you are the man!  I think you cracked the nut.  Applying your
> > suggested equation for frequency works great.
>
> Ah great! :-) You know what, it simplifies some more for f = 0..Pi:
>
>   freq = sqrt(2-2*COS(f))
>        = 2*SIN(f/2)

Even better!  I've added this simplifcation to my code.

> > One other question: do you know why the amplitude grows near the Nyquist
> > rate?  The amplitude is essentially constant at frequencies up to about
> > 1/4 Nyquist.  Then it grows fairly rapidly, becoming about 30dB larger
> > right near Nyquist.
>
> Let's find out! We concluded that the cos() part was fine:
>
>   cos(t) = COS(t*f)                                         (3)
>
> The question then is, is sin() at the right phase to cos()? We would
> expect that:
>
>   sin(t) = SIN(t*f)
>
> Recall that from (2) we can solve sin() in terms of cos():
>
>   sin(t) = (cos(t)-cos(t-1)) / freq                    (4)
>
> We can promptly see that sin() is of the same frequency as cos(), but
> perhaps the phase is off. If we use the trigonometric definitions of
> sin() and cos() and the latest definition of freq, we can easily confirm
> that something is wrong by assigning some random values, say t=1 and f=Pi:
>
>   SIN(t*f) = (COS(t*f)-COS(t*f-f)) / (SIN(f/2)*2)
>    SIN(Pi) = (COS(Pi)-COS(0)) / (SIN(Pi/2)*2)
>          0 = (-1 - 1) / (1*2)
>          0 = -1
>
> Hence, initializing the algorithm by cos(0)=1 and sin(0)=0 will fail as
> sin() is not in quadrature phase with cos(). Plugging into (4) the
> trigonometric definitions of cos() and freq hints that setting sin(0) =
> SIN(f/2) gives the correct amplitude (the phase will still be off,
> though). If this is the case, the error of setting sin(0)=0 is, as you had
> found out, worse for higher frequencies.
>
> -olli

This didn't quite work out for me.  I was setting sin(0) = 1 and cos(0) = 0.
Experimentally I found that setting sin(0) = cos(f/2) and cos(0) = 0 gave me
the right amplitude.  Maybe that is equivalent to what you were suggesting,
but it does seems to work.

However, it sounds like for dynamic frequency changes I am going to have to
dynamically "fix" the state variables, which is a bit of a hassle.  My
current solution of running over-sampled seems to work fine because the
amplitude error becomes negligibly small.

Thanks again for the insights Olli (and Clay)!
-Jon

Clay,

After reading your paper, it looks like my oscillator doesn't exactly match
any of the ones you list.  It is very similar to the Equi-amplitude,
staggered update oscillator of Fig. 5, but it differs by a minus sign
(assuming I'm interpretting this correctly.).  When I change the signs to
make it match Fig. 5 it doesn't seem to oscillate at all.  So either I'm
doing something wrong or maybe there is an error in that figure?

To summarize, the equations I am using are:
  sin(t) = sin(t-1)-(freq*cos(t-1))                  (1)
  cos(t) = cos(t-1)+(freq*sin(t))                    (2)
This seems to work as expected.

When I change them to match Fig. 5 they become:
  sin(t) = sin(t-1)-(freq*cos(t-1))                  (1)
  cos(t) = cos(t-1)-(freq*sin(t))                    (2b)
This doesn't seem to oscillate.

Any thoughts?

-Jon

"Clay S. Turner" <physics@bellsouth.net> wrote in message
news:1ae4d220.0309180716.643d437@posting.google.com...
> other stuff snipped.
>
> Hello Olli,
>
> There is a much simpler way to see that this is not a quadrature
oscillator.
> Just write it in matrix form and see if the upper left and lower right
> matrix elements aren't the same. Since for this example, they are not,
then
> this is not a quadrature osc.
>
> Jon had determined the iteration as
>
>
>
> [c' ]     [ 1-w^2        w  ] [ c ]
> [   ]  =  [                 ] [   ]
> [s' ]     [   -w         1  ] [ s ]
>
> In the limit of small w, then the upper left and lower right elements
become
> the same so the relative phase between c and s tends to become 90 degrees,
> but for non small w, it is not quadrature.
>
> Clay
>
>
>
>
>
>
>
> >
> >   SIN(t*f) = (COS(t*f)-COS(t*f-f)) / (SIN(f/2)*2)
> >    SIN(Pi) = (COS(Pi)-COS(0)) / (SIN(Pi/2)*2)
> >          0 = (-1 - 1) / (1*2)
> >          0 = -1
> >
> > Hence, initializing the algorithm by cos(0)=1 and sin(0)=0 will fail as
> > sin() is not in quadrature phase with cos(). Plugging into (4) the
> > trigonometric definitions of cos() and freq hints that setting sin(0) =
> > SIN(f/2) gives the correct amplitude (the phase will still be off,
> > though). If this is the case, the error of setting sin(0)=0 is, as you
had
> > found out, worse for higher frequencies.
> >
> > -olli