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

1、1、考虑两个谐波信号x(t)和y(t),其中Mf) = Acos(wct +, y(t) = Bcos(vv/)式中4和久为正的常数,。为均匀分布的随机变量，其概率密度函数为/) =0,其他而8是一个具有零均值和单位方差的标准高斯随机变量，即其分布函数为fB(b)=4=exp(-/72 / 2), yc b 8Q2兀(1)求W)的均值“)、方差b、自相关函数凡和自协方差函数q(r)0(2)若。与6为相互统计独立的随机变量，求.*)和.*)的互相关函数弋/)与互协方差函 数 J ”)。解：x(f)的均值，为： TOC o 1-5 h z 1 r”1=0o储=E(j(/) = E(Acos(vvf

2、z + 0)= Acos(wet + 族卜/0 =+ ) 2tt 为2兀方差。为:.,.,AAAa； = E(x2(t) = E(A2 cosz(vv.r + 姆)=E(l + cos(2vv./ + 20)= + E(cos(2卬/ + 2。)= 自相关函数6(r)为：/?)A cos(wet+wcr + )= A:E(cos(h;7 +(/) cos(wct+wc T +(/)A 二A2AA=,(cos(2vv;/4-2vvrr + 2 姆 + cos(u;.r)= cos(H；,r) + (cos(2vvrZ+2iv,r + 2)= cos(vvr) 自协方差函数耳(r)为：c*) =

3、R) = cos(ivr)y(f)的均值为：uy (t)=吗(y(7)= EB(Bcos(wct)= EB)cos(vv/) = 0，所以 E(B)=O 由互相关函数的定义可知：R. (r) = E(Acos(h;/ + Q)Bcos(町f- wcr) 由题意知道。与8为相互统计独立的随机变量，所以有R、(t)= *(4cos(w/ + e)X8cos(*f-+T) = 4cos(w/ +e)E(8)cos(w/ 卬/) =A x 0 + 0 x cos(vv/ vv. r) = 0 互协方差函数c,(D%(r)=q(r)=。.接收信号由下式给出：y = Acos(qi + d) + q,i

4、= L2,.,N ,式中弓%。1)即是零均 值和单位方差的高斯噪声，为载波角频率，而8是未知的相位。假设如牡，。相互独立, 求未知相位的最大似然估计。解：由于如在，。相互独立，所以)，.以也相互独立并且服从高斯分布，可以得到 曲.以与8的联合戒率密度函数分布由此，可以得到似然函数N N_L = -ln(2r)- -22yr- 4 cos(g/ + 0)22 i-i该似然函数对e求偏导，并令该偏导函数为零，即可得到如下公式: 丸 n=人 cos(例 i + e)sin(i+e)= O该函数为非线性方程不容易求解,Z A cos(、i + Oml) sin(卬 + 6ML) = Z 凹 sm(、i

5、 + Oml )因此，最大似然估计6M.为上述函数的零点值。 则若忽略双倍频率2 ,则可简化到如下式子:NAZysin(Qj+e)= o根据三角公式分解得到如下式子:A jVsm6ML Z凹 cos(如)=-cos% g)：由此，可以得到如下公式N- Zxsm&itan &ml = Z y cos(Q.i)1-1所以相位的最大似然估计如下：NAZ E sm(ei)般=arctan(Z cosgi).离散时间的二阶AR过程由差分方程x() = 4M-1) +3(-2) +叭)描述，式中贝) 是一零均值、方差为。:的白噪声。证明式)的功率谱为bp(f)=；1 + a； + 6/； 24 (1 -

6、% ) cos(2兀/) - 2a2 cos(41 f )证明：由AR过程的功率谱公式知(/)=-_廿-其中1 -azej4JT, =(-aYeiZ!tf azeJ4irf)(1 -aleJ2ir/ - a?/)=+ a； + a； _ a二(/ 4- e41Tf ) + alazeiZftf + ala2e)2JTf al(ej2r/ + eiZxf )=1 + a：-2a1(l-a2)cos(2/)-2az cos(4f)将其带入第一个公式可得：a2?=；21 + a； + a； - 2% (1 - % ) cos(2兀/) - 2a2 cos(41 f )信 号 的 函 数 表 达 式

7、为：X(/) = sin(2100/) +1.5sm(2300/) + A(t)sin(2200/) + dn(t) + n(t),其中，4(f)为一随时间变 化的随机过程，疝为经过390-410HZ带通滤波器后的高斯白噪声，(，)为高斯白噪声， 采样频率为1kHz,采样时间为2.048s。分别利用周期图谱、ARMA、Buig最大燧方法估计 信号功率谱，其中ARMA方法需要讨论定阶的问题。解：由题意知采样点数一共为：1000 x2.048=2048个数据点。A(f)为一随时间变化的随 机过程，由于随机过程有很多类型，如维纳过程、正态随机过程，本文采用了均值为0,方 差为1的正态随机过程来作为演

8、示，来代替A。)，高斯白噪声采用强度为2的高斯白噪声 代替(/)，其带通滤波后为力小)。其中滤波器采用的是契比雪夫数字滤波器。 可得到x(t)如下图所示：1、周期图法matlab中的周期图功率谱法原理是通过计算采样信号的FFT,获得离散点的幅度，再根 据幅度与功率之间的关系，转换为离散点的功率，再通过坐标变换将离散点的功率图转换为 连续功率谱密度。S/cp/:计算采样信号x(ii)的DFT,使用FFT方法来计算。如果此处将及频率处的幅度对 称到物理实际频率，得到的就是单边谱，否则就是双边谱根据正余弦信号功率与幅度的关系以及直流功率与幅度的关系，将幅度转换为离 散功率谱。对横纵坐标进行转换，横坐

9、标乘以频率分辨率转换为实际连续物理频率，纵坐标 除以频率分辨率转换为功率谱密度。调用 MATLAB 中自带的 matlab 中Pxx.f=periodog】am(x,window.nfitfs)函数可得计算结 果如下：周期图法求功率谱 TOC o 1-5 h z 21.5 -玄0050100150200250300350400450500f/fs5.附件中表sheet1为某地2008年4月28 口凌晨12点至2008年5月4 口凌晨12点的 电力系统负荷数据，采样时间间隔为1小时，利用Kalman方法预测该地5月5 口的电力系 统负荷，并给出预测误差(5月5 口的实际负荷数据如表sheet2)

10、。解：卡尔曼滤波是以最小均方误差作为估计的最佳准则，来寻求一套递推估计的算法, 其基本思想是：采用信号与噪声的状态空间模型，利用前一时刻地估计值和现在时刻的观测 值来更新对状态变量的估计，求得出现时刻的估计值。它适合于实时处理和计算机运算。现设线性时变系统的离散状态防城和观测方程为：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时刻观测噪声。

11、 卡尔曼滤波的算法流程为： 1、预估计X(k) =F(kl)X(k-l)2、计算预估计协方差矩阵e(k)=F(k,k-l)xC(k)xF(k,k-D+T(k,k-l)xQ(k)xT(k,k-l)，Q(k)=U(k)xU(kY3、计算卡尔曼增益矩阵K(kC(k) xH(k)rxH(k)x C(k) xH(k)f+R(k)-R(k尸N(k)xN(k)f4、更新估计X(k)=X(k) +K(k)xY(k)-H(k)x X(k)5、计算更新后估计协方差矩阵C(k) = LK(k)xH(k)x C(k) xi-K(k)xH(k)T+K(k)xR(k)xK(k)，X(k+l)=X(k)C(k+l)=C(k

12、)6、重更以上步骤最终可以获得如下结果：使用Kalrmm对电力系统负荷数据进行颈测35003000250020001500100050020406080100120140160180时间点数350030002500使用Kalman对电力系统负荷数据进行预测300200 1000 100预测值与真实值之误差168170172174176178180182184186188190192时间点数本题将表中的作为观测数据，图中横坐标为1表示2008.4.28 1时刻数据,2表示2008.4.28 口 2时刻的数据，一次类推，168表示200855 口 1时刻的数据。从表中可以看出预测误差 的最大值为3

13、00。预测误差的大小与代码中的R、Q值得设置有关。Q越大预测误差越小， 但是同时也表明系统内的噪声很大。本题中取得Q、R值均为高斯分布的协方差。代码见附录。6,设某变压器内部短路后，故障电流信号分解得到卜.式:y(0= 20et/x + 20sin (ent + 60 ) + 12sin(2cnt + 45 ) + lOsin(3(nt + 30 ) + 6sin(4(ot + 22.5 ) + 5sin(5wt + 36 )式中，fn = 27ll0, T 30ms,o = 5OH/分别利用小波变换、短时傅里叶变换和维格纳威利分布分析故障电流信号的时频特性。解：（1）小波变换：连续小波变换的

14、定义：。啊(“,s) =匚/(0 炉(彳)力计算连续时间小波变换的4个步骤:选取一个小波，然后将其和待分析信号从起点开始的一部分进行相乘积分。 计算相关系数C。将小波向右移，重复U和2的步骤直到分析完整个信号。将小波进行尺度伸缩后再重及1, 2, 3步骤，直至完成所有尺度的分析。（2）短时傅里叶变换短时傅里叶变换定义如下:STF/ (,) = /(/), g“=匚 /(f)gQ - )产力STFTf （,9=/OgQ37（助（出一9*吐文/口J24j（3）维格纳威利分布变换 维格纳威利分布定义如下：WDt（/.Q）= （： / + 弁r在MATLAB中没有维格纳威利分布变换的相关函数，需要安装

15、一个MATLAB版本的时频 分析工具箱。调用里面的函数即可。小波变换和短时傅里叶变换MATLAB均自带了相关的 函数。程序见附录。代码运行结果结果如下：小波时频图501001502002500.10.20,3040.50.60.70 80.9时间t/S70605040302010妞时停单叶交换结果50000.10.20.30.40.50.60.70.80.91时间t/s7.假定一电力系统谐波与间谐波信号的函数表达式如下：y () = 0.001 cos(2 乃 x 10 + 疚)+ cos(2 乃 x 50 + 4)+ 0.1 cos(2;r x 150 + 内)+ 0.002 cos (2乃

16、 x 50 + 四)+()其中，采样频率为1024Hz，相位-内为独立的均匀分布-,+;4)为一噪声信号，信噪比取为20dB ,分别采用三种现代信号处理方法进行谐波与间谐波频率提取与谱估计。解:本题目采用的频率提取的三种方法为小波变换、短时傅里叶变换和维格纳威利分布。采用周期图法、MUSIC法、Burg法进行谱估计。确定出谐波的频率为50Hz和150Hz。原始侑号200.20,4060.811.21.41.61.8时向（U秒） 小波时频图00.20.40.60.811.21.41.61.8时间t/S 500 400 300 200 100 0Wignor-Villo timo-froquonc

17、y distribution500450400350z 300三2502001501OO500.20.40.60.811.2位同Vs1.41.61.80.5周期图法求功率谱4 3 2 1 o o o oqp、w-a50100 150 200 250 300 350 400 450 500 550 600f/HzMUSIC方法o o o o4 2 2- (gp)与-40050100150 200 250 300350 400 450500频率(f/Hz)Burg法谱估计 0.0150.010.0050 050 100 150 200 2500060附录代码：第四题：clc;clear;fs=10

18、00;%采样频率T=2.048;%采样时间t=0:l/fs:T;A = normmd(0,l JJength(t);%方差为 1,均值为。的高斯分布N=wgn(l ,length(t),2);%强度为2的高斯白噪 声300 350 400 450 500 550 600f/fsDn=bandp(N,390,410,200,450,0.1,30.fs);figure。)；subplot(211);plot(t,N);title(,原始高斯白噪声工subplot(212);plot(t.Dn);titleC带通滤波后高斯白噪声)；Sig=siii(2*pi* 100. *t)+1.5 *sin(2

19、 *pi*300.*t)+A.*siii(2 *pi*200.*t)+Dn+N;figure。)；plot(t,Sig);titlef原始输入信号)axis(O 2.1 -7 7);%周期图谱Pxx,f=periodogiam(Sig,lengtli(t),fs);% 周 期图法figuie(3);plot(f,Pxx);出le(，周期图法求功率谱工xlabelCf/Hz); ylabelC 功率/db);% ARMA谱估计z=iddata(Sig);%#信号转化为matlab接受的 格式RecoidAIC=;fbrp=l:20%自回归对此PACF,给定滞后长 度上限p和qfor q=l:20

20、%移动平均对应ACF m=ainiax(z(l :lengtli(t).p,q);AIC = aic(m); %armax(p,q)选择对 应FPE最小，AIC值最小模型RecordAIC=RecordAIC;p q AIC;end endfbr k=l:size(RecoidAIC,l)ifRecoidAIC(k,3)=niui(RecoidAIC(：3) %选 择AIC最小模型pa_AIC=RecordAIC(k, 1);qa_AIC=RecordAIC(k,2); break;end endmAIC=aimax(z( 1 :length(t). pa_AIC.qa_AIC );Pxx2,

21、f2=freqz(niAIC.C4iiAIC.a,fs);P2=(abs(Pxx2).* 1).A2;P2tol=10*logl0(P2);figuie(4);plot(f2/pi*fs/2,P2tol); title(ARMA 法(AIC 准 则)Xabul(Hz);ylabulC 振幅/dB);plot(RecordAIC(: ,3);ylabel(AIC(p,q),);% burg法计算Pxx,F = pburg(Sig,60,length(t),fs);%burg 法 figure(6);plot(EPxx);titlefBuig法谱估计工xlabel(,f/fs,);%X 轴坐标名称

22、ylabelC功率谱/dB);%Y轴坐标名称 %fiinction y=bandp(x,f Lf3,fsl,Ah,ip,rs.Fs)%带通滤波%使用注意事项：通带或阻带的截止频率与 采样率的选取范围是不能超过采样率的一 半%即,fl,B,fsl,fsh,的值小于 Fs/2%x:需要带通滤波的序列% f 1 ：通带左边界% f 3：通带右边界%fsl：衰减截止左边界% fsh：衰变截止右边界%ip：边带区衰减DB数设置%rs：截止区衰减DB数设置%FS：序列x的采样频率% fl=3OO;f3=5OO;%通带截止频率上下限% %=200；fsh=600;%阻带截止频率上下限% ip=0.1;rs=

23、30;%通带边衰减DB值和阻带 边衰减DB值%Fs=2000;% 采样率 %wpl=2*pi*fl/Fs;wp3=2*pi*f3/Fs;wsl=2*pi*fsLTs;wsh=2*pi*fsh/Fs;wp=wpl wp3;ws=wsl wsh;%设计切比雪夫滤波器；n,wn=cheb 1 ord(ws/pi,wp/pi,rps);bz 1 ,az 1 =cheby 1 (n,rp,wp/pi);%查看设计滤波器的曲线hjv=fieqz(bz 1 ,az 1,256,Fs);h=20*log 10(abs(h);y=filter(bzl.azLx);end第5题%本题目需要提醒一点：给的数据为观测

24、数 据Z而不是X clc;clear;xl=xlsiead(/负荷数据.xlsYshee); xl=xl(:,2);x2=xlsieadC./负荷数据.xls?sheet2);x2=x2(:,2);x=xl;x2;Nl=length(xl);N=length(x);A=l；B=0;H=l；w=nomund(0,1000,1 ,N);% 这里随便取值 v=noimrnd(0.1000,1 ,N);P(1尸16;%随便取值Z=x;X(l)=24;%随便取值R=cov(v);Q=cov(w);fbr 1=2 :Ntempx=A*X(i-1 );%+B *u ； TempP=A*P(i-1 )* A+

25、Q;K(i)=TempP*H,* 1 /(H*TempP*Hr+R); X(i)=X(i-1 )+K(i)*(Z(i)-tempx);P(i)=( 1 -K(i)*H)*TempP;endt=l:length(Z);figure;plot(t,Z,b;t,X(t),T);titleC使用 Kahuan对电力 系统负荷数据进行预测工xlabelf时间点数);ylabe(电力系统负荷 *);axis tight;legend(，负荷真实值？Kalman 预 测值工 figure;subplot(2,l,l);t=length(x 1): length(x);使用 Kalman 对电 力系统负荷数据

26、进行预测力xlabelC时间点数 );ylabel(电力系统负荷力axis tight;legend(,负 荷真实值？Kalman预测值)； set(gca,XTick,length(xl):2:length(x); subplot(2,l,2);enoi-Z-X*;plot(t,enor(t);titleC预测值与真实值之误差 力xlabelC时间点数) set(gca/XTicklength(x l):2:length(x);ylabel(,5月5 口预测值与真实值误差axis tight;第六题：%小波变换clc;clear;close all;f=50;%信号频率ouniiga=2*p

27、i*f;N_sample=2048;%总采样点数Fs=1000;%采样频率 t=O:l/Fs:l;Tao=0.03;A=l;%信号幅度X=20*exp(-t/Tao)+20*sm(ouniiga*t+pi/3)+12*si n(2 *ounuga *t+pi/4)+10*siii(3 *ounuga*t+piz6 )+6*sin(4*oumiga *t+pi/8)+5 * sm(5 *ouiniga *t +pi/5); % 信号函数表达式figure;plot(t,x);titl原始信号工xlabelfWfu t/sTontSize14);ylabel(幅值；FontSiz-M);%原信号函

28、数wavenaine=,cmor3 -3 totalscal=256;Fc=centfrq(waveiiaine); %小波中心频率 c=2*Fc*totalscal;scals=c./(l :totalscal);fscal2fiq(scals,wavename,l/Ts); % 将尺度 转换为频率coefs=cwl(x,scals,wavename); % 求连续小 波系数figure;unagesc(tXabs(coefs);colorbar;xlabelfWfu t/sTontSize14);ylabel(频率 f HzTontSize 14);title (,小波时频图；FontSi

29、ze, 16);axis(0 1 0 300);%短时傅里叶变换S,F,T,P=spectiogiam(x,256,250,256,Fs);figure;surf(TE10*logl0(P)/edgecolor7none,)； axis tight;view(0,90);xlabelfRt 间/s); ylabeiq频率/Hz);titleC短时傅里叶变换结果力% Wigner-Ville time-frequency distnbution. X=hilbert(xf);scals=c./(l :totalscal);fscal2fiq(scals,wavenameJ/Ts); % 将尺度

30、转换为频率coefs=cwt(x 1,seals,wavename); % 求连续小波系数figure;unagesc(n,f,abs(coefs);colorbar;xlabelC时间 t/s,FontSize, 14);ylabel(濒率 f;HzTontSize 14);titleC 小波时频图；FontSize： 16);figure;contour(t/Fs,#Fs,abs(tft);xlabelC时间 t/sr);ylabelC频率 fHz*);title(,Wigner-Villetiine-fiequencydistribution); axis(0 1 0 300) % 第七题： clc;clear; close all;%参数设置Fs= 1024; %采样频率 n = 0:LFs201;%采样时间 N = length(n); % 采样点Wl=0.001 *cos(2 *pi*n*l O+unifhid(-pi,pi)+c os(2*pi*50*n+unifrnd(-pi,pi)+0.1 *cos(2*pi*