function eval_lfo2()

L = 80000;
fs = 48000;
fstart = 10;
fend = 1000;
df = (fend-fstart)/L;

ip = 0.0 % phase of the first output sample in radians

% Init

% Loop
f = fstart;
fchg = 0;

ylast = 0;
a = 0.5;
b = 2.0 * sin(f*pi / fs);
y0 = a*cos(2*pi*ip);
y1 = a*sin(2*pi*ip);

for n=1:L,

y0 = y0 - b*y1;
y1 = y1 + b*y0;

fchg = (mod(n, 1000) == 0);

if fchg
    b = 2.0 * sin(f*pi / fs);
end

ylast = y1;

lfo(n) = y0;
f = f + df;

close all

end;

figure;
plot(1:L, lfo); grid;

wavwrite(0.5*lfo, fs, 'lfo.wav');