% parameters for signals
o1   = 59*pi/200;  phi1 = pi/4;  a1   = 1;
o2   = 139*pi/200; phi2 = pi/3;  a2   = .5;
o3   = 141*pi/200; phi3 = -pi/3; a3   = .5;
% parameters for noise
order = 1; sigma = 1.5;
a     = [1.0000   -.85];
% generate data
n     = 150; t=0:n-1;
T     = tf([1 zeros(1,order)],a,-1);
randn('seed',101); noise = sigma*randn(n,1);
sig   = sin(o1*t+phi1)+.5*sin(o2*t+phi2)+.5*sin(o3*t+phi3);
y     = lsim(T,noise)+sig';
% estimate
freqNum = 3; N = 200;
[omega,mag,phi,a,sigma_est]=L1AR(y,freqNum,order,[0 pi],N);
% plot results
w     = linspace(0,pi,500);
ptrue = abs(freqresp(T,w)).^2;
T_est = tf([1 zeros(1,order)],a,-1);
pest  = abs(freqresp(T_est,w)).^2*sigma_est^2;
figure(1); periodogram(y);
figure(2);
    subplot(2,1,1);
        plot(w,vec(ptrue),'r','LineWidth',1.2);hold on;
        plot(o1,2*pi*a1,'r^','MarkerSize',12);
        legend('true spectrum of colored noise','true spectral lines');
        set(gca,'YScale','log','XLim',[0 pi]);
        arrow([o1 o2 o3],(2*pi*[a1 a2 a3]),12);
    subplot(2,1,2);
        plot(w,vec(pest),'b','LineWidth',1.2);hold on;
        plot(omega(1),2*pi*mag(1),'b^','MarkerSize',12);
        legend('estimaged spectrum of colored noise','estimated spectral lines');
        set(gca,'YScale','log','XLim',[0 pi]);
        arrowb(omega,(2*pi*mag),12);