解：
C1）
2二
=0
0
x（t）的均值Ux（t）为：
12応1
山=E（x（t））=E（Acos（Wct+©））=云『Acos（Wct+©）d©Asin（Wet+0）
方差；「x（t）为:
2222g=E（x2（t））二E（A2cos2（wct•J）=AE（1-cos（2wct2J）=——E（cos（2wct-2'））=一
2222自相关函数Rx（）为：
&（.）=E（Acos（wct4）Acos帆t+wc?
：
「；））=A2E（cos（wctq）cos（wct+wc2：
；））
一一一一
E（cos（2wct+2wcx"2）cos（wc.））cos（wc.）E（cos（2wct+2w^_2））cos（wc.）
2222
自协方差函数c<（）为：
A2
q（.）=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（WctJBcos（Wct「Wc.））
由题意知道••与B为相互统计独立的随机变量，所以有
Rxy（0=E（Acos（wct•）E（Bcos（wct—wc.））=AE（cos（wct•）E（B）cos（wct—wc.））=A00cos（wct—wc.）=0
互协方差函数cxy（.）
2•接收信号由下式给出:
值和单位方差的高斯噪声，
Gy（）二Rxy（）=0
y二Acos（ci二）d=12…，N，式中】~N（0,1）即r是零均
.c为载波角频率，而是未知的相位。
假设「J6…去相互独立，
求未知相位的最大似然估计TU。
解：
由于-1/'2,../'N相互独立，所以yi,..yN也相互独立并且服从高斯分布，可以得到
yi,..yN与二的联合概率密度函数分布
1龙
f（y1,..yNI&）=（松）ne*
由此，可以得到似然函数
N1N2
LIn（2二）[yi—Acos（，ci•R]
22izt
该似然函数对二求偏导，并令该偏导函数为零，即可得到如下公式：
一=E[y—Acosfeici+6）]sing+8）=0
:
-Via
A
因此，最大似然估计JML为上述函数的零点值。
则
NAANA
Acos（訂^ML）sin（fl）='ysin（fl）
i-1i1
该函数为非线性方程，不容易求解，若忽略双倍频率2,c，则可简化到如下式子:
NA
Eysin（mJ+8）=0
i土根据三角公式分解得到如下式子：
ANAN
sinTml送ycos⑥J）=一costU瓦ysin（轨i）
i2i1
由此，可以得到如下公式一一
N
Ayisin•ci
tan寸ml=_
yicos（・ci）
iz1所以相位的最大似然估计如下：
N
A'ysin（「ci）
^ml二-arctan（）
yicos（ci）
i土
3•离散时间的二阶AR过程由差分方程x（n）二ax（n—1）*2X（n「2）w（n）描述，式中w（n）
■a2-a2
是一零均值、方差为；「W的白噪声。
证明x（n）的功率谱为
—2a（1—a2）cos（2二f）—2a2cos（4二f）
证明：
由AR过程的功率谱公式知
巳（f）
-a?
e
其中
F-aei—a2ej4”f）（1-ae*fpe"）
=1q2・a2-a2（e*：
f飞口）乜啓心飞怒中_a（e*：
f飞心）
22
=1aia2-2a,（1-a2）cos（2二f）-2a?
cos（4二f）
将其带入第一个公式可得：
Px（f）2-
1+a+a2-2ai（1-a?
）cos（2兀f）一2a?
cos（4兀f）
4、信号的函数表达式为:
xtj=sin（2二100t）1.5sin（2二300t）Atsin（2二200t）dntnt，其中，At为一随时间变化的随机过程，dnt为经过390-410HZ带通滤波器后的高斯白噪
声，nt为高斯白噪声，采样频率为1kHz，采样时间为2.048s。
分别利用周期图谱、ARMA、Burg最大熵方法估计信号功率谱，其中ARMA方法需要讨论定阶的问题。
解：
由题意知采样点数一共为：
1000X2.048=2048个数据点。
At为一随时间变化的随
机过程，由于随机过程有很多类型，如维纳过程、正态随机过程，本文采用了均值为0,方
差为1的正态随机过程来作为演示，来代替At，高斯白噪声采用强度为2的高斯白噪声
代替nt，其带通滤波后为dnt。
其中滤波器采用的是契比雪夫数字滤波器。
可得到x（t）如下图所示：
1、周期图法
matlab中的周期图功率谱法原理是通过计算采样信号的FFT，获得离散点的幅度，再根
据幅度与功率之间的关系，转换为离散点的功率，再通过坐标变换将离散点的功率图转换为连续功率谱密度。
Stepl:
计算采样信号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）之间总有如下输入输出关系：
pq
x（n）=-、akx（n—k）亠二Ru（n—k）
k±kzfi：
x（n）-7h（k）u（n-k）
对以上两个式子两边分别取
k_0
Z变换，并假定bo=l，可得
H（z）罟
A（z）
为:
ARMA阶数确定：
本题目采用AIC准则确定ARMA的阶数。
分别计算p、q从1到20阶数的计算出AIC（p,q）,女吓图所示，当横坐标大概为230左右时，AIC（p,q）取得最小，将此时的p,q作为带入到模
型即可。
ARMA法谱估计结果:
3C
2b
3、Burg最大熵法
Burg算法的具体实现步骤：
步骤1计算预测误差功率的初始值和前、后向预测误差的初始值，并令m=1。
1N2Po石'x（n）
Nn土
fo（n）=g°（n）=x（n）
步骤2求反射系数
-扛fm£（n）gml（n-1）
_
m一N2~-[fm丄（n）-gm2（n-1）]
2n_m1
步骤3计算前向预测滤波器系数
步骤4计算预测误差功率
步骤5计算滤波器输出
2
Pm二（1_Km|）Pm丄
fm（n）二fm_1（n）■Kmgm^（n-1）gm（n）*£丄（n）gm」（n-1）
步骤6令m—m+1，并重复步骤2至步骤5，直到预测误差功率Pm不再明显减小。
最后，再利用Levinson递推关系式估计AR参数，继而得到功率谱估计
0M1001SG200260350400450SOO
唯
5•附件中表sheetl为某地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、预估计
乂（k）=F（k,k-1）X（k-1）
2、计算预估计协方差矩阵
C（k）=F（k,k-1）C（k）采（k,k-1）'+T（k,k-1）Q（k）XT（k,k-1）'
Q（k）=U（k）U（k）'
3、计算卡尔曼增益矩阵
A
K（k）=C（k）>H（k）'[H（k）滋（k）XH（k）'+R（k）]-
R（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）]x?
（k）>[I-K（k）H（k）]'+K（k）R（k）>K（k）'
X（k+1）=X（k）
C（k+1）=C（k）
6、重复以上步骤
最终可以获得如下结果：
岐用也Im胡对鹿力系推负花数克逬廿硕遐
图中横坐标为1表示2008428日1时刻数据，2表示2008428
168表示2008.5.5日1时刻的数据。
从表中可以看出预测误差
的最大值为300。
预测误差的大小与代码中的R、Q值得设置有关。
Q越大预测误差越小,
但是同时也表明系统内的噪声很大。
本题中取得Q、R值均为高斯分布的协方差。
代码见附
录。
6•设某变压器内部短路后，故障电流信号分解得到下式：
分别利用小波变换、短时傅里叶变换和维格纳威利分布分
式中，，
析故障电流信号的时频特性。
解：
（1）小波变换：
连续小波变换的定义：
计算连续时间小波变换的4个步骤：
1、选取一个小波，然后将其和待分析信号从起点开始的一部分进行相乘积分。
2、计算相关系数c。
3、将小波向右移，重复1和2的步骤直到分析完整个信号。
4、将小波进行尺度伸缩后再重复1，2，3步骤，直至完成所有尺度的分析。
（2）短时傅里叶变换
短时傅里叶变换定义如下：
STFTf（u,'）=：
■：
.f（t）,gu,（t^；'.'f（t）g（t-u）e±tdt
STFTf（u,）=-f（t）g（t-u）e丄Pt二丄：
？
（J?
（.-）eiu^）d.
2
（3）维格纳威利分布变换维格纳威利分布定义如下：
WDx（tQ）=丘xt-x*e叮d.
在MATLAB中没有维格纳威利分布变换的相关函数，需要安装一个MATLAB版本的时频
分析工具箱。
调用里面的函数即可。
小波变换和短时傅里叶变换MATLAB均自带了相关的
函数。
程序见附录。
代码运行结果结果如下：
IIIIII11II
-ro
ddudJ0-40.5&7o.a>柑上
时间诒
•SCSI
4L-J
7•假定一电力系统谐波与间谐波信号的函数表达式如下：
yn=0.001cos（2c.：
1On••cos（2x50n亠「2），0.1cos（2H50n亠：
；3）0.002cos2二50n「％厂打n
其中，采样频率为1024Hz，相位丄'4为独立的均匀分布U丨-二，；I；'n为一噪声信号，信噪比取为20dB。
分别采用三种现代信号处理方法进行谐波与间谐波频率提取与谱估计。
解：
本题目采用的频率提取的三种方法为小波变换、短时傅里叶变换和维格纳威利分布。
采用周期图法、MUSIC法、Burg法进行谱估计。
确定出谐波的频率为50Hz和150Hz。
岸Q
O.ft11.2
时
小波刘频图
5W
45D
3OT
咖
1M
0.J
n4
D.fl
Wi（jnwf8-Vill«tirnv-Fraqu・iiu#cli»fribwlion
讥a
400
400
250
2C0
1*30
1QO
50
Q
D.2
OBI
1L21.41.B
1fl
0.5
4
-2
0_
-20
-40
050100150200250300350400450500
频
BurgftiH怙计
0.015
N=wgn（1,length（t）,2）;%强度为2的高斯白噪
1，均
附录代码：
第四题：
clc;
clear;
fs=1000;%采样频率
T=2.048;%采样时间
t=0:
1/fs:
T;
A=normrnd（0,1,1,length（t））;%方差为
值为0的高斯分布
Dn=bandp（N,390,410,200,450,0.1,30,fs）;figure
（1）;
subplot（211）;
plot（t,N）;
title（'原始高斯白噪声'）;subplot（212）;
plot（t,Dn）;
title（'带通滤波后高斯白噪声’）；
Sig=sin（2*pi*100.*t）+1.5*sin（2*pi*300.*t）+A.*sin（2*pi*200.*t）+Dn+N;
figure
（2）;
plot（t,Sig）;
title（'原始输入信号'）;
axis（[O2.1-77]）;
%%周期图谱
[Pxx,f]=periodogram（Sig,[],length（t）,fs）;%周期图法
figure（3）;
plot（f,Pxx）;
title（'周期图法求功率谱'）;xlabel（'f/Hz'）;ylabel（'功率/db'）;
%%ARMA谱估计
z=iddata（Sig'）;%将信号转化为matlab接受的格式
RecordAIC=[];
forp=1:
20%自回归对应PACF，给定滞后长度上限p和q
forq=1:
20%移动平均对应ACFm=armax（z（1:
length（t））,[p,q]）;AIC=aic（m）;%armax（p,q）选择对应FPE最小，AIC值最小模型
RecordAIC=[RecordAIC;pqAIC];end
end
fork=1:
size（RecordAIC,1）
if
RecordAIC（k,3）==min（RecordAIC（:
3））%选
择AIC最小模型
pa_AIC=RecordAIC（k,1）;qa_AIC=RecordAIC（k,2）;break;
end
end
mAIC=armax（z（1:
length（t））,[pa_AIC,qa_AIC]
）;
[Pxx2,f2]=freqz（mAIC.c,mAIC.a,fs）;
P2=（abs（Pxx2）.*1）.A2;
P2tol=10*log10（P2）;
figure（4）;
plot（f2/pi*fs/2,P2tol）;title（'ARMA法（AIC准则）'）;xlabel（'f/Hz'）;ylabel（'振幅/dB'）;plot（RecordAIC（:
3））;
ylabel（'AIC（p,q）'）;
%%burg法计算
[Pxx,F]=pburg（Sig,60,length（t）,fs）;%burg法figure（6）;
plot（F,Pxx）;
title（'Burg法谱估计'）;
xlabel（'f/fs'）;%X轴坐标名称
ylabel（'功率谱/dB'）;%Y轴坐标名
称
%%
functiony=bandp（x,f1,f3,fsl,fsh,rp,rs,Fs）
%带通滤波
%使用注意事项：
通带或阻带的截止频率与采样率的选取范围是不能超过采样率的一半
%即，f1,f3,fs1,fsh,的值小于Fs/2
%x:
需要带通滤波的序列
%f1:
通带左边界
%f3:
通带右边界
%fs1:
衰减截止左边界
%fsh:
衰变截止右边界
%rp:
边带区衰减DB数设置
%rs:
截止区衰减DB数设置
%FS:
序列x的采样频率
%f1=300;f3=500;%通带截止频率上下限
%fsl=200;fsh=600;%阻带截止频率上下限
%rp=0.1;rs=30;%通带边衰减DB值和阻带边衰减DB值
%Fs=2000;%采样率
%
wp1=2*pi*f1/Fs;
wp3=2*pi*⑶Fs;
wsl=2*pi*fsl/Fs;
wsh=2*pi*fsh/Fs;
wp=[wp1wp3];
ws=[wslwsh];
%
%设计切比雪夫滤波器；
[n,wn]=cheb1ord（ws/pi,wp/pi,rp,rs）;
[bz1,az1]=cheby1（n,rp,wp/pi）;
%查看设计滤波器的曲线
[h,w]=freqz（bz1,az1,256,Fs）;
h=20*log10（abs（h））;
y=filter（bz1,az1,x）;
end
第5题
%本题目需要提醒一点：
给的数据为观测数
据Z而不是X
clc;
clear;
x1=xlsread（'./负荷数据.xls','sheet1'）;x1=x1（:
2）;
x2=xlsread（'./负荷数据.xls','sheet2'）;x2=x2（:
2）;
x=[x1;x2];
N1=length（x1）;
N=length（x）;
A=1;
B=0;
H=1;
w=normrnd（0,1000,1,N）;%这里随便取值v=normrnd（0,1000,1,N）;
P
（1）=16;%随便取值
Z=x;
X
（1）=24;%随便取值
R=cov（v）;
Q=cov（w）;
fori=2:
N
tempx=A*X（i-1）;%+B*u（i）;
TempP=A*P（i-1）*A'+Q;
K（i）=TempP*H'*1/（H*TempP*H'+R）;X（i）=X（i-1）+K（i）*（Z（i）-tempx）;
P（i）=（1-K（i）*H）*TempP;
end
t=1:
length（Z）;
figure;
plot（t,Z,'b',t,X（t）,'r'）;title（'使用Kalman对电力系统负荷数据进行预测'）;
xlabel（'时间点数'）;ylabel（'电力系统负荷'）;axistight;legend（'负荷真实值','Kalman预测值'）；
figure;
subplot（2,1,1）;
t=length（x1）:
length（x）;
plot（t,x（t）,'b',t,X（t）,'r'）;title（'使用Kalman对电力系统负荷数据进行预测'）;xlabel（'时间点数'）;ylabel（'电力系统负荷'）;axistight;legend（'负荷真实值','Kalman预测值'）;
set（gca,'XTick',length（x1）:
2:
length（x））;subplot（2,1,2）;
error=Z-X';
plot（t,error（t））;title（'预测值与真实值之误差'）;xlabel（'时间点数'）;
set（gca,'XTick',length（x1）:
2:
length（x））;
ylabel（'5月5日预测值与真实值误差'）;axistight;
第六题：
%%小波变换
clc;
clear;
closeall;
f=50;%信号频率
oumiga=2*pi*f;
N_sample=2048;%总采样点数
Fs=1000;%采样频率
t=0:
1/Fs:
1;
Tao=0.03;
A=1;%信号幅度
x=
20*exp（-t/Tao）+20*sin（oumiga*t+pi/3）+12*sin（2*oumiga*t+pi/4）+10*sin（3*oumiga*t+pi/6
）+6*sin（4*oumiga*t+pi/8）+5*sin（5*oumiga*t+pi/5）;%信号函数表达式
figure;
plot（t,x）;
title（'原始信号'）;
xlabel（'时间t/s','FontSize',14）;
ylabel（'幅值','FontSize',14）;
%原信号函数
wavename='cmor3-3';
totalscal=256;
Fc=centfrq（wavename）;%小波中心频率c=2*Fc*totalscal;
scals=c./（1:
totalscal）;
f=scal2frq（scals,wavename,1/Fs）;%将尺度转换为频率
coefs=cwt（x,scals,wavename）;%求连续小
波系数
figure;
imagesc（t,f,abs（coefs））;colorbar;
xlabel（'时间t/s','FontSize',14）;ylabel（'频率f/Hz','FontSize',14）;title（'小波时频图','FontSize',16）;axis（[010300]）;
%%短时傅里叶变换
[S,F,T,P]=spectrogram（x,256,250,256,Fs）;figure;
surf（T,F,10*log10（P）,'edgecolor','none'
升级会员 