figure(1) subplot(4,1,1)
plot(t,x1);title('原始信号：x(t)');xlabel('t/s'); ylabel('幅值');axis([0 2.05 -6 6]); subplot(4,1,2)
plot(t,x0);title('期望的无噪声信号：x(t)');xlabel('t/s'); ylabel('幅值');axis([0 2.05 -6 6]);
%%%%%%Wiener滤波%%%%%% od=12;
%计算互相关函数
phixs=xcorr(x1,x0); for i=1:od
rxs(i)=phixs(i+N); end
%计算自相关函数
phixx=xcorr(x1,x1); for i=1:od for j=1:od
Rxx(i,j)=phixx(i-j+N); end end
%由维纳-霍夫方程得到滤波器的最优解 h1=(inv(Rxx))*rxs'; fx=filter(h1,1,x1); subplot(4,1,3) plot(t,fx,'r');
title('Wiener滤波后信号'); % % xlabel('f/Hz'); ylabel('幅值');
axis([0 2.06 -6 6]);
%%%%%%%%%Kalman滤波去噪%%%%%%%%%%%
a1=-1.352;a2=1.338;a3=-0.662;a4=0.240; A=[-a1 -a2 -a3 -a4; 1 0 0 0; 0 1 0 0;
0 0 1 0];%状态转移矩阵 H=[1 0 0 0];%观测矩阵 Q=[1 0 0 0; 0 0 0 0; 0 0 0 0;
0 0 0 0 ];%状态噪声方差阵 R=1;%观测噪声方差阵
X(:,1)=[x1(4);x1(3);x1(2);x1(1)]; p(:,:,1)=[10 0 0 0; 0 1 0 0; 0 0 1 0;
0 0 0 1 ];%一步预测误差阵 for k=2:N
p1(:,:,k)=A*p(:,:,k-1)*A'+Q;
K(:,k)=p1(:,:,k)*H'/(H*p1(:,:,k)*H'+R); X(:,k)=A*X(:,k-1)+K(:,k)*[x1(k)-H*A*X(:,k-1)]; p(:,:,k)= p1(:,:,k)-K(:,k)*H*p1(:,:,k); end
subplot(414) plot(t,X(1,:),'r'); title('Kalman滤波后信号'); xlabel('t /s'),ylabel('幅值'); axis([0 2.06 -6 6]);
滤波前后信号图如下所示：
题5：附件中表sheet1 为某地2008年4月28日凌晨12点至2008年5月4日凌
晨12点的电力系统负荷数据，采样时间间隔为1小时，利用ARMA方法预测该地5月5日的电力系统负荷，并给出预测误差（5月5日的实际负荷数据如表
sheet2）。
解：ARMA模型算法：
ARMA模型时间序列分析法简称为时序分析法，是一种利用参数模型对有序随机振动响应数据进行处理，从而进行模态参数识别的方法。参数模型包括AR自回归模型、MA滑动平均模型和ARMA自回归滑动平均模型。
N个自由度的线性系统激励与响应之间的关系可用高阶微分方程来描述，在离散时间域内，该微分方程变成由一系列不同时刻的时间序列表示的差分方程，即ARMA时序模型方程：
2N2N?k?0akxt?k??bkft?k （5-1）
k?0 式（5-1）表示响应数据序列xt与历史值xt?k的关系，其中等式的左边称为自回归差分多项式，即AR模型，右边称为滑动平均差分多项式，即MA模型。N为自回归模型和滑动均值模型的阶次，ak、bk分别表示待识别的自回归系数和滑动均值系数，ft表示白噪声激励。当k＝0时，设a0?b0?1。 由于ARMA过程{xt}具有唯一的平稳解为
xt??hift?i （5-2）
i?0?式中：hi为脉冲响应函数。 xt的相关函数为
R??E[xtxt??]??i?0???k?0hihkE[ft?ift???k] （5-3）
ft是白噪声，故
??2k???i （5-4） E[ft?ift???k]??0other?式中：?2为白噪声方差。
将此结果代人式（5-3），即可得
R???2??i?0hihi?? （5-5）
因为线性系统的脉冲响应函数hi，是脉冲信号?，激励该系统时的输出响应，故由ARMA过程定义的表达式为
2N2N?k?0akht?k??bk?t?k?bt （5-6）
k?0 利用式（5-5）和式（5-6），可以得出：
?k?02NakRl?k??hi?ak?i?l?k??i?0k?0?2N2?i?0?hibi?l （5-7）
对于一个ARMA过程，当是大于其阶次2N时，参数bk＝0。故当l>2N时，式（5-7）恒等于零，于是有
Ri??akRl?k?0,l?2N （5-8）
k?02N或写成
?k?02NakRl?k??Ri,l?2N （5-9）
设相关函数的长度为L，并令M＝2N。对应不同的l值，由代人以上公式可得一组方程：
?a1RM?a2RM?1??a1RM?1?a2RM????aR?aR?2L?2?1L?1??aMR1?RM?1?aMR2?RM?2 （5-10）
?aMRL?M?RL或缩写为矩阵形式:
[R](L?M)?M{a}M?1?{R?}(L?M)?1 （5-11）
式（11）为推广的Yule-walker方程。一般情况下，由于L比2N大得多，采用伪逆法可求得方程组的最小二乘解，即
{a}?([R]T[R])?1([R]T{R?}) （5-12）
由此求得自回归系数ak(k?1,2,,2N)。
滑动平均模型系数bk(k?1,2,,2N)可通过以下非线性方程组来求解：
2b02?b12??bM?c0b0b1??bM?1bM?c1 （5-13）
b0bM?cM其中
ck??i?02N?j?02NaiajCk?i?j,k?0,1,2,,2N （5-14）
式中：Ck为响应序列xt的自协方差函数。
当求得自回归系数ak和滑动均值系数bk后，可以通过ARMA模型传递函数的表达式计算系统的模态参数，ARMA模型的传递函数为
2NH(z)???k?0k?02Nbkz?k （5-15）
akz?k