%根據100年間每年12個月記錄到的太陽黑子出現次數的平均值求太陽黑子的活動周期
clear all;
%100年間記錄的太陽黑子出現平均數
x=[101 82 66 35 31 7 20 92 154 125 85 68 38 23 10 24 83 132 131 118 90 67 60 47 41 21 16 6 4 7 14 34 45 43 48 42 28 10 8 2 0 1 5 12 14 35 46 41 30 24 16 7 4 2 8 17 36 50 62 67 71 48 28 8 13 57 122 138 103 86 63 37 24 11 15 40 62 98 124 96 66 64 54 39 21 7 4 23 55 94 96 77 59 44 47 30 16 7 37 74];
N=length(x);

%畫原數據
figure;
subplot(2,1,1);
plot(x);
title('記錄數據');
%可以根據原數據看出大概11年是一個周期
%畫功率譜
p=abs(fft(x,100).^2)/N;%計算數據的功率譜
subplot(2,1,2);
plot(p);
axis([0 50 0 100000]);
title('功率譜');
%可以看到n=23時有峰值，但是旁瓣較多

%AR模型  
%Levinson-Durbin算法
%3階
M=3;%M為階數
[Pxx]=LD(M,N,x);
figure;
subplot(3,1,1);
plot(Pxx);
title('3階Levinson-Durbin算法得到的功率譜曲線');
%10階
M=20;
[Pxx]=LD(M,N,x);
subplot(3,1,2);
plot(Pxx);
axis([0 50 0 100000]);
title('20階Levinson-Durbin算法得到的功率譜曲線');
subplot(3,1,3);
plot(Pxx);
axis([8 13 0 20000]);
title('在8~13區間觀察功率譜');
%可以明顯的看到n在11時有峰值

%burg算法
%3階
M=3;
[Pxx]=BG(M,N,x);
%Pxx=pburg(x,M,256);
figure;
subplot(3,1,1);
plot(Pxx);
title('3階burg算法得到的功率譜曲線');
%10階
M=10;
[Pxx]=BG(M,N,x);
subplot(3,1,2);
% Px=pburg(x,M,256);
plot(Pxx);
axis([0 50 0 100000]);
title('20階burg算法得到的功率譜曲線');
subplot(3,1,3);
plot(Pxx);
axis([8 13 0 20000]);
title('在8~13區間觀察功率譜');