# IIR oscillator produces stable output in MATLAB, unstable on chip, how do I tune it?

I'm trying to make an oscillator using IIR filter and I'm having problems tuning it so that it's stable on the platform I'm using (ARM Cortex M4F and CMSIS DSPLib's arm_biquad_cascade_df1_f32). Note that this is my first time doing this and that I don't know much about control theory, so I could be missing something very obvious here.

Looking through some textbooks, I found the following formula for a transfer function of an oscillator:

$$H(z)=\frac{Asin(\omega_0)}{1-(2rcos(\omega_0))z^{-1}+r^2z^{-2}}$$ where, as far as I understand it, A is the amplitude of the oscillator, $\omega_0$ is the normalized angular frequency, which is calculated as $\omega_0=\frac{2f_{real}}{f_{sampling}}$ and the r is parameter which sets the poles distance to the unit circle and decides if the system's output will be dampened, opposite of dampened (can't remember the word at the moment) or constant.

I designed an oscillator using following parameters I got from some other requirements: sampling frequency $F_s=759493.7$, oscillator frequency $f1=25630$, amplitude $A=1$ (I plan to fine-tune this later), $r=1$.

This gives me then $\omega_0=\frac{2f_1}{F_s}=0.0674923$ and the transfer function of: $$H(z)=\frac{0.0674411}{1-1.9954464z^{-1}+1z^{-2}}$$

Looking at the impulse response in MATLAB, I get the expected result:

So far so good. I then tried to implement this using arm_biquad_cascade_df1_f32. Here's the block diagram of the filter according to documentation:
Since I only have a single second order section in my oscillator, I'm using only one section as shown on the picture. I set the coefficients to be:
$b_0=0.0674923$, $b_1=0$, $b_2=0$
$a_1=-1.9954464$, $a_2=1$

When I start executing this filter, the output is unstable and quickly goes to infinity. Manually experimenting with the parameter $r$, I managed to get a slowly decaying response with $r=0.5$. I suspect that what I'm seeing are the results of the MCU being only able to work with 32bit floating point numbers and lacking precision.

So my question is: How do I model the oscillator while taking into account the numerical precision errors?

-
Also any hints on how to improve this question are very welcome! –  AndrejaKo May 12 at 10:20
actually, $A$ is not the output level of the oscillator, but is a transfer gain from input to output. in fact, for an oscillator, there is no input and $b_0$, $b_1$, and $b_2$ are all zero. what gets the oscillator going is that your initial states are not zero. –  robert bristow-johnson May 12 at 16:54
@robert bristow-johnson Interesting! I haven't analyzed this area before, so my understanding isn't as good as it should be. –  AndrejaKo May 12 at 17:06
check the sign on $a_1$ and $a_2$. i think you got it wrong. if they're added in the denominator of your transfer function, they are subtracted in the difference equation (but your diagram shows them as added). –  robert bristow-johnson May 12 at 19:06
@robert bristow-johnson Well my difference equation is $y[n]=(2cos(\omega_0)y[n-1]-y[n-2]+Asin(\omega_0)\delta[n]$. The diagram posted is from the documentation of the function and the function should do addition by default. From what I understand, if I need subtraction, I should enter negative coefficients and then let them be added. This should give me a negative value for $a_1$ and a positive value for $a_2$. –  AndrejaKo May 12 at 21:18

i'll try to simplify what you have here. so, being an oscillator, there is no input, only output. but you can think of it as a filter (with zero states) with an impulse as input. but that impulse essentially sets the states of your filter and the rest of the input is zero.

the feedback equation (from your diagram, setting $b_0$, $b_1$, and $b_2$ to $0$) is

$$y[n] \ = \ 2 \cos(\omega_0) \ y[n-1] \ - \ y[n-2]$$

$\omega_0 = 2 \pi \frac{f_0}{f_s}$ and $f_s$ is the sampling frequency.

this equation sets the frequency to $\omega_0$ but does not set the amplitude or phase. you set the amplitude and phase by the judiciously setting the initial values of the two states, $y[n-1]$ and $y[n-2]$. assuming your initial output sample is $y[0]$, your initial states are $y[-1]$ and $y[-2]$.

so if you want your output to be

$$y[n] = A \cos(\omega_0 n + \phi)$$

for a given $A$ and $\phi$, then you just walk the equation back to times $n=-1$ and $n=-2$ to get your initial states $y[-1]$ and $y[-2]$.

it's that easy.

-
BTW, if you have numerical problems with very small $\omega_0$, use this trig identity: $$\cos(\omega_0) = 1 - 2 \sin^2\left(\frac{\omega_0}{2} \right)$$ that changes the recursion equation to $$y[n] \ = \ 2y[n-1] \ - \ y[n-2] \ - \ 4 \sin^2\left(\frac{\omega_0}{2} \right) y[n-1]$$ this works pretty well for both fixed-point and single-precision floating-point implementations. –  robert bristow-johnson May 12 at 17:57

A first thing: I think there is a $\pi$ missing in your conversion between frequencies and angular frequency. Your impulse response shows a period of about 90, while your period should be about 29 samples - there's clearly something wrong here.

What you are seeing on your target hardware is due to the difference in precision between matlab's implicit use of double precision (64-bit), and the use of 32-bit for the routines you chose to use.

Even if this is not very apparent through the transfer function or the biquad structure, what this "filter" with no input does is simply rotate its initial conditions. With 32-bit precision, there is no guarantee that the magnitude of the input vector will be preserved by rotation - it can either grow (instability) or decay.

If you intend to use this as an oscillator (that is to say, if the input will always be null except for the initial conditions):

• Regularly - say every 100th sample - mess with your state variables to readjust the norm of the vector to 1.
• Or use another simpler and more stable technique like a phase accumulator and a sine lookup table.
-
See this answer for an efficient example on how to stabilize the state variables. dsp.stackexchange.com/a/9868/3997 –  Hilmar May 12 at 11:57
Thanks for the $\pi$ factor. It really does look like it's missing. I'm not so sure about the numerical problems anymore. I mean, the results go from around 1 to say around 10e6 in maybe 20 or 30 samples. Could that big effect be due to that? –  AndrejaKo May 12 at 15:00
No, it shouldn't diverge that fast in as few samples - the drift takes place over thousands of samples. This points to an implementation error ; maybe you are not calling arm_biquad_cascade_df1_f32 with the arguments in the right format or layout? –  pichenettes May 12 at 15:03
pichen, probably the 2nd-order critically-stable feedback filter is the simplest way to generate a sinusoid. more simple than LUT and phase-accumulator. now if the waveform is more "sophisticated", then the LUT and phase-accumulator are likely the simplest. –  robert bristow-johnson May 12 at 17:32
also, i must credit James McCartney of Supercollider, for straightening me out on this (i.e. a decade ago, i had the same misconception), but as long as $r=1$ or $a_2 = -1$ and you can represent $1$ perfectly in your numerical scheme (which shouldn't be too hard), the amplitude should neither grow nor decay in time. –  robert bristow-johnson May 12 at 18:30