ganhaohao
2009-03-10, 10:35
下面的数据是从1770年至1869年这100年间每年12个月所记录到的太阳黑子出现次数的平均值。
a. 用周期图法做该数据的功率谱，从该曲线您是否可以看出太阳黑子活动的周期？
b. 对该数据建立一个三阶的AR模型，分别用Levinson-Durbin法和Burg法求解AR模型的参数，并得到功率谱曲线，从该曲线来发现太阳黑子的活动周期。改变AR模型的阶次，看不同阶次下所估计出的周期的差别。
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];
下面是matlab程序：
clear;
clc;
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);
R=ZXG(x);
%han=zeros(1,1024);
%定义一个窗长为N的汉明窗
for i=1:N
han(i)=100*(0.54-0.46*cos(2*pi*i/N));
end
%用周期图法求数据的功率谱估计P
X=abs(fft(x,1024));
for i=1:1024
P(i)=X(i)*X(i)/N;
end
plot(P);
%改进之后的周期图法,u为汉明窗
for i=1:99
u(i)=R(i)*han(i);
end
UU=abs(fft(u,1024));
figure(2);
plot(UU);
axis([0 100 0 1000000]);
grid on
%三阶AR模型参数得其功率谱
a_11=-R(2)/R(1);
D_1=(1-abs(a_11)*abs(a_11))*R(1);
a_22=-(R(3)+a_11*R(2))/D_1;
a_21=a_11+a_22*a_11;
D_2=(1-abs(a_22)*abs(a_22))*D_1;
a_33=-(R(4)+a_21*R(3)+a_22*R(2))/D_2;
a_32=a_22+a_33*a_21;
a_31=a_21+a_33*a_22;
D_3=(1-abs(a_33)*abs(a_33))*D_2;
A=[a_31 a_32 a_33];
D=abs(fft(A,256));
for i=1:256
PP(i)=D_3/((1+D(i))*(1+D(i)));
end
figure(3);
plot(PP);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%因为AR谱估计高阶算法与三阶算法类似， %
%%%而这种方法与线性预测谱估计等效，故 %
%%%可以直接直接调用求线性预测的系统函数LPC(X,N)
L=LPC(x,70);
LL=abs(fft(L,256));
for i=1:256
PPP(i)=D_3/((1+LL(i))*(1+LL(i)));
end
figure(4);
plot(PPP);
%axis([0 128 0 5]);
此程序主要用于输入数据的功率谱估计，该程序由三部分组成：burg_Kp.m，Mainfun.m，ZXG.m.其中：Mainfun.m为主函数，burg_Kp.m为Burg算法函数，ZXG.m为求输入信号自相关的函数(估计)。
问题在于，如何在功率谱密度图中添加一条信度值为0.9的置信曲线？
a. 用周期图法做该数据的功率谱，从该曲线您是否可以看出太阳黑子活动的周期？
b. 对该数据建立一个三阶的AR模型，分别用Levinson-Durbin法和Burg法求解AR模型的参数，并得到功率谱曲线，从该曲线来发现太阳黑子的活动周期。改变AR模型的阶次，看不同阶次下所估计出的周期的差别。
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];
下面是matlab程序：
clear;
clc;
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);
R=ZXG(x);
%han=zeros(1,1024);
%定义一个窗长为N的汉明窗
for i=1:N
han(i)=100*(0.54-0.46*cos(2*pi*i/N));
end
%用周期图法求数据的功率谱估计P
X=abs(fft(x,1024));
for i=1:1024
P(i)=X(i)*X(i)/N;
end
plot(P);
%改进之后的周期图法,u为汉明窗
for i=1:99
u(i)=R(i)*han(i);
end
UU=abs(fft(u,1024));
figure(2);
plot(UU);
axis([0 100 0 1000000]);
grid on
%三阶AR模型参数得其功率谱
a_11=-R(2)/R(1);
D_1=(1-abs(a_11)*abs(a_11))*R(1);
a_22=-(R(3)+a_11*R(2))/D_1;
a_21=a_11+a_22*a_11;
D_2=(1-abs(a_22)*abs(a_22))*D_1;
a_33=-(R(4)+a_21*R(3)+a_22*R(2))/D_2;
a_32=a_22+a_33*a_21;
a_31=a_21+a_33*a_22;
D_3=(1-abs(a_33)*abs(a_33))*D_2;
A=[a_31 a_32 a_33];
D=abs(fft(A,256));
for i=1:256
PP(i)=D_3/((1+D(i))*(1+D(i)));
end
figure(3);
plot(PP);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%因为AR谱估计高阶算法与三阶算法类似， %
%%%而这种方法与线性预测谱估计等效，故 %
%%%可以直接直接调用求线性预测的系统函数LPC(X,N)
L=LPC(x,70);
LL=abs(fft(L,256));
for i=1:256
PPP(i)=D_3/((1+LL(i))*(1+LL(i)));
end
figure(4);
plot(PPP);
%axis([0 128 0 5]);
此程序主要用于输入数据的功率谱估计，该程序由三部分组成：burg_Kp.m，Mainfun.m，ZXG.m.其中：Mainfun.m为主函数，burg_Kp.m为Burg算法函数，ZXG.m为求输入信号自相关的函数(估计)。
问题在于，如何在功率谱密度图中添加一条信度值为0.9的置信曲线？