LFO (Part 2)

Hi,

I started a new thread because the posts take between 3 to 9 hours to
appear published in www.deja.com, which is what I use to participate
in comp.dsp.

It's related to the post about LFO implementation (in ADSP-21065L). I
wrote a simple m-file to check what's wrong with the original algo.
Here we go:

% LfoGen
%
% Generate oscillator
%
% Usage:
% [y,t] = lfoGen(Fo, Fs, tTot)
%   Fo  :   Oscillator frequency
%   Fs  :   Sampling frequency
%   tTot:   Total time
%
%   y   :   Output signal
%   t   :   Time signal

function [y,t] = lfoGen(Fo, Fs, tTot)

delta = 2*pi*Fo/Fs;
cosCoef = 2*cos(delta);
lastCos = 1.0;
currentCos = cosCoef;
t = (0:1/Fs:tTot)';
y = zeros(size(t));
y(1) = lastCos;
y(2) = currentCos;
for k = 3:length(y) 
    y(k) = cosCoef*y(k-1)-y(k-2);
end

figure, plot(t,y)

Then, for example, in the command window, let's generate a 1Hz LFO,
sampled at 44.1 kHz, running time is 4 secs:

>> y = lfoGen(1, 44100, 4);

The output signal is 1Hz, but with amplitude around 7000 (not 1) and
looks more like a sine signal than cosine (although phase info is no
very relevant for the app).

Let's try at 100 Hz:

>> y = lfoGen(100, 44100, 4);

Now, the signal is 100 Hz (ok) but amplitude is around 70 (not 1) and
looks more like a sine signal than cosine again (although phase info
is no very relevant for the app).

At 1kHz, the output frequency is effectively 1kHz, but amplitude is
around 7 (not 1). At 10kHz, the output frequency is effectively 10kHz,
but amplitude is around 1.011 (still not not 1!!!).

Could someone please help me to understand why, and correct it?????

Regards,

JaaC

Hello Jaime,

Thanks for the compliments on my paper. You have R. Lyons to thank, as he
was the one who got me to write it. The problem with extremely low frequency
osc comes from the differencing of nearly equal numbers so all of the bits
of precision are lost and the error accumulates. So you can implement the
osc by using two nearly equal in frequency quadrature oscillators (with equi
amplitude outputs) and then complex multiply one's output with the other.
For example with a sampling rate of 44.1kHz, pick -5500 Hz and 5502 Hz for
the two frequencies and you will get a 2 Hz output. The negative is used so
the complex multiply becomes a downshift in frequency. When the standard
quadrature osc operates at 1/8 the sample rate (stepping 45 degrees per
sample), the matrix becomes orthonormal which means its condition number is
one! The advantage of this is the error accumulates at the lowest rate. So
the amplitude stabilization can be done every once in a while as opposed to
every sample.  While this may seem like a lot of work, for some apps the two
oscs can be implemented with relatively low precision.

I hope this helps.
Clay


-- 
Clay S. Turner, V.P.
Wireless Systems Engineering, Inc.
Satellite Beach, Florida 32937
(321) 777-7889
www.wse.biz
csturner@wse.biz



"Jaime Andres Aranguren Cardona" <jaime.aranguren@ieee.org> wrote in message
news:14a86f87.0402281939.421e992e@posting.google.com...
> Hi,
>
> I started a new thread because the posts take between 3 to 9 hours to
> appear published in www.deja.com, which is what I use to participate
> in comp.dsp.
>
> It's related to the post about LFO implementation (in ADSP-21065L). I
> wrote a simple m-file to check what's wrong with the original algo.
> Here we go:
>
> % LfoGen
> %
> % Generate oscillator
> %
> % Usage:
> % [y,t] = lfoGen(Fo, Fs, tTot)
> %   Fo  :   Oscillator frequency
> %   Fs  :   Sampling frequency
> %   tTot:   Total time
> %
> %   y   :   Output signal
> %   t   :   Time signal
>
> function [y,t] = lfoGen(Fo, Fs, tTot)
>
> delta = 2*pi*Fo/Fs;
> cosCoef = 2*cos(delta);
> lastCos = 1.0;
> currentCos = cosCoef;
> t = (0:1/Fs:tTot)';
> y = zeros(size(t));
> y(1) = lastCos;
> y(2) = currentCos;
> for k = 3:length(y)
>     y(k) = cosCoef*y(k-1)-y(k-2);
> end
>
> figure, plot(t,y)
>
> Then, for example, in the command window, let's generate a 1Hz LFO,
> sampled at 44.1 kHz, running time is 4 secs:
>
> >> y = lfoGen(1, 44100, 4);
>
> The output signal is 1Hz, but with amplitude around 7000 (not 1) and
> looks more like a sine signal than cosine (although phase info is no
> very relevant for the app).
>
> Let's try at 100 Hz:
>
> >> y = lfoGen(100, 44100, 4);
>
> Now, the signal is 100 Hz (ok) but amplitude is around 70 (not 1) and
> looks more like a sine signal than cosine again (although phase info
> is no very relevant for the app).
>
> At 1kHz, the output frequency is effectively 1kHz, but amplitude is
> around 7 (not 1). At 10kHz, the output frequency is effectively 10kHz,
> but amplitude is around 1.011 (still not not 1!!!).
>
> Could someone please help me to understand why, and correct it?????
>
> Regards,
>
> JaaC

Hi J.
It seems to me that there are significant probs related to wordlength
precision for trig-based LF osc's.
I would be inclined to use a sine lookup table WITH interpolation.

The 'oscillator' is really a sawtooth ramp gen., followed by sine table
interpolation routine.  The frequency may be changed at any time by
adjusting variable 'step size'

    +-----------+
    | step size |-----+
    +-----------+     |
                    +---+   +-----------+          (RAMP)
                    |ADD|---|Index Limit|--------> 0.0 <= Ramp out < N
                    +---+   +-----------+    |      
    +-----------+     |                      V
+-->|Table Index|-----+                      |
|   +-----------+                            |
|                                            |
+--------------------------------------------+

Frequency accuracy is easily improved to any precision by using multiple
word fixed point arithmetic, you can generate a period of years.  Floating
point will give significant rounding errors when adding a small step size
of 0.001 to a much larger Table Index number like 1023.

If you make the number of entry points on the sine table large enough (>256)
the distortion quickly drops into the noise.  Creating an arbitrary size
sine table is relatively easy.  Runtime execution is fast and unaffected by
table size.  You may do this in fixed or floating point.

                         +----------+
                         |  Sine    |
                         |  Lookup  |
                         |  Table   |
                         |          |
(RAMP)   +---------+  I  |          |                         +-------+
---->-+--|Integer()|---->|  sin[I]  |--+----------------------|(+) Add|-->
      |  +---------+     |          |  |                   +--|(+)    |
      |                  |          |  |  +-------+        |  +-------+
      |                  |          |  +--|(-) Add|  Diff  |          Sine
      |                  | sin[I+1] |-----|(+)    |--+     |          O/P
      |                  |          |     +-------+  |     |
      |                  |          |                |     |
      |                  +----------+              +----+  |
      |                                            |Mult|--+ Scaled diff
      |                                            +----+
      |  +----------+  Frac                          |
      +--|Fraction()|--------------------------------+
         +----------+

Hope I got the pikky rite, working from soggy memory (mine is implemented in
unrefreshed dynamic RAM).

The Sine O/P will vary between -1 and approx +1. It is stable and unaffected
by frequency.  I have used this from below 0.001Hz up to Fnyquist.  It is
not really very smart, but I only want an oscillator to do two things:
1. Work
2. Work all the time.

Jim Adamthwaite

> Hello Jaime,
> 
> Thanks for the compliments on my paper. You have R. Lyons to thank, as he
> was the one who got me to write it. 

Sure! Rick, thank you very much for your efforts in keeping up the DSP
Tips & Tricks section in IEEE's Signal Processing Magazine (at least
up to now).

BTW, what has been the result of the process of getting the comments
from the readers of that sectin of the magazine, from comp.dsp, to
send to the Magazine editors ????

I think mine was a good example of how useful for real-life use by DSP
engineers is that section of the magazine. And for sure there are lots
of other stories with other articles, which we have not heared abou
yet. The good things must be said, too!

Regards,

JaaC

"Clay S. Turner" <CSTurner@WSE.Biz> wrote in message news:<0FJ0c.10915$6e7.432@bignews1.bellsouth.net>...
> [..] The problem with extremely low frequency
> osc comes from the differencing of nearly equal numbers so all of the bits
> of precision are lost and the error accumulates. So you can implement the
> osc by using two nearly equal in frequency quadrature oscillators (with equi
> amplitude outputs) and then complex multiply one's output with the other.
> For example with a sampling rate of 44.1kHz, pick -5500 Hz and 5502 Hz for
> the two frequencies and you will get a 2 Hz output. 

Hello, Clay.

But do you think that for example -5500 Hz and 5500.01 Hz it would
work? I know it theoretically should, but maybe due to that precision
issue, 5500.01 Hz becomes 5500 Hz and we don't achieve what we wanted
to. I'm just guessing, I don't know, but maybe you do.

> The negative is used so
> the complex multiply becomes a downshift in frequency. When the standard
> quadrature osc operates at 1/8 the sample rate (stepping 45 degrees per
> sample), the matrix becomes orthonormal which means its condition number is
> one! The advantage of this is the error accumulates at the lowest rate. So
> the amplitude stabilization can be done every once in a while as opposed to
> every sample.  While this may seem like a lot of work, for some apps the two
> oscs can be implemented with relatively low precision.

Mmhh... if you say so! Anyway, fortunately I was able to get it done
by means of implementing one quadrature oscillator at Fo as low as
0.01 Hz, with 32 bit fixed point arithmetic, and the 'correction
factor', and it seems to be working fine. It could be a good excercise
to try to do the same (0.01 Hz osc sampled at 44.1 kHz) in a 16bit
DSP, using your approach. Has anyone done that?

Regards,

JaaC

Hello Jaime,

comments below:



"Jaime Andres Aranguren Cardona" <jaime.aranguren@ieee.org> wrote in message
news:14a86f87.0403020527.c1972d1@posting.google.com...
> "Clay S. Turner" <CSTurner@WSE.Biz> wrote in message
news:<0FJ0c.10915$6e7.432@bignews1.bellsouth.net>...
> > [..] The problem with extremely low frequency
> > osc comes from the differencing of nearly equal numbers so all of the
bits
> > of precision are lost and the error accumulates. So you can implement
the
> > osc by using two nearly equal in frequency quadrature oscillators (with
equi
> > amplitude outputs) and then complex multiply one's output with the
other.
> > For example with a sampling rate of 44.1kHz, pick -5500 Hz and 5502 Hz
for
> > the two frequencies and you will get a 2 Hz output.
>
> Hello, Clay.
>
> But do you think that for example -5500 Hz and 5500.01 Hz it would
> work? I know it theoretically should, but maybe due to that precision
> issue, 5500.01 Hz becomes 5500 Hz and we don't achieve what we wanted
> to. I'm just guessing, I don't know, but maybe you do.

Certainly 0.01 Hz at 44.1kHz sampling poses some problems. The relative
error puts you at something like 19 bits below full scale. Moving out to the
sweet spot (frequency wise) helps but it won't solve this on a 16 bitter by
itself.

>
> > The negative is used so
> > the complex multiply becomes a downshift in frequency. When the standard
> > quadrature osc operates at 1/8 the sample rate (stepping 45 degrees per
> > sample), the matrix becomes orthonormal which means its condition number
is
> > one! The advantage of this is the error accumulates at the lowest rate.
So
> > the amplitude stabilization can be done every once in a while as opposed
to
> > every sample.  While this may seem like a lot of work, for some apps the
two
> > oscs can be implemented with relatively low precision.
>
> Mmhh... if you say so! Anyway, fortunately I was able to get it done
> by means of implementing one quadrature oscillator at Fo as low as
> 0.01 Hz, with 32 bit fixed point arithmetic, and the 'correction
> factor', and it seems to be working fine. It could be a good excercise
> to try to do the same (0.01 Hz osc sampled at 44.1 kHz) in a 16bit
> DSP, using your approach. Has anyone done that?
>
> Regards,
>
> JaaC

I'm glad the correction factor worked well here. I've used that trick for
years with much success. I'm still working on a two osc approach that can
reduce the precision requirements. Assuming infinite precision math, the two
oscillator solution can be written as a single osc. But I'm trying
topologies that reduce sensitivity to quantization errors. And so I'm
working on factoring a single osc into two where the error terms cancel out.
When I get some good reportable results I'll publish them.

Another approach to your problem would be to generate the sine with a much
lower sample rate and then use an FIR interpolation filter to upsample it.
How much harmonic distortion can you stand?

Of course another way that was mentioned by some is a simple phase
accumlator (a simple counter) followed by trigonometric evalution. Generally
today's DSPs are fast enough that most of these schemes are workable.

-- 
Clay S. Turner, V.P.
Wireless Systems Engineering, Inc.
Satellite Beach, Florida 32937
(321) 777-7889
www.wse.biz
csturner@wse.biz

On 2 Mar 2004 05:12:58 -0800, jaime.aranguren@ieee.org (Jaime Andres
Aranguren Cardona) wrote:

>> Hello Jaime,
>> 
>> Thanks for the compliments on my paper. You have R. Lyons to thank, as he
>> was the one who got me to write it. 
>
>Sure! Rick, thank you very much for your efforts in keeping up the DSP
>Tips & Tricks section in IEEE's Signal Processing Magazine (at least
>up to now).

Hi Jaime,
  you're most welcome.

>BTW, what has been the result of the process of getting the comments
>from the readers of that sectin of the magazine, from comp.dsp, to
>send to the Magazine editors ????
>
>I think mine was a good example of how useful for real-life use by DSP
>engineers is that section of the magazine. And for sure there are lots
>of other stories with other articles, which we have not heared abou
>yet. The good things must be said, too!
>
>Regards,
>
>JaaC

I've compiled and formatted the E-mails 
I received from you guys here on 
comp.dsp.  I've been so God-awful busy that,
as yet, I haven't composed the letter to the 
managers of the IEEEE Sig. Proc. magazine 
asking them to allocate more magazine pages 
to the "DSP Tips & Tricks" column.

Somehow, I have to convince them that working 
DSP engineers (practitioners) should be considered 
just as important to the IEEE as university professors. 

That letter will have to be *very* carefully 
written.  I hope to write it in the next month and 
have my Associate-Editor-partner, Dr. Amy Bell, 
review it.  Then I'll send the letter, and all the 
testimonials from you guys, to the 
magazine's managers (they're all university 
professors).

Take care,
[-Rick-]