work.m

上载的内容为随机信号处理的作业,具体是:其中W(t)为均值为零

N=500;       %%%序列的点数
A= rand(1,N);         ....... a:=  1.73*randn(500,1)
%%产生500个随机值（噪声），其均值在0.5左右，方差在0.08左右，因此要将其调整为均值为0，方差为3的序列
average=mean(A);        %%求均值
for i=1:N
    A(i)=A(i)-average;      %%根据原均值，将噪声均值调整到0
end
variance =var(A);      %%求方差
for i=1:N
    A(i)=A(i)*sqrt( 3 / variance);     %%根据方差公式变换，将其方差调整至3
end
%%至此，噪声序列点已得到
%%进一步根据方程X(t)-1.5X(t-1)+2.5X(t-2)=A(t)+0.8A(t-1)
%%方程变形得 X(t)=A(t)+0.8A(t-1)+1.5X(t-1)-2.5X(t-2)

X=zeros(1,N);
X(1)=1;   %%两个初始点
X(2)=0.2;
%%得到完整的样点
for t=3:N
    X(t)=A(t)+0.8*A(t-1)+1.5*X(t-1)-0.5*X(t-2); 
end
subplot(211);
plot(X);


Rk=zeros(1,N);     %%计算样点的协方差
Rk0=0;
for i=1:N
    Rk0=Rk0+X(i)*X(i);  %%计算Rk0
end
Rk0=Rk0/N;
for k=1:N-1                %%
    for i=1:N-k
        Rk(k)= Rk(k)+ X(i)*X(i+k); 
    end
    Rk(k)=Rk(k)/N;
end
%%%%%%%%%计算出自相关函数
P=zeros(1,500);      %%
for i=1:500                 
    P(i)=Rk(i)/Rk0;
end
%%subplot(211);plot(P);
%%TT=xcorr(A);
%%subplot(212);plot(TT(501:999));

%%%%%%%%%%计算偏相关函数
n=N-2;B=0;C=0;
T=zeros(n);
%for e=1;N-1
  %  p(e)=p(e+1);
%end
T(1,1)=P(2);
for k=1:n;
    for i=1:k
        for j=1:k
            B=B+P(k+2-j)*T(k,j);
            C=C+P(j+1)*T(k,j);
        end
        C=1/(1-C);
        T(k+1,k+1)=(P(k+2)-B)*C;
        T(k+1,k+1-i)=T(k,k+1-i)-T(k+1,k+1)*T(k,i);
        B=0;
        C=0;
    end
end
L=ones(1,N);        %%偏相关函数
for i=1:n
    L(i)=T(i,i);
end
%%subplot(212);
%%plot(L);


%%%%%求ARMA(2,1)的参数

f=inv([1,P(1); P(1),1])*[P(1);P(2)];
f
e2=Rk0*(1-L(1)*P(1)-L(2)*P(2));

r0=e2*(1+P(1)^2);
r1=e2*(-P(1)^2);
th=-r1/e2

%%重建模型
Y=zeros(1,N);
Y(1)=1;   %%两个初始点
Y(2)=0.2;
%%得到完整的样点
for t=3:N
    Y(t)=A(t)+th*A(t-1)+f(1)*Y(t-1)+f(2)*Y(t-2);    
end
subplot(212);
plot(Y);