Matlab Sinusoidal Oscillator Implementations

This section provides Matlab/Octave program listings for the sinusoidal resonator/oscillator algorithms discussed above:

• Planar 2D Rotation (2DR) (``complex multiply'')
• Modified Coupled Form (MCF) (``Magic Circle'')
• Digital Waveguide Resonator (DWR)

The test program computes an impulse response of each resonator, and plots them overlaid for comparison. Next, it optionally plays each one as a sound and plots them individually.

% Filter test program in matlab

%N = 300;   % Number of samples to generate
N = 3000;   % Number of samples to generate
f = 100;   % Desired oscillation frequency (Hz)
fs = 8192; % Audio sampling rate (Hz)
%R = .99;   % Decay factor (1=never, 0=instant)
R = 1;   % Decay factor (1=never, 0=instant)
b1 = 1;    % Input gain to state variable x(n)
b2 = 0;    % Input gain to state variable y(n)

%nd = 16;   % Number of significant digits to use
nd = 4;   % Number of significant digits to use
base = 2;  % Mantissa base (two significant figures)
           % (see 'help chop' in Matlab)

u = [1,zeros(1,N-1)]; % driving input test signal
theta = 2*pi*f/fs; % resonance frequency, rad/sample

% ================================================
% 2D PLANAR ROTATION (COMPLEX MULTIPLY)

x1 = zeros(1,N); % Predeclare saved-state arrays
y1 = zeros(1,N);

x1(1) = 0;      % Initial condition
y1(1) = 0;      % Initial condition

c = chop(R*cos(theta),nd,base); % coefficient 1
s = chop(R*sin(theta),nd,base); % coefficient 2

for n=1:N-1,
  x1(n+1) = chop(   c*x1(n) - s*y1(n) + b1*u(n), nd,base);
  y1(n+1) = chop(   s*x1(n) + c*y1(n) + b2*u(n), nd,base);
end

% ================================================
% MODIFIED COUPLED FORM ("MAGIC CIRCLE")
% (ref: http://ccrma.stanford.edu/~jos/wgo/ )

x2 = zeros(1,N); % Predeclare saved-state arrays
y2 = zeros(1,N);

x2(1) = 0.0;     % Initial condition
y2(1) = 0.0;     % Initial condition

e = chop(2*sin(theta/2),nd,base); % tuning coefficient

for n=1:N-1,
  x2(n+1) = chop(R*(x2(n)-e*y2(n))+b1*u(n),nd,base);
  y2(n+1) = chop(R*(e*x2(n+1)+y2(n))+b2*u(n),nd,base);
end

% ================================================
% DIGITAL WAVEGUIDE RESONATOR (DWR)

x3 = zeros(1,N); % Predeclare saved-state arrays
y3 = zeros(1,N);

x3(1) = 0;      % Initial condition
y3(1) = 0;      % Initial condition

g = R*R;        % decay coefficient
t = tan(theta); % toward tuning coefficient
cp = sqrt(g/(g + t^2*(1+g)^2/4 + (1-g)^2/4)); % exact
%cp = 1 - (t^2*(1+g)^2 + (1-g)^2)/(8*g); % good approx
b1n = b1*sqrt((1-cp)/(1+cp)); % input scaling

% Quantize coefficients:
cp = chop(cp,nd,base);
g = chop(g,nd,base);
b1n = chop(b1n,nd,base);

for n=1:N-1,
  gx3 = chop(g*x3(n), nd,base); % mpy 1
  y3n = y3(n);
  temp = chop(cp*(gx3 + y3n), nd,base); % mpy 2
  x3(n+1) = temp - y3n + b1n*u(n);
  y3(n+1) = gx3 + temp + b2*u(n);
end

% ================================================
% playandplot.m

figure(4);
title('Impulse Response Overlay');
ylabel('Amplitude');
xlabel('Time (samples)');
alt=1;%alt = (-1).^(0:N-1); % for visibility
plot([y1',y2',(alt.*y3)']); grid('on');
legend('2DR','MCF','WGR');
title('Impulse Response Overlay');
ylabel('Amplitude');
xlabel('Time (samples)');

saveplot(sprintf('eps/tosc%d.ps',nd));

if 1
  playandplot(y1,u,fs,'2D rotation',1);
  playandplot(y2,u,fs,'Magic circle',2);
  playandplot(y3,u,fs,'WGR',3);
end

function playandplot(y,u,fs,ttl,fnum);
% sound(y,fs,16);
  figure(fnum);

  unwind_protect # protect graph state
  subplot(211);
  axis("labely");
  axis("autoy");
  xlabel("");
  title(ttl);
  ylabel('Amplitude');
  xlabel('Time (ms)');
  t = 1000*(0:length(u)-1)/fs;
  timeplot(t,u,'-');
  legend('input');
  subplot(212);
  ylabel('Amplitude');
  xlabel('Time (ms)');
  timeplot(t,y,'-');
  legend('output');
  unwind_protect_cleanup # restore graph state
    grid("off");
    axis("auto","label");
    oneplot();
  end_unwind_protect
endfunction

Note that the chop utility only exists in Matlab. In Octave, the following ``compatibility stub'' can be used:

function [y] = chop(x,nd,base)
y = x;