load yearssn

N=256;
t=yearssn(1:N,1);
x=yearssn(1:N,2);
plot(t,x);

figure
subplot(2,2,1);
plot(t,x);
xlabel('Year');
ylabel('Sunspot number');
hold on

subplot(2,2,2);
plot(hanning(N));
title('Hanning window');
xlabel('Time');

[p,ff]=psd(x,[],Fs);
subplot(2,2,3);
plot(ff,p);
xlabel('Frequency');
ylabel('Power');
title('Modified periodogram');

subplot(2,2,4);
plot(ff,log(p));
xlabel('Frequency');
ylabel('Log Power');
title('Log modified periodogram');

figure
subplot(2,2,1);
plot(ff,log(p));
grid on
xlabel('Frequency');
ylabel('Log Power');
title('Log modified periodogram');

[pw,fw]=pwelch(x,[],[],[],Fs);
subplot(2,2,2);
plot(fw,log(pw));
grid on
xlabel('Frequency');
ylabel('Log Power');
title('Log Welch');

nw=4;
w=nw/N;
[pm1,fm1]=pmtm(x,nw,[],Fs);
subplot(2,2,3);
plot(fm1,log(pm1));
grid on
xlabel('Frequency');
ylabel('Log Power');
title(sprintf('Log Multitaper - SpecRes=%1.2f',w));

nw=2;
w=nw/N;
[pm2,fm2]=pmtm(x,nw,[],Fs);
subplot(2,2,4);
plot(fm1,log(pm2));
grid on
xlabel('Frequency');
ylabel('Log Power');
title(sprintf('Log Multitaper - SpecRes=%1.2f',w));

figure
nw=4;
[e,v]=dpss(N,nw);
for i=1:4,
    subplot(2,2,i);
    plot(e(:,i));
    xlabel('Time');
    title(sprintf('Slepian sequence %d',i));
end