gonglvpu.m

计算太阳黑子的活动周期

clear;
close all;
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];
Fs=1000;
%周期图法求功率谱
fftn=1024;
X=fft(x,fftn);
P=abs(X).^2/length(x);
k=0:(fftn/2-1);
Pu=P(k+1);
f=k*Fs/fftn;                   %模拟频率f=w*Fs/(2*pi)
figure;
plot(f,Pu)

%Levinson算法
p=20;
[a,e2]=Levinson(x,p); 
s=0;
w=0:0.00001:pi;
for k=1:p
    s=s+(a(p,k).*exp(-j*k*w)); 
end
P=e2(p)./((abs(1+s)).^2);
f=w*Fs/(pi*2);
hold on;
plot(f,abs(P),'r--');


%burg算法
[a,e2]=burg(x,p);
s=0;
w=0:0.00001:pi;
for k=1:p
    s=s+(a(p,k).*exp(-j*k*w));
end
P=e2(p+1)./((abs(1+s).^2));
f=w*Fs/(pi*2);
hold on;
plot(f,abs(P),'k--');
axis([30 200 0 25000]);;grid on;
legend('周期图法','Levinson算法','burg算法');