Vicanek's quadrature oscillator

quadoscil.m

close all;
clear all;

N = 1e7;
f = 1e-2;

times=100;

y=zeros(1,N);
z=zeros(1,N);

y(1) = 1;
z(1) = 0;

errork1 = 1e-5;
errork2 = 1e-6;
r1 = -1 + 2*rand(1);
r2 = -1 + 2*rand(1);

errorcomp = 1e-6;

k1 = tan(f/2);
%k1 = k1*(1+r1*errork1);
k1 = k1+r1*errork1;
k2 = 2*k1/(1+k1^2);
%k2 = k2*(1+r2*errork2);
k2 = k2+r2*errork2;

r = -1 + 2*rand(3,N-1);

for rep=1:times
    for n=1:N-1
        t = y(n) - k1*z(n);
        %t = t*(1+errorcomp*r(1,n));
        t = t + errorcomp*r(1,n);
        z(n+1) = z(n) + k2*t;
        %z(n+1) = z(n+1)*(1+errorcomp*r(2,n));
        z(n+1) = z(n+1) + errorcomp*r(2,n);
        y(n+1) = t - k1*z(n+1);
        %y(n+1) = y(n+1)*(1+errorcomp*r(3,n));
        y(n+1) = y(n+1) + errorcomp*r(3,n);
    end
    if rep ~= times
        y(1) = y(N);
        z(1) = z(N);
    end
end

x = exp(1i*f*(1:N));

plot(y,z);
axis([-1 1 -1 1]);
axis square;

figure;
freq = round(f*N/(2*pi));
dev = 1e3;
dev = max(dev,1);
dev = min(dev,N);
fft1 = db(fft(y+1i*z)/N);
plot(fft1(freq-dev:freq+dev));
hold on;
fft2 = db(fft(x)/N);
plot(fft2(freq-dev:freq+dev),'r');

figure;
plot(fft1(N-freq-dev:N-freq+dev));
hold on;
plot(fft2(N-freq-dev:N-freq+dev), 'r');