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

题1：(1) 错误！未找到引用源。是随错误！未找到引用源。变化的随机信号，因此错误！未找到引用源。=错误！未找到引用源。. 所以谐波信号x(t)的均值为 错误！未找到引用源。 =错误！未找到引用源。

由于谐波信号x(t)的均值等于零，故其方差等于二阶矩，既有错误！未找到引用源。

错误！未找到引用源。

所以x(t)的方差为 错误！未找到引用源。 谐波信号x(t)的自相关函数错误！未找到引用源。 又错误！未找到引用源。 所以错误！未找到引用源。

由于x(t)的均值为0，故所以错误！未找到引用源。 (2) y(t)是随B变化的随机信号，因此错误！未找到引用源。

B是标准高斯随机变量，所以错误！未找到引用源。，所以错误！未找到引用源。. 由于错误！未找到引用源。统计独立，故有 错误！未找到引用源。

而x(t)和y(t)的均值均为0，所以错误！未找到引用源。

题2：令错误！未找到引用源。，由于错误！未找到引用源。是零均值、方差为错误！未找到引用源。的高斯随机过程，错误！未找到引用源。和错误！未找到引用源。是确定的过程，所以x(n)也是一高斯随机过程，其均值错误！未找到引用源。是时间的函数.所以x(n)的概率密度函数是

f(x;?)??n?1N12??2exp{?12?212?N2[x(n)?A?Bn]2}

1=

(2??2)N/2exp{??n?1[x(n)?A?Bn]2}

在多个未知参数的情况下，Cramer-Rao不等式变为矩阵不等式 :

??J?1(?)

其中错误！未找到引用源。无偏估计子错误！未找到引用源。的协方差矩阵，而错误！未找到引用源。是Fisher信息矩阵J的逆矩阵，而信息矩阵错误！未找到引用源。的构成元素为

错误！未找到引用源。

本题中，计算得

错误！未找到引用源。 错误！未找到引用源。=错误！未找到引用源。 错误！未找到引用源。 错误！未找到引用源。

错误！未找到引用源。 错误！未找到引用源。=错误！未找到引用源。 =错误！未找到引用源。

求逆矩阵可得

所以可得

2(2N?1)?2 varA(错误！未找到引用源。 )?N(N?1)?错误！未找到引用源。

所以估计子错误！未找到引用源。的估计方差的Cramer-Rao下界分别是错误！

12?22未找到引用源。和N(N?1)错误！未找到引用源。

题3.信号的函数表达式为： x(t)?sin(2?100t)?1.5?sin(2?300t)?A(t)sin(2?200t)?dn(t)?n(t)，

其中，A(t)为一随时间变化的随机过程，dn(t)为经过390-410Hz带通滤波器后的高斯白噪声，n(t)为高斯白噪声，采样频率为1kHz，采样时间为2.048s。分别利用Wiener滤波和Kalman滤波进行去噪。

解：维纳滤波原理：

维纳(Wiener)是用来解决从噪声中提取信号的一种过滤（或滤波）方法。这种线性滤波问题，可以看做是一种估计问题或一种线性估计问题。一个线性系统，如果它的单位样本响应为h(n)，当输入一个随机信号x(n)，且：

x(n)?s(n)?v(n)

其中：x(n)表示信号，v(n)表示噪声，则输出y(n)为：

y(n)??h(n)x(n?m)

m我们希望x(n)通过线性系统h(n)后得到的y(n)尽量接近于s(n)，因此称y(n)为s(n)的估计值，用s(n)表示，即：

y(n)?s(n)

则维纳滤波器的输入—输出关系可用下面图1表示。

图3-1 维纳滤波器的输入—输出关系

实际上,上式所示的卷积形式可以理解为从当前和过去的观察值

x(n),x(n?1),,x(n?m)，来估计信号的当前值s(n)。因此，用h(n)进行过滤问题

实际上是一种统计估计问题。

一般地，从当前的和过去的观察值x(n),x(n?1),,x(n?m)，估计当前的信号值y(n)?s(n)成为过滤或滤波；从过去的观察值，估计当前的或者将来的信号值

y(n)?s(n?N)(N?0)称为外推或预测；从过去的观察值，估计过去的信号值y(n)?s(n?N)(N?0)称为平滑或内插。因此维纳滤波器又常常被称为最佳线性过

滤与预测或线性最优估计。这里所谓的最佳与最优是以最小均方误差为准。

如果我们分别以s(n)与s(n)表示信号的真实值与估计值，而用e(n)表示他们之间的误差，即：

e(n)?s(n)?s(n)

显然e(n)可能是正值，也可能是负值，并且它是一个随机变量。因此，用它的均方误差来表达误差是合理的，所谓均方误差最小即它的平方的统计期望最小。

卡尔曼滤波算法原理：

卡尔曼滤波是基于状态空间方法的一套递推滤波算法，在状态空间方法中，引入了状态变量的概念。卡尔曼滤波的模型包括状态空间模型和观测模型。状态模型是反映状态变化规律的模型，通过状态方程来描写相邻时刻的状态转移变化规律；观测模型反映了实际观测量与状态变量之间的关系。Kalman滤波问题就是联合观测信息及状态转移规律来得到系统状态的最优估计。

假设动态系统的状态空间模型为

X(t?1)??X(t)??W(t) Y(t)?HX(t)?V(t)

其中，X(t)为系统在时刻t的状态；Y(t)为对状态的观测值；W(t)为系统噪声，方差阵为Q；V(t)为观测噪声，方差阵为R；?为状态转移矩阵；H为观测矩阵；

?为系统噪声驱动矩阵。

卡尔曼滤波的计算流程为： 计算状态估计值：

X?(t?1|t?1)?X?(t?1|t)?K(t?1)?(t?1) 计算状态一步预测：

X?(t?1|t)??X?(t|t) 计算新息：

?(t?1)?Y(t?1)?HX?(t?1|t) 计算卡尔曼滤波增益：

K(t?1)?P(t?1|t)HT[HP(t?1|t)HT?R]?1

计算一步预测均方误差：

P(t?1|t)??P(t|t)?T??Q?T

计算一步预测估计均方误差：

P(t?1|t?1)?[In?K(t?1)H]P(t?1|t)

下面给出卡尔曼滤波的系统模型框图：

根据上述要求，编写程序如下：

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]; %归一化频率f1 A=[0 0 1 1 0 0]; %设置带通1或带阻0 weigh=[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: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;%观测噪声方差阵

图4 使用ARMA模型对电力系统负荷进行预测曲线

图4 预测误差曲线

由仿真结果分析得，使用ARMA模型对电力系统负荷进行预测，其预测结果与实际结果相比较而言，预测误差较小，在可接受范围内，所以预测结果可以采用。

Matlab程序如下： clc; clear;

x1=xlsread('C:\\Users\\dell\\Desktop\\现代信号\\负荷数据.xls','Sheet1'); x1=x1(:,2);

x2=xlsread('C:\\Users\\dell\\Desktop\\现代信号\\负荷数据.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('功率');

%计算自相关函数、偏相关函数 figure

subplot(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]; end end

for k = 1:size(test,1)

if test(k,3) == min(test(:,3)) %选择AIC值最小的模型 p_test = test(k,1); q_test = test(k,2); break; end end

%建立arma模型检验

m = armax(z(1:168),[p_test q_test]); %armax(p,q),[p_test q_test]对应AIC值最小 figure

e = resid(m,z(1:168)); %拟合做残差分析 plot(e);

figure

subplot(2,1,1)

autocorr(e.OutputData) %置信水平0.95，检验残差的自相关和偏相关函数 subplot(2,1,2)

parcorr(e.OutputData)

[Pr,DWr] = dwtest(e.OutputData,z.OutputData); %检验线性回归残差是否相互独立

if Pr<0.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:192

X(i)=x(i-1)+po(i-1);

end

X=X+mean(po);

t=[1:1:192]; figure;

plot(t,x,'r',t,X,'b');title('使用ARMA模型对电力系统负荷数据进行预测');xlabel('时间');ylabel('电力系统负荷');axis tight;legend('负荷真实值','ARMA预测值'); figure;

subplot(2,1,1); t=[168:1:192];

plot(t,x(168:192),'r',t,X(168:192),'b');title('使用ARMA模型对电力系统负荷数据进行预测');xlabel('时间');ylabel('电力系统负荷');axis tight;legend('负荷真实值','ARMA预测值'); subplot(2,1,2); error=x-X';

plot(t,error(168:192));

title('预测值与真实之误差');xlabel('时间');ylabel('5月5日预测值与真实值误差');axis tight;

题6

本文来源：https://www.bwwdw.com/article/e6d6.html