1、西南交通大学信号处理期末作业1、考虑两个谐波信号 x(t)和y(t)，其中x(t) =Acos(Wct), y(t) =Bcos(Wct)式中A和Wc为正的常数，为均匀分布的随机变量，其概率密度函数为1,0 _ _2二f ( ) 2- ,9其他而B是一个具有零均值和单位方差的标准高斯随机变量，即其分布函数为1 2fB(b) exp(_b /2), -二:b :(1 )求x(t)的均值Ux(t)、方差c2x(t)、自相关函数Rx( )和自协方差函数G(.)。(2)若与B为相互统计独立的随机变量， 求x(t)和x(t)的互相关函数Rxy( )与互协方差函数 cy(.)。解：C 1)2 二=00x(
2、t)的均值Ux(t)为：1 2応 1山=E(x(t) =E(Acos(Wct +)=云Acos(Wct +)d Asin(Wet +0)方差；x(t)为:2 2 2 2 g =E(x2(t)二E(A2cos2(wct J) =AE(1 - cos(2wct 2 J)= E(cos(2wct - 2 )=一2 2 2 2 自相关函数Rx()为：&( .) =E(Acos(wct 4) A cos 帆 t+wc? ：；)= A2 E(cos(wct q)cos(wct+wc2 ：；)一 一 一 一E(cos(2wct+2wc x 2 ) cos(wc .) cos(wc.) E(cos(2wct+
3、2w _2 ) cos(wc.)2 2 2 2自协方差函数cH(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表示2008428168表示2008.5.5日1时刻的数据。从表中可以看出预
4、测误差的最大值为300。预测误差的大小与代码中的 R、Q值得设置有关。Q越大预测误差越小,但是同时也表明系统内的噪声很大。 本题中取得Q、R值均为高斯分布的协方差。 代码见附录。6设某变压器内部短路后，故障电流信号分解得到下式：分别利用小波变换、短时傅里叶变换和维格纳威利分布分式中， ，析故障电流信号的时频特性。解：(1)小波变换：连续小波变换的定义：计算连续时间小波变换的 4个步骤：1、 选取一个小波，然后将其和待分析信号从起点开始的一部分进行相乘积分。2、 计算相关系数c。3、 将小波向右移，重复 1和2的步骤直到分析完整个信号。4、 将小波进行尺度伸缩后再重复 1，2，3步骤，直至完成所
5、有尺度的分析。 (2 )短时傅里叶变换短时傅里叶变换定义如下：STFTf(u, )=：. f(t),gu, (t ；.f(t)g(t-u)etdtSTFTf (u, )= - f(t)g(t -u)e丄Pt 二丄 ：？( J?( . - )eiu)d .2(3)维格纳威利分布变换 维格纳威利分布定义如下：WDx(tQ)=丘x t - x* e 叮d.在MATLAB中没有维格纳威利分布变换的相关函数，需要安装一个 MATLAB版本的时频分析工具箱。调用里面的函数即可。小波变换和短时傅里叶变换 MATLAB均自带了相关的函数。程序见附录。代码运行结果结果如下：III II I 11 I I-rod
6、 du dJ 0-4 0.5 & 7 o.a 柑上时间诒SCSI4L-J7假定一电力系统谐波与间谐波信号的函数表达式如下：y n =0.001cos（2 c. ：1 On cos（2x 50n 亠2），0.1cos（2H50n 亠：；3） 0.002cos 2二 50n 厂打 n其中，采样频率为1024Hz，相位 丄4为独立的均匀分布 U丨-二，；I； n为一噪声信号， 信噪比取为20dB。分别采用三种现代信号处理方法进行谐波与间谐波频率提取与谱估计。解：本题目采用的频率提取的三种方法为小波变换、 短时傅里叶变换和维格纳威利分布。采用周期图法、MUSIC法、Burg法进行谱估计。确定出谐波的频
7、率为 50Hz和150Hz。岸QO.ft 1 1.2时小波刘频图5W45D3OT咖1M0.Jn 4D.flWi (jn wf8-Vill tirnv-Fraquiiu# clif ribwlion讥a4004002502C01*301QO50QD.2OBI1L2 1.4 1.B1 fl0.54-20_-20-400 50 100 150 200 250 300 350 400 450 500频BurgftiH 怙计0.015N=wgn(1,length(t),2);%强度为2的高斯白噪1，均附录代码：第四题：clc;clear;fs=1000;%采样频率T=2.048;%采样时间t=0:1/f
8、s:T;A = normrnd(0,1,1,length(t);% 方差为值为0的高斯分布Dn=ba ndp(N,390,410,200,450,0.1,30,fs); figure(1);subplot(211);plot(t,N);title(原始高斯白噪声); subplot(212);plot(t,D n);title(带通滤波后高斯白噪声)；Sig=si n(2*pi*100.*t)+1.5*si n(2*pi*300.*t)+ A.*si n( 2*pi*200.*t)+D n+N;figure(2);plot(t,Sig);title(原始输入信号);axis(O 2.1-7 7
9、);%周期图谱Pxx,f=periodogram(Sig,le ngth(t),fs);% 周 期图法figure(3);plot(f,Pxx);title(周期图法求功率谱); xlabel(f/Hz); ylabel(功率 /db);% ARMA谱估计z=iddata(Sig);%将信号转化为matlab接受的 格式RecordAIC=;for p=1:20 %自回归对应PACF，给定滞后长 度上限p和qfor q=1:20%移动平均对应 ACF m=armax(z(1:le ngth(t),p,q); AIC = aic(m); %armax(p,q)选择对 应FPE最小，AIC值最小模
10、型RecordAIC=RecordAIC;p q AIC; endendfor k=1:size(RecordAIC,1)ifRecordAIC(k,3)=mi n(RecordAIC(:,3) % 选择AIC最小模型pa_AIC=RecordAIC(k,1); qa_AIC=RecordAIC(k,2); break;endendmAIC=armax(z(1:le ngth(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
11、*fs/2,P2tol); title(ARMA 法(AIC 准 则);xlabel(f/Hz);ylabel(振幅 /dB); plot(RecordAIC(:,3);ylabel(AIC(p,q);% burg法计算Pxx,F = pburg(Sig,60,le ngth(t),fs);%burg 法 figure(6);plot(F,Pxx);title(Burg 法谱估计);xlabel(f/fs); %X 轴坐标名称ylabel(功率谱/dB); %Y轴坐标名称%fun cti on y=ba ndp(x,f1,f3,fsl,fsh,rp,rs,Fs)%带通滤波%使用注意事项：通带或
12、阻带的截止频率与 采样率的选取范围是不能超过采样率的一 半%即，f1,f3,fs1,fsh,的值小于 Fs/2%x:需要带通滤波的序列% f 1 :通带左边界% f 3 :通带右边界% 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
13、*pi*fsl/Fs;wsh=2*pi*fsh/Fs;wp=wp1 wp3;ws=wsl wsh;%设计切比雪夫滤波器；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而不是Xclc;clear;x1=xlsread(./ 负荷数据.xls,sheet1); x1=x1(:,2);x2=xlsread(./ 负荷数据.x
14、ls,sheet2); x2=x2(:,2);x=x1;x2;N1= le ngth(x1);N=le ngth(x);A=1;B=0;H=1;w=normrnd(0,1000,1,N);% 这里随便取值 v=n ormrnd(0,1000,1,N);P(1)=16;%随便取值Z=x;X(1)=24;%随便取值R=cov(v);Q=cov(w);for i=2:Ntempx=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
15、)*TempP;endt=1:le ngth(Z);figure;plot(t,Z,b,t,X(t),r);title(使用 Kalman 对电力 系统负荷数据进行预测);xlabel(时间点数);ylabel(电力系统负荷 );axis tight;legend(负荷真实值,Kalman 预 测值)；figure;subplot(2,1,1);t=le ngth(x1):le ngth(x);plot(t,x(t),b,t,X(t),r);title(使用 Kalman 对电 力系统负荷数据进行预测);xlabel(时间点数 );ylabel(电力系统负荷);axis tight;legen
16、d(负 荷真实值,Kalman预测值);set(gca,XTick,le ngth(x1):2:le ngth(x); subplot(2,1,2);error=Z-X;plot(t,error(t);title(预测值与真实值之误差 );xlabel(时间点数);set(gca,XTick,le ngth(x1):2:le ngth(x);ylabel(5月5日预测值与真实值误差);axis tight;第六题：%小波变换clc;clear;close all;f=50;%信号频率oumiga=2*pi*f;N_sample=2048;%总采样点数Fs=1000;%采样频率t=0:1/Fs:
17、1;Tao=0.03;A=1;%信号幅度x =20*exp(-t/Tao)+20*si n( oumiga*t+pi/3)+12*si n(2*oumiga*t+pi/4)+10*si n( 3*oumiga*t+pi/6)+6*s in (4*oumiga*t+pi/8)+5*si n(5*oumiga*t +pi/5); % 信号函数表达式figure;plot(t,x);title(原始信号);xlabel(时间 t/s,FontSize,14);ylabel(幅值,FontSize,14);%原信号函数wave name=cmor3-3;totalscal=256;Fc=centfrq
18、(wavename); % 小波中心频率 c=2*Fc*totalscal;scals=c./(1:totalscal);f=scal2frq(scals,wave name,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(0 1 0 300);%短时傅里叶变换S,F,T,P=spectrogram(x,256,250,256,Fs); figure;surf(T,F,10*log10(P),edgecolor, non e
