西南交通大学信号处理期末作业

1、考虑两个谐波信号()xt和()yt，其中()cos()cxtAwt,()cos()cytBwt式中A和cw为正的常数,为均匀分布的随机变量，其概率密度函数为1,02()20,f其他,而B是一个具有零均值和单位方差的标准高斯随机变量，即其分布函数为21()exp(/2),2Bfbbb（1）求()xt的均值()xut、方差2()xt、自相关函数()xR和自协方差函数()xc。（2）若与B为相互统计独立的随机变量，求()xt和()xt的互相关函数()xyR与互协方差函数()xyc。解：（1）()xt的均值()xut为：220011(())(cos())cos()sin()022xcccuExtEAwtAwtdAwt方差2()xt为：22222222(())(cos())(1cos(22))(cos(22))2222xcccAAAAExtEAwtEwtEwt自相关函数()xR为：22222()(cos()cos(+))(cos()cos(+))(cos(2+22)cos())cos()(cos(2+22))cos()2222xcccccccccccccREAwtAwtwAEwtwtwAAAAEwt自协方差函数()xc为：2()()cos()2xxcAcRw（2）()yt的均值为：()(())(cos())()cos()0yBBccutEytEBwtEBwt,所以()=0EB由互相关函数的定义可知：()(cos()cos())xycccREAwtBwtw由题意知道与B为相互统计独立的随机变量，所以有()(cos()(cos())(cos()()cos())00cos()0xyccccccccREAwtEBwtwAEwtEBwtwAwtw互协方差函数()xyc()()0xyxycR2.接收信号由下式给出：cos(),1,2,...,iciyAiiN，式中~(0,1)iN即i是零均值和单位方差的高斯噪声，c为载波角频率，而是未知的相位。假设12,,...N相互独立，求未知相位的最大似然估计^ML。解：由于12,,...N相互独立，所以1,..Nyy也相互独立并且服从高斯分布，可以得到1,..Nyy与的联合概率密度函数分布21[cos()]211(,..|)(2)NiciyAiNNfyye由此，可以得到似然函数211ln(2)[cos()]22NiciNLyAi该似然函数对求偏导，并令该偏导函数为零，即可得到如下公式：1[cos()]sin()0NicciLyAii因此，最大似然估计^ML为上述函数的零点值。则^^^11cos()sin()sin()NNMLMLMLcciciiAiiyi该函数为非线性方程，不容易求解，若忽略双倍频率2c，则可简化到如下式子：^1sin()0Niciyi根据三角公式分解得到如下式子：^^11sincos()cossin()NNMLicMLiciiyiyi由此，可以得到如下公式^11sintancos()NiciMLNiciyiyi所以相位的最大似然估计如下：^11sin()arctan()cos()NiciMLNiciyiyi3.离散时间的二阶AR过程由差分方程12()(1)(2)()xnaxnaxnwn描述，式中()wn是一零均值、方差为2w的白噪声。证明()xn的功率谱为22212122()12(1)cos(2)2cos(4)wxPfaaaafaf证明：由AR过程的功率谱公式知222412()1xjfjfPfaeae其中2242424121212224422221221212122121221(1)(1)1()()12(1)cos(2)2cos(4)jfjfjfjfjfjfjffjfjfjfjfaeaeaeaeaeaeaaaeeaaeaaeaeeaaaafaf将其带入第一个公式可得：22212122()12(1)cos(2)2cos(4)xPfaaaafaf4、信号的函数表达式为：sin(2100)1.5sin(2300)sin(2200)xtttAttdntnt，其中，At为一随时间变化的随机过程，dnt为经过390-410Hz带通滤波器后的高斯白噪声，nt为高斯白噪声，采样频率为1kHz，采样时间为2.048s。分别利用周期图谱、ARMA、Burg最大熵方法估计信号功率谱，其中ARMA方法需要讨论定阶的问题。解：由题意知采样点数一共为：1000×2.048=2048个数据点。At为一随时间变化的随机过程，由于随机过程有很多类型，如维纳过程、正态随机过程，本文采用了均值为0，方差为1的正态随机过程来作为演示，来代替At，高斯白噪声采用强度为2的高斯白噪声代替nt，其带通滤波后为dnt。其中滤波器采用的是契比雪夫数字滤波器。可得到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)之间总有如下输入输出关系：10()()()pqkkkkxnaxnkbunk0()()()kxnhkunk对以上两个式子两边分别取Z变换，并假定b0=1，可得()()()BzHzAz其中1()1pkkkAzaz,1()1qkkkBzbz,0()()kkHzhkz。为了保证H(z)是稳定的最小相位系统，A(z)和B(z)的零点都应该在单位圆内。假定u(n)是一个方差为2的白噪声序列，由随机信号通过线性系统的理论可知，输出序列X(n)的功率谱为：222**2()()()()()()()jwjwjwjwxjwjwjwBeBeBePeAeAeAeARMA阶数确定：本题目采用AIC准则确定ARMA的阶数。分别计算p、q从1到20阶数的计算出AIC（p,q）,如下图所示，当横坐标大概为230左右时，AIC(p,q)取得最小，将此时的p,q作为带入到模型即可。ARMA法谱估计结果：3、Burg最大熵法Burg算法的具体实现步骤：步骤1计算预测误差功率的初始值和前、后向预测误差的初始值，并令m=1。210)(1NnnxNP)()()(00nxngnf步骤2求反射系数NmnmmNmnmmmngnfngnfK12121111])1()([21)1()(步骤3计算前向预测滤波器系数),()()(11imaKiaiammmm1,...,1mimmKma)(步骤4计算预测误差功率12)1(mmmPKP步骤5计算滤波器输出)1()()(11ngKnfnfmmmm)1()()(11ngnfKngmmmm步骤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、预估计ˆX(k)=F(k,k-1)·X(k-1)2、计算预估计协方差矩阵ˆC(k)=F(k,k-1)×C(k)×F(k,k-1)'+T(k,k-1)×Q(k)×T(k,k-1)'Q(k)=U(k)×U(k)'3、计算卡尔曼增益矩阵K(k)=ˆC(k)×H(k)'×[H(k)׈C(k)×H(k)'+R(k)]-1R(k)=N(k)×N(k)'4、更新估计X(k)=ˆX(k)+K(k)×[Y(k)-H(k)׈X(k)]5、计算更新后估计协方差矩阵C(k)=[I-K(k)×H(k)]׈C(k)×[I-K(k)×H(k)]'+K(k)×R(k)×K(k)'X(k+1)=X(k)C(k+1)=C(k)6、重复以上步骤最终可以获得如下结果：本题将表中的作为观测数据，图中横坐标为1表示2008.4.28日1时刻数据，2表示2008.4.28日2时刻的数据，一次类推，168表示2008.5.5日1时刻的数据。从表中可以看出预测误差的最大值为300。预测误差的大小与代码中的R、Q值得设置有关。Q越大预测误差越小，但是同时也表明系统内的噪声很大。本题中取得Q、R值均为高斯分布的协方差。代码见附录。6.设某变压器内部短路后，故障电流信号分解得到下式：y(t)=20e-tτ⁄+20sin(ωt+60°)+12sin(2ωt+45°)+10sin(3ωt+30°)+6sin(4ωt+22.5°)+5sin(5ωt+36°)式中，ω=2πf0，τ=30ms，f0=50Hz分别利用小波变换、短时傅里叶变换和维格纳威利分布分析故障电流信号的时频特性。解：（1）小波变换：连续小波变换的定义：*,1(,)(),()()()fustuCWTusfttftdtss计算连续时间小波变换的4个步骤：1、选取一个小波，然后将其和待分析信号从起点开始的一部分进行相乘积分。2、计算相关系数c。3、将小波向右移，重复1和2的步骤直到分析完整个信号。4、将小波进行尺度伸缩后再重复1，2，3步骤，直至完成所有尺度的分析。（2）短时傅里叶变换短时傅里叶变换定义如下：,(,)(),()()()itfuSTFTuftgtftgtuedt()1ˆˆ(,)()()()()2itiufSTFTuftgt