%根據(jù)100年間每年12個月記錄到的太陽黑子出現(xiàn)次數(shù)的平均值求太陽黑子的活動周期
clear all;
%100年間記錄的太陽黑子出現(xiàn)平均數(shù)
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);

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

%AR模型  
%Levinson-Durbin算法
%3階
M=3;%M為階數(shù)
[Pxx]=LD(M,N,x);
figure;
subplot(2,1,1);
plot(Pxx);
%axis([0 150 0 20000]);
title('3階Levinson-Durbin算法得到的功率譜曲線');
%20階
M=10;
[Pxx]=LD(M,N,x);
subplot(2,1,2);
plot(Pxx);
axis([8 13 0 15000]);
title('20階Levinson-Durbin算法得到的功率譜曲線');
%可以明顯的看到n大概在23時有峰值，T=256/n=11，11年是一個周期

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