% E04_6.m            28apr17   06mar23
% sinusoidal oscillator
clear all

A = 1.3;  ome0 = .093*pi;
ze = [];  po = [exp(j*ome0)  exp(-j*ome0)]';
% ze = [];  po = [1.02*exp(j*ome0)  1.02*exp(-j*ome0)]';
[b,a] = zp2tf(ze,po,1);
% fvtool(b,a);      % test the system
return

%%
figure(2);
y = A * sin(ome0*[0:199]);
subplot(511); plot(y);       % desired signal

x = [A*sin(ome0), zeros(1,199)];
y = filter(1,a,x);
subplot(512); plot(y);       % simulation

clear y;
y(1) = - A * sin(ome0);
y(2) = 0;
for n = 3:200
    y(n) = 2*cos(ome0)*y(n-1) - y(n-2);
    end;
subplot(513); plot(y);       % realization

yn(1) = - A * sin(ome0);
yn(2) = 0;
for n = 3:200
    eps = (rand(1,1)-0.5)*.001;  % additive uniform noise
    yn(n) = 2*cos(ome0)*yn(n-1) + eps - yn(n-2);
    end;
subplot(514); plot(yn);       % realization w. quantization error
subplot(515); plot(yn-y);

return;

