1、考虑两个谐波信号x(t)和y(t)，其中x(t)?Acos(wct??),y(t)?Bcos(wct)式中A和wc为正的常数,?为均匀分布的随机变量，其概率密度函数为
?1?,0???2?, f(?)??2??0,其他?而B是一个具有零均值和单位方差的标准高斯随机变量，即其分布函数为
1fB(b)?exp(?b2/2),???b??
2?（1）求x(t)的均值ux(t)、方差?2x(t)、自相关函数Rx(?)和自协方差函数cx(?)。
（2）若?与B为相互统计独立的随机变量，求x(t)和x(t)的互相关函数Rxy(?)与互协方差函数cxy(?)。
解： （1）
x(t)的均值ux(t)为：
1ux?E(x(t))?E(Acos(wct??))?2??2?01Acos(wct??)d??Asin(wct??) ?0
2?02?方差?2x(t)为：
A2A2A2A2 ??E?(x(t))?E(Acos(wct??))?E(1?cos(2wct?2?))??E(cos(2wct?2?))?2222自相关函数Rx(?)为：
2x222Rx(?)?E(Acos(wct??)Acos(wct+wc???))?A2E(cos(wct??)cos(wct+wc???)) A2A2A2A2?E(cos(2wct+2wc??2?)?cos(wc?))?cos(wc?)?E(cos(2wct+2wc??2?))?cos(wc?)2222自协方差函数cx(?)为：
A2cx(?)?Rx(?)?cos(wc?)
2（2）y(t)的均值为：
uy(t)?EB(y(t))?EB(Bcos(wct))?E(B)cos(wct)?0,所以E(B)=0
由互相关函数的定义可知：
Rxy(?)?E(Acos(wct??)Bcos(wct?wc?))
由题意知道?与B为相互统计独立的随机变量，所以有
Rxy(?)?E(Acos(wct??)E(Bcos(wct?wc?))?AE(cos(wct??)E(B)cos(wct?wc?))?A?0?0?cos(wct?wc?)?0
互协方差函数cxy(?)
cxy(?)?Rxy(?)?0
2.接收信号由下式给出：yi?Acos(?ci??)??i,i?1,2,...,N，式中?i~N(0,1)即?i是零均
?c为载波角频率，值和单位方差的高斯噪声，而?是未知的相位。假设?1,?2,...?N相互独立，求未知相位的最大似然估计?ML。
解：由于?1,?2,...?N相互独立，所以y1,..yN也相互独立并且服从高斯分布，可以得到y1,..yN与?的联合概率密度函数分布
f(y1,..yN|?)?1(2?)N?^e?[yi?Acos(?ci??)]22i?1N
由此，可以得到似然函数
N1NL??ln(2?)??[yi?Acos(?ci??)]
22i?12该似然函数对?求偏导，并令该偏导函数为零，即可得到如下公式：
?LN??[yi?Acos(?ci??)]sin(?ci??)?0 ??i?1因此，最大似然估计?ML为上述函数的零点值。 则
^?Acos(?ci??ML)sin(?ci??ML)??yisin(?ci??ML)
i?1i?1N^^N^该函数为非线性方程，不容易求解，若忽略双倍频率2?c，则可简化到如下式子：
?yisin(?ci??)?0
i?1N^根据三角公式分解得到如下式子：
sin?ML?yicos(?ci)??cos?ML?yisin(?ci)
i?1i?1^N^N由此，可以得到如下公式
tan?ML??^i?1N?ysin?iicicN
?ycos(?i)i?1所以相位的最大似然估计如下：
?ML??arctan(^?ysin(?i)icN?ycos(?i)ici?1i?1N)
3.离散时间的二阶AR过程由差分方程x(n)?a1x(n?1)?a2x(n?2)?w(n)描述，式中w(n)2是一零均值、方差为?w的白噪声。证明x(n)的功率谱为
2?w21?a12?a2?2a1(1?a2)cos(2?f)?2a2cos(4?f)Px(f)?
证明：
由AR过程的功率谱公式知
Px(f)?2???j4?f21?a1e?j2?f?a2e
其中
1?a1e?j2?f?a2ej4?f2?(1?a1e?j2?f?a2ej4?f)(1?a1e?j2?f?a2ej4?f)2?1?a12?a2?a2(e?j4?f?e4?f)?a1a2ej2?f?a1a2e?j2?f?a1(e?j2?f?ej2?f) 2?1?a12?a2?2a1(1?a2)cos(2?f)?2a2cos(4?f)将其带入第一个公式可得：
Px(f)?2??21?a12?a2?2a1(1?a2)cos(2?f)?2a2cos(4?f)
x?t??sin(2?100t)?1.5sin(2?300t)?A?t?sin(2?200t)?dn?t??n?t?，4、信号的函数表达式为：
其中，A?t?为一随时间变化的随机过程，dn?t?为经过390-410Hz带通滤波器后的高斯白噪声，n?t?为高斯白噪声，采样频率为1kHz，采样时间为2.048s。分别利用周期图谱、ARMA、Burg最大熵方法估计信号功率谱，其中ARMA方法需要讨论定阶的问题。
解：由题意知采样点数一共为：1000×2.048=2048个数据点。A?t?为一随时间变化的随机过程，由于随机过程有很多类型，如维纳过程、正态随机过程，本文采用了均值为0，方差为1的正态随机过程来作为演示，来代替A?t?，高斯白噪声采用强度为2的高斯白噪声代替n?t?，其带通滤波后为dn?t?。其中滤波器采用的是契比雪夫数字滤波器。 可得到x(t)如下图所示：
1、周期图法
matlab中的周期图功率谱法原理是通过计算采样信号的FFT，获得离散点的幅度，再根据幅度与功率之间的关系，转换为离散点的功率，再通过坐标变换将离散点的功率图转换为连续功率谱密度。
Step1:计算采样信号x(n)的DFT，使用FFT方法来计算。如果此处将复频率处的幅度对称到物理实际频率，得到的就是单边谱，否则就是双边谱
Step2:根据正余弦信号功率与幅度的关系以及直流功率与幅度的关系，将幅度转换为离散功率谱。
Step3:对横纵坐标进行转换，横坐标乘以频率分辨率转换为实际连续物理频率，纵坐标除以频率分辨率转换为功率谱密度。
调用MATLAB中自带的matlab中[Pxx,f]=periodogram(x,window,nfft,fs)函数可得计算结果如下：
2、ARMA方法
参数模型估计的思想是:
?假定研究的过程X(n)是一个输入序列u(n)激励一个线性系统H(z)的输出。 ?有已知的X(n)，或其自相关函数来估计H(z)的参数。 ?由H(z)的参数来估计X(n)的功率谱。
不论X(n)是确定性信号还是随机信号，u(n)与X(n)之间总有如下输入输出关系：
x(n)???akx(n?k)??bku(n?k)
k?1k?0pqx(n)??h(k)u(n?k)
k?0?对以上两个式子两边分别取Z变换，并假定b0=1，可得
H(z)?pq?kB(z) A(z)?其中A(z)?1??akz,B(z)?1??bkz?k,H(z)??h(k)z?k。
k?1k?1k?0为了保证H(z)是稳定的最小相位系统，A(z)和B(z)的零点都应该在单位圆内。假定u(n)是一个方差为?2的白噪声序列，由随机信号通过线性系统的理论可知，输出序列X(n)的功率谱为：
Px(ejw)?2?2B*(ejw)B(ejw)A*(ejw)A(ejw)??2B(ejw)A(e)jw2
ARMA阶数确定：
本题目采用AIC准则确定ARMA的阶数。分别计算p、q从1到20阶数的计算出AIC（p,q）,如下图所示，当横坐标大概为230左右时，AIC(p,q)取得最小，将此时的p,q作为带入到模型即可。
ARMA法谱估计结果：
3、Burg最大熵法
Burg算法的具体实现步骤：
步骤1 计算预测误差功率的初始值和前、后向预测误差的初始值，并令m = 1。
1P0?N2
f0(n)?g0(n)?x(n)
n?1?x(n)N步骤2 求反射系数
?Km?Nn?m?1?fNm?1?(n)gm?1(n?1)122[fm?1(n)?gm?1(n?1)]?2n?m?1
步骤3 计算前向预测滤波器系数
?am(i)?am?1(i)?Kmam,...,m?1 ?1(m?i),i?1am(m)?Km
步骤4 计算预测误差功率
Pm?(1?Km)Pm?12
步骤5计算滤波器输出
fm(n)?fm?1(n)?Kmgm?1(n?1)
?gm(n)?Kmfm?1(n)?gm?1(n?1)
步骤6 令m←m+1，并重复步骤2至步骤5，直到预测误差功率Pm不再明显减小。 最后，再利用Levinson递推关系式估计AR参数，继而得到功率谱估计。 Burg最大熵法谱估计结果如下图：
5.附件中表sheet1为某地2008年4月28日凌晨12点至2008年5月4日凌晨12点的电力系统负荷数据，采样时间间隔为1小时，利用Kalman方法预测该地5月5日的电力系统负荷，并给出预测误差（5月5日的实际负荷数据如表sheet2）。
解：卡尔曼滤波是以最小均方误差作为估计的最佳准则，来寻求一套递推估计的算法，其基本思想是：采用信号与噪声的状态空间模型，利用前一时刻地估计值和现在时刻的观测值来更新对状态变量的估计，求得出现时刻的估计值。它适合于实时处理和计算机运算。
现设线性时变系统的离散状态防城和观测方程为：
X(k)=F(k,k-1)X(k-1)+T(k,k-1)U(k-1)
Y(k) = H(k)·X(k)+N(k)
其中：X(k)和Y(k)分别是k时刻的状态矢量和观测矢量，F(k,k-1)为状态转移矩阵，U(k)为k时刻动态噪声，T(k,k-1)为系统控制矩阵，H(k)为k时刻观测矩阵，N(k)为k时刻观测噪声。 卡尔曼滤波的算法流程为： 1、预估计
?=F(k,k-1)·X(k-1) X(k)2、计算预估计协方差矩阵
?=F(k,k-1)×C(k)×F(k,k-1)'+T(k,k-1)×Q(k)×T(k,k-1)' C(k)