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

Question

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?

Answer 1

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.

Answer 2

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.