西南交大现代信号处理部分答案.doc

题1：(1) 是随变化的随机信号，因此=.所以谐波信号的均值为 =由于谐波信号的均值等于零，故其方差等于二阶矩，既有 所以x(t)的方差为 谐波信号的自相关函数又所以由于x(t)的均值为0，故所以(2) y(t)是随B变化的随机信号，因此B是标准高斯随机变量，所以，所以.由于统计独立，故有而x(t)和y(t)的均值均为0，所以题2：令，由于是零均值、方差为的高斯随机过程，和是确定的过程，所以x(n)也是一高斯随机过程，其均值是时间的函数.所以x(n)的概率密度函数是=在多个未知参数的情况下，Cramer-Rao不等式变为矩阵不等式 :其中无偏估计子的协方差矩阵，而是Fisher信息矩阵J的逆矩阵，而信息矩阵J的构成元素为本题中，计算得 = = =求逆矩阵可得所以可得 所以估计子的估计方差的Cramer-Rao下界分别是和题3.信号的函数表达式为：，其中，A(t)为一随时间变化的随机过程，为经过390-410Hz带通滤波器后的高斯白噪声，为高斯白噪声，采样频率为1kHz，采样时间为2.048s。分别利用Wiener滤波和Kalman滤波进行去噪。解：维纳滤波原理：维纳(Wiener)是用来解决从噪声中提取信号的一种过滤（或滤波）方法。这种线性滤波问题，可以看做是一种估计问题或一种线性估计问题。一个线性系统，如果它的单位样本响应为，当输入一个随机信号，且：其中：表示信号，表示噪声，则输出为： 我们希望通过线性系统后得到的尽量接近于，因此称为的估计值，用表示，即：则维纳滤波器的输入输出关系可用下面图1表示。图3-1 维纳滤波器的输入输出关系实际上,上式所示的卷积形式可以理解为从当前和过去的观察值，来估计信号的当前值。因此，用进行过滤问题实际上是一种统计估计问题。一般地，从当前的和过去的观察值，估计当前的信号值成为过滤或滤波；从过去的观察值，估计当前的或者将来的信号值称为外推或预测；从过去的观察值，估计过去的信号值称为平滑或内插。因此维纳滤波器又常常被称为最佳线性过滤与预测或线性最优估计。这里所谓的最佳与最优是以最小均方误差为准。如果我们分别以与表示信号的真实值与估计值，而用表示他们之间的误差，即：显然可能是正值，也可能是负值，并且它是一个随机变量。因此，用它的均方误差来表达误差是合理的，所谓均方误差最小即它的平方的统计期望最小。卡尔曼滤波算法原理：卡尔曼滤波是基于状态空间方法的一套递推滤波算法，在状态空间方法中，引入了状态变量的概念。卡尔曼滤波的模型包括状态空间模型和观测模型。状态模型是反映状态变化规律的模型，通过状态方程来描写相邻时刻的状态转移变化规律；观测模型反映了实际观测量与状态变量之间的关系。Kalman滤波问题就是联合观测信息及状态转移规律来得到系统状态的最优估计。假设动态系统的状态空间模型为其中，X(t)为系统在时刻t的状态；Y(t)为对状态的观测值；W(t)为系统噪声，方差阵为Q；V(t)为观测噪声，方差阵为R；为状态转移矩阵；H为观测矩阵；为系统噪声驱动矩阵。卡尔曼滤波的计算流程为：计算状态估计值：计算状态一步预测：计算新息：计算卡尔曼滤波增益：计算一步预测均方误差：计算一步预测估计均方误差：下面给出卡尔曼滤波的系统模型框图：根据上述要求，编写程序如下：fs = 1000; %采样频率t = 0:1/fs:2.047;N = length(t); %采样点数randn(state,0);W1=randn(1,length(t); %高斯白噪声% 一维数字滤波 %f2=390;f3=410; %带通滤波带宽设置wc1=2*f2/fs;wc2=2*f3/fs;f1=0 wc1-0.05 wc1 wc2 wc2+0.05 1; %归一化频率f1A=0 0 1 1 0 0; %设置带通1或带阻0weigh=1 1 1 ; %权重b=remez(50,f1,A,weigh); %传函分子W2=filter(b,1,W1); %一维数字滤波B=normrnd(0,1,1,N); %随机过程A（t）x0= sin(2*pi*t*100) +1.5*sin(2*pi*t*300)+B.*sin(2*pi*t*200);x1 = sin(2*pi*t*100) +1.5*sin(2*pi*t*300)+B.*sin(2*pi*t*200)+W1+W2; % 原始信号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:odrxs(i)=phixs(i+N);end%计算自相关函数phixx=xcorr(x1,x1);for i=1:odfor j=1:odRxx(i,j)=phixx(i-j+N);endend%由维纳-霍夫方程得到滤波器的最优解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:Np1(:,:,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);endsubplot(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时序模型方程： （5-1） 式（5-1）表示响应数据序列与历史值的关系，其中等式的左边称为自回归差分多项式，即AR模型，右边称为滑动平均差分多项式，即MA模型。N为自回归模型和滑动均值模型的阶次，、分别表示待识别的自回归系数和滑动均值系数，表示白噪声激励。当k0时，设。 由于ARMA过程具有唯一的平稳解为 （5-2）式中：为脉冲响应函数。 的相关函数为 （5-3） 是白噪声，故 （5-4）式中：为白噪声方差。 将此结果代人式（5-3），即可得 （5-5） 因为线性系统的脉冲响应函数，是脉冲信号，激励该系统时的输出响应，故由ARMA过程定义的表达式为 （5-6） 利用式（5-5）和式（5-6），可以得出： （5-7） 对于一个ARMA过程，当是大于其阶次2N时，参数0。故当l2N时，式（5-7）恒等于零，于是有 （5-8）或写成 （5-9） 设相关函数的长度为L，并令M2N。对应不同的l值，由代人以上公式可得一组方程： （5-10）或缩写为矩阵形式: （5-11） 式（11）为推广的Yule-walker方程。一般情况下，由于L比2N大得多，采用伪逆法可求得方程组的最小二乘解，即 （5-12）由此求得自回归系数。滑动平均模型系数可通过以下非线性方程组来求解： （5-13）其中 （5-14）式中：为响应序列的自协方差函数。 当求得自回归系数和滑动均值系数后，可以通过ARMA模型传递函数的表达式计算系统的模态参数，ARMA模型的传递函数为 （5-15）由以上分析，利用Matlab仿真，ARMA模型图形如下：图1 电力系统负荷数据图2 平稳化处理后电力系统负荷数据图3 电力系统负荷数据自相关和偏相关函数曲线图4 使用ARMA模型对电力系统负荷进行预测曲线图4 预测误差曲线由仿真结果分析得，使用ARMA模型对电力系统负荷进行预测，其预测结果与实际结果相比较而言，预测误差较小，在可接受范围内，所以预测结果可以采用。Matlab程序如下：clc;clear;x1=xlsread(C:UsersdellDesktop现代信号负荷数据.xls,Sheet1);x1=x1(:,2);x2=xlsread(C:UsersdellDesktop现代信号负荷数据.xls,Sheet2);x2=x2(:,2);x=x1;x2; % 电力系统负荷数据，前168个点为历史数据，后24个点为待预测数据N1=length(x1);N=length(x);t=1:N;figure;plot(x);title(电力系统负荷数据);xlabel(时间);ylabel(电力系统负荷);axis tight; %对原始数据进行去趋势处理，即零均值化、平稳化处理z1=diff(x); %差分Y=z1-mean(z1); %零均值H,PValue,TestStat,CriticalValue = adftest(Y); %检验是否为时间平稳序列 figure;plot(Y,r)title(平稳化处理后电力系统负荷数据);xlabel(时间),ylabel(功率); %计算自相关函数、偏相关函数figuresubplot(2,1,1)autocorr(Y); %计算置信度为95%的acf，并画出其自相关函数曲线；subplot(2,1,2) parcorr(Y); %计算置信度为95%的pacf，并画出其偏自相关函数曲线； %自相关系数和偏相关系数均拖尾性，所以选择arma模型对电力系统负荷数据进行预测 %估计arma模型参数，以AIC标准来定阶z=iddata(Y);test = ; for p = 1:10 %自回归对应PACF,给定滞后长度上限p和q for q = 1:10 %移动平均对应ACF m = armax(z,p q); AIC = aic(m); %armax(p,q),选择对应FPE最小，AIC值最小的模型 test = test;p q AIC; endendfor k = 1:size(test,1) if test(k,3) = min(test(:,3) %选择AIC值最小的模型 p_test = test(k,1); q_test = test(k,2); break; endend %建立arma模型检验m = armax(z(1:168),p_test q_test); %armax(p,q),p_test q_test对应AIC值最小figuree = resid(m,z(1:168); %拟合做残差分析plot(e); figuresubplot(2,1,1)autocorr(e.OutputData) %置信水平0.95，检验残差的自相关和偏相关函数subplot(2,1,2)parcorr(e.OutputData) Pr,DWr = dwtest(e.OutputData,z.OutputData); %检验线性回归残差是否相互独立if Pr0.05 disp(can not use this model);else disp(can use this model);end %预测5月5日电力系统负荷数据p=predict(m,z,1);t=1:1:191;po = p.OutputData;figure;plot(t,z,r,t,po,b); %显示预测值（进行反差分、加平均值）X(1)=x(1);for i=2:192X(i)=x(i-1)+po(i-1); endX=X+mean(po); t=1:1:192;figure;plot(t,x,r,t,X,b);title(使用ARMA模型对电力系统负荷数据进行预测);xlabel(时间);ylabel(电力系统负荷);axis tight;l