answer2.asv

利用burg和Levinson-Durbin算法

%根据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,256).^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(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算法得到的功率谱曲线');