Quadrature oscillator question

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.

If you expand the cos term by substituting for sin, you get:

cos = cos_last*(1-freq^2) + freq*sin_last

Questions:

1. Have you seen this form of oscillator before?  Does it go by any
particular name?

2. The notes I have say to set freq = (desired freq in Hz) *
2*pi/SampleRate.  This works at lower frequencies, but is off by
increasingly larger values near Nyquist.  It turns out the largest value you
can set for frequency is SampleRate/pi which actually generates a frequency
of SampleRate/2.
How can I change the formula for freq to linearize this so the generated
frequencies are correct?  I figured some sort of tangent function should do
it, but I can't figure it out.  The standard pre-warping equation doesn't
seem to be right (though it helps quite a bit).

"Jon Harris" <jon_harrisTIGER@hotmail.com> wrote in message
news:3f64b47d$1_1@newsfeed.slurp.net...
> 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.
>
> If you expand the cos term by substituting for sin, you get:
>
> cos = cos_last*(1-freq^2) + freq*sin_last
>
> Questions:
>
> 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

It looks, though, like yours is going backwards in time, unless you exchange
the signs.

> 2. The notes I have say to set freq = (desired freq in Hz) *
> 2*pi/SampleRate.  This works at lower frequencies, but is off by
> increasingly larger values near Nyquist.  It turns out the largest value
you
> can set for frequency is SampleRate/pi which actually generates a
frequency
> of SampleRate/2.
> How can I change the formula for freq to linearize this so the generated
> frequencies are correct?  I figured some sort of tangent function should
do
> it, but I can't figure it out.  The standard pre-warping equation doesn't
> seem to be right (though it helps quite a bit).

Well, it is only valid in the limit dt --> 0, so the smaller the better.

If it isn't small, not only will the frequency be wrong, but the amplitude
will grow exponentially.

-- glen

Thanks for the insight Glen, some comments below.

"Glen Herrmannsfeldt" <gah@ugcs.caltech.edu> wrote in message
news:TVo9b.462422$YN5.304813@sccrnsc01...
>
> "Jon Harris" <jon_harrisTIGER@hotmail.com> wrote in message
> news:3f64b47d$1_1@newsfeed.slurp.net...
> > 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.
> >
> > If you expand the cos term by substituting for sin, you get:
> >
> > cos = cos_last*(1-freq^2) + freq*sin_last
> >
> > Questions:
> >
> > 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).

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.
More on this important point below...

> It looks, though, like yours is going backwards in time, unless you
exchange
> the signs.

Could very well be.  Since it is just a tone generator, forwards/backwards
doesn't really matter.

> > 2. The notes I have say to set freq = (desired freq in Hz) *
> > 2*pi/SampleRate.  This works at lower frequencies, but is off by
> > increasingly larger values near Nyquist.  It turns out the largest value
> you
> > can set for frequency is SampleRate/pi which actually generates a
> frequency
> > of SampleRate/2.
> > How can I change the formula for freq to linearize this so the generated
> > frequencies are correct?  I figured some sort of tangent function should
> do
> > it, but I can't figure it out.  The standard pre-warping equation
doesn't
> > seem to be right (though it helps quite a bit).
>
> Well, it is only valid in the limit dt --> 0, so the smaller the better.

I have gotten around this up to this point by oversampling the sine
generation so that dt remains fairly small.

> If it isn't small, not only will the frequency be wrong, but the amplitude
> will grow exponentially.

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.

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.

If that's true, the real question is, is there a better oscillator I should
be using?  My requirements are as follows:

1. Is stable in amplitude and frequency over the audio frequency range
20-20k with a sample rate of 48K.
2. Has extremely low distortion/harmonic components (maybe 100dB down) when
implemented on a SHARC DSP
3. Is fairly computationally efficient on a SHARC DSP (ideally 25
instruction cycles or less) and doesn't require significant amounts of
memory (ideally < 10 words, so no LUTs).
4. The frequency can be swept rapidly and smoothly without having audible
glitches or needing to reset internal variables.

It requirement #4 that is usually problematic.  Other methods I've tried
work fine at a certain frequency, but don't allow for easy and clean
on-the-fly adjustments to the frequency.  Note that I DON'T need the
quadrature output, that just comes along for free with this technique.

Any ideas out there?

Jon Harris wrote:
> 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.
> 
> If that's true, the real question is, is there a better oscillator I should
> be using?  My requirements are as follows:
> 
> 1. Is stable in amplitude and frequency over the audio frequency range
> 20-20k with a sample rate of 48K.
> 2. Has extremely low distortion/harmonic components (maybe 100dB down) when
> implemented on a SHARC DSP
> 3. Is fairly computationally efficient on a SHARC DSP (ideally 25
> instruction cycles or less) and doesn't require significant amounts of
> memory (ideally < 10 words, so no LUTs).
> 4. The frequency can be swept rapidly and smoothly without having audible
> glitches or needing to reset internal variables.
> 
> It requirement #4 that is usually problematic.  Other methods I've tried
> work fine at a certain frequency, but don't allow for easy and clean
> on-the-fly adjustments to the frequency.  Note that I DON'T need the
> quadrature output, that just comes along for free with this technique.

In order to meet all of your requirements, I would recommend a DDS (one
number sets frequency, phase continuity on changes) followed by a polynomial
approximation to the sine function (no LUT, computationally efficient).
Add as many polynomial terms as needed to achieve the distortion spec.

FWIW, here is the quadrature oscillator code that I've used in the past
on a 2181-based project (not SHARC). IIRC, this had a distortion level of
about -78dB (measured on an Audio Precision), which was all I needed for
a fixed-frequency (1004 Hz) test tone in an audio transmission system.
If you port this to the SHARC, the additional precision that's available
might get you to the distortion levels you need, but there's still the
question of frequency agility.

Note that there's a hidden "trick" not mentioned in the comments: The
coefficients in sine_table and cosine_table have all been rounded up from
their exact values, never down. This insures that the sum of the squares
will always be slightly greater than one and the correction factor will
always be less than one, simplifying the code that applies the correction.

-- Dave Tweed

/*---------------------------------------------------------------------------*/

..var test_sine;
..var test_cosine;
..var tone_sine;
..var tone_cosine;
..var tone_sample_hi;
..var tone_sample_lo;

/* Three sets of coefficients to generate a 1004-Hz tone at 32, 44.1 and 48 kHz
 * sample rates. One pair of numbers is copied into test_sine and test_cosine
 * whenever the sample rate is changed.
 */
..var sine_table[3]   =  6418,  4671,  4294;
..var cosine_table[3] = 32133, 32433, 32485;

/* A recursion algorithm to implement a quadrature oscillator (sine and cosine
 * of same frequency and amplitude, with 90 degrees phase difference) The
 * equations used are:
 *   C(x+1) = C(x)cos(a) - S(x)sin(a) and
 *   S(x+1) = C(x)sin(a) + S(x)cos(a),
 * where a = 2*pi*desired_freq/sampling_freq
 * cos(a) and sin(a) are pre-computed constants for the three sample
 * frequencies. The initial conditions are C(0)=full scale, and S(0)=0.
 *
 * The disadvantage of this simple calculation is the error building up over
 * time. The creeping can be eliminated by using the MAC's saturation
 * capabilities. tw_sin and tw_cos are designed to be a full range +32767 to
 * -32768 signal. a_cos and a_sin are rounded UP so the creeping is in the
 * direction of a larger signal. To get the desired test output of -12 dBfs,
 * tone_sine is attenuated before output to the buffer.
 */
generate_tone:
        dis m_mode;
        mx0 = dm(test_cosine);
        mx1 = dm(test_sine);
        my0 = dm(tone_cosine);
        my1 = dm(tone_sine);

        mr = mx0 * my0 (ss);
        mr = mr - mx1 * my1 (ss);
        if mv sat mr;
        dm(tone_cosine) = mr1;

        mr = mx1 * my0 (ss);
        mr = mr + mx0 * my1 (ss);
        if mv sat mr;
        dm(tone_sine) = mr1;

        /* Compute sin**2 + cos**2, which should equal one.
         */
        mx0 = dm(tone_sine);
        my0 = dm(tone_sine);
        mx1 = dm(tone_cosine);
        my1 = dm(tone_cosine);
        mr = mx0 * my0 (ss);
        mr = mr + mx1 * my1 (ss);

        /* Correct the sine and cosine values by multiplying by the inverse
         * of this value. However, since the error is expected to be small,
         * we use the approximation 1/x = 2-x to compute the correction factor.
         * Note that the sum and the correction factor (its inverse) are
         * handled as unsigned numbers.
         */
        ay0 = 0;
        ar = ay0 - mr1;
        mr = ar * my1 (us);
        if mv sat mr;
        dm(tone_cosine) = mr1;

        mr = ar * my0 (us);
        if mv sat mr;
        dm(tone_sine) = mr1;

        /* Compute final output amplitude (-12dBfs)
         */
        my0 = 8231;
        mr = mr1 * my0 (rnd);
        dm(tone_sample_hi) = mr1;
        dm(tone_sample_lo) = mr0;

        rts;

/*---------------------------------------------------------------------------*/

"Jon Harris" <jon_harrisTIGER@hotmail.com> wrote in message
news:3f64db19$1_1@newsfeed.slurp.net...
> Thanks for the insight Glen, some comments below.
>
> "Glen Herrmannsfeldt" <gah@ugcs.caltech.edu> wrote in message
> news:TVo9b.462422$YN5.304813@sccrnsc01...
> >
> > "Jon Harris" <jon_harrisTIGER@hotmail.com> wrote in message
> > news:3f64b47d$1_1@newsfeed.slurp.net...
> > > 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.
> > >
> > > If you expand the cos term by substituting for sin, you get:
> > >
> > > cos = cos_last*(1-freq^2) + freq*sin_last
> > >
> > > Questions:
> > >
> > > 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.

> 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.

Which variables are integer, and which are 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.

-- glen

"Jon Harris" <jon_harrisTIGER@hotmail.com> wrote in message
news:3f64b47d$1_1@newsfeed.slurp.net...
> 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.
>
> If you expand the cos term by substituting for sin, you get:
>
> cos = cos_last*(1-freq^2) + freq*sin_last
>
> Questions:
>
> 1. Have you seen this form of oscillator before?  Does it go by any
> particular name?
>
> 2. The notes I have say to set freq = (desired freq in Hz) *
> 2*pi/SampleRate.  This works at lower frequencies, but is off by
> increasingly larger values near Nyquist.  It turns out the largest
value you
> can set for frequency is SampleRate/pi which actually generates a
frequency
> of SampleRate/2.
> How can I change the formula for freq to linearize this so the
generated
> frequencies are correct?  I figured some sort of tangent function
should do
> it, but I can't figure it out.  The standard pre-warping equation
doesn't
> seem to be right (though it helps quite a bit).
>
>
>

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

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

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

"Dave Tweed" <dtweed@acm.org> wrote in message
news:3F663AB8.9E6C2735@acm.org...
> Jon Harris wrote:
> > 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.
> >
> > If that's true, the real question is, is there a better oscillator I
should
> > be using?  My requirements are as follows:
> >
> > 1. Is stable in amplitude and frequency over the audio frequency range
> > 20-20k with a sample rate of 48K.
> > 2. Has extremely low distortion/harmonic components (maybe 100dB down)
when
> > implemented on a SHARC DSP
> > 3. Is fairly computationally efficient on a SHARC DSP (ideally 25
> > instruction cycles or less) and doesn't require significant amounts of
> > memory (ideally < 10 words, so no LUTs).
> > 4. The frequency can be swept rapidly and smoothly without having
audible
> > glitches or needing to reset internal variables.
> >
> > It requirement #4 that is usually problematic.  Other methods I've tried
> > work fine at a certain frequency, but don't allow for easy and clean
> > on-the-fly adjustments to the frequency.  Note that I DON'T need the
> > quadrature output, that just comes along for free with this technique.
>
> In order to meet all of your requirements, I would recommend a DDS (one
> number sets frequency, phase continuity on changes) followed by a
polynomial
> approximation to the sine function (no LUT, computationally efficient).
> Add as many polynomial terms as needed to achieve the distortion spec.

Thanks for the idea, Dave.  The DDS approach would certainly be able to deal
well with smooth frequency changes.  I'd have to see how many polynomial
terms I'd need to get the distortion to the -100dB or better level.  I'm
guesing quite a few.

> FWIW, here is the quadrature oscillator code that I've used in the past
> on a 2181-based project (not SHARC). IIRC, this had a distortion level of
> about -78dB (measured on an Audio Precision), which was all I needed for
> a fixed-frequency (1004 Hz) test tone in an audio transmission system.
> If you port this to the SHARC, the additional precision that's available
> might get you to the distortion levels you need, but there's still the
> question of frequency agility.

The frequency agility is usually the gotcha!

<snip code>

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