How to implement a digital oscillator?

Question

I have a floating-point digital signal processing system that operates at a fixed sample rate of $f_s = 32768$ samples per second implemented using an x86-64 processor. Assuming that the DSP system is synchronously locked to whatever matters, what is the best way to implement a digital oscillator at some frequency $f$?

Specifically, I want to generate the signal: $$y(t) = \sin(2\pi f t)$$ where $t=n/f_s$ for sample number $n$.

One idea is to keep track of a vector $(x,y)$ which we rotate by an angle $\Delta\phi = 2\pi f/f_s$ on each clock cycle.

As a Matlab pseudocode implementation (the real implementation is in C):

%% Initialization code

f_s = 32768;             % sample rate [Hz]
f = 19.875;              % some constant frequency [Hz]

v = [1 0];               % initial condition     
d_phi = 2*pi * f / f_s;  % change in angle per clock cycle

% initialize the rotation matrix (only once):
R = [cos(d_phi), -sin(d_phi) ; ...
     sin(d_phi),  cos(d_phi)]

Then, on each clock cycle, we rotate the vector around a little bit:

%% in-loop code

while (forever),
  v = R*v;        % rotate the vector by d_phi
  y = v(1);       % this is the sine wave we're generating
  output(y);
end

This allows the oscillator to be computed with only 4 multiplications per cycle. However, I'd worry about phase error and amplitude stability. (In simple tests I was surprised that the amplitude did not die or explode immediately--maybe the sincos instruction guarantees $\sin^2+\cos^2=1$?).

What is the right way to do this?

Answer 1

You're right that the strictly recursive approach is vulnerable to the accumulation of error as the number of iterations increases. One more robust way this is typically done is to use a numerically-controlled oscillator (NCO). Basically, you have an accumulator that keeps track of the instantaneous phase of the oscillator, updated as follows:

$$ \delta = \frac{2\pi f}{f_s} $$

$$ \phi[n] = (\phi[n-1] + \delta) \mod 2\pi $$

At each time instant, then, you're left with converting the accumulated phase in the NCO to the desired sinusoidal outputs. How you do this depends on your requirements for computational complexity, accuracy, etc. One obvious way is to just calculate the outputs as

$$ x_c[n] = \cos(\phi[n]) $$

$$ x_s[n] = \sin(\phi[n]) $$

using whatever implementation of sine/cosine you have available. In high-throughput and/or embedded systems, the mapping from phase to sine/cosine values is often done via a lookup table. The size of the lookup table (i.e. the amount of quantization you do on the phase argument to sine and cosine) can be used as a tradeoff between memory consumption and approximation error. The nice thing is that the amount of computations required is typically independent of the table size. In addition, you can limit your LUT size if needed by taking advantage of the symmetry inherent in the cosine and sine functions; you only really need to store one-fourth of a period of the sampled sinusoid.

If you need higher accuracy than a reasonably-sized LUT can give you, then you can always look at interpolation between table samples (linear or cubic interpolation, for example).

Another benefit of this approach is that it is trivial to incorporate frequency or phase modulation with this structure. The frequency of the output can be modulated by varying $\delta$ accordingly, and phase modulation can be implemented by simply adding to $\phi[n]$ directly.

Answer 2

What you have is a very good and efficient oscillator. The potential numerical drift problem can actually be solved. Your state variable v has two parts, one is eseentially the real part and the other the imaginary part. Let's call then r and i. We know that r^2+i^2 = 1. Over time this may drift up and down, however that can easily be corrected by multiplication with a gain correction factor like this $$g = \frac{1}{\sqrt{r^{2}+i^{2}}}$$

Obviously this is very expensive, however we know that the gain correction is very close to unity and we can approximate this with a simple Taylor expansion to $$g = \frac{1}{\sqrt{r^{2}+i^{2}}}\approx \frac{1}{2}\cdot \left ( 3-(r^{2}+i^{2}) \right )$$

Moreover we don't need to do this on every single sample, but once every 100 or 1000 samples is more than enough to keep this stable. This is particularly useful if you do frame based processing. Updating once per frame is just fine. Here is a quick Matlab calculates 10,000,000 samples.

%% seed the oscillator
% set parameters
f0 = single(100); % say 100 Hz
fs = single(44100); % sample rate = 44100;
nf = 1024; % frame size

% initialize phasor and state
ph =  single(exp(-j*2*pi*f0/fs));
state = single(1 + 0i); % real part 1, imaginary part 0

% try it
x = zeros(nf,1,'single');
testRuns = 10000;
for k = 1:testRuns
  % overall frames
  % sample: loop
  for i= 1:nf
    % phasor multiply
    state = state *ph;
    % take real part for cosine, or imaginary for sine
    x(i) = real(state);
  end
  % amplitude corrections through a taylor exansion aroud
  % abs(state) very close to 1
  g = single(.5)*(single(3)-real(state)*real(state)-imag(state)*imag(state) );
  state = state*g;
end
fprintf('Deviation from unity amplitude = %f\n',g-1);

Answer 3

You can avoid the unstable magnitude drift if you do not make it recursively update the vector v. Instead, rotate your prototype vector v to the current output phase. This still requires some trig functions, but only once per buffer.

No magnitude drift and arbitrary frequency

pseudocode:

init(freq)
  precompute Nphasor samples in phasor
  phase=0

gen(Nsamps)
    done=0
    while done < Nsamps:
       ndo = min(Nsamps -done, Nphasor)
       append to output : multiply buf[done:done+ndo) by cexp( j*phase )
       phase = rem( phase + ndo * 2*pi*freq/fs,2*pi)
       done = done+ndo

You can do away with the multiply, the trig functions required by cexp, and the modulus remainder over 2pi if you can tolerate a quantized frequency translation. e.g. fs/1024 for a 1024 sample phasor buffer.