资源描述
1、考虑两个谐波信号 和 ，其中 , 式中 和()xtyt()cos()xtAwt(cos()ytBwtA为正的常数, 为均匀分布的随机变量，其概率密度函数为cw,1,02()f其 他而 是一个具有零均值和单位方差的标准高斯随机变量，即其分布函数为B 21()exp(/),Bfbb（1）求 的均值 、方差 、自相关函数 和自协方差函数 。()xt()xut2xtxR()xc（2）若 与 为相互统计独立的随机变量，求 和 的互相关函数 与互协方差 ()t yR函数 。xyc解：（1）的均值 为：()t()xut 220 011(cos)cos()sin() 2cEAwtAwtdAwt方差 为：2()xt2 22(cs)(1cs)(os2)c c cAtEtEt自相关函数 为：)xR22 2 2()os(o(+)(o)s(+)c+2)cscsc2cos()xc cEAwttwAwttww 自协方差函数 为：()x 2()os()xxcR（2） 的均值为：()yt,所以()()cos()()0yBBcutEtwtEBt()=0EB由互相关函数的定义可知： osxyccRAw由题意知道 与 为相互统计独立的随机变量，所以有()cos()(cos)()(s)00xy ccccRAwtttt 互协方差函数 xy ()0xyxyR2.接收信号由下式给出： ，式中 即 是零cos,12,.i iAiN(0,1)ii均值和单位方差的高斯噪声， 为载波角频率，而 是未知的相位。假设 相互2.N独立，求未知相位的最大似然估计 。ML解：由于 相互独立，所以 也相互独立并且服从高斯分布，可以得到12,.N1,.Ny与 的联合概率密度函数分布1,.Ny 21cos()1(,.|)(2)iyAiNNfye由此，可以得到似然函数21ln(2)cos()NiiLyAi该似然函数对 求偏导，并令该偏导函数为零，即可得到如下公式：1cos()sin()0Ni ciy 因此，最大似然估计 为上述函数的零点值。ML则 1 1cos()sin()sin()N NMLMLMLc ci iAi y 该函数为非线性方程，不容易求解，若忽略双倍频率 ，则可简化到如下式子：2c1si()0Nciy根据三角公式分解得到如下式子： 1 1sinos()ssin()NNMLicMLciiy由此，可以得到如下公式 1sitanco()NiMLiy所以相位的最大似然估计如下： 1sin()arct(oNciMLiy3.离散时间的二阶 AR 过程由差分方程 描述，式中12()()(xnaxnw是一零均值、方差为 的白噪声。证明 的功率谱为()wn2w22112()()cos()cos(4)wxPfaafaf证明：由 AR 过程的功率谱公式知 2241()xjfjfPfae其中 22424241 114 2221122()() (coscosjfjf jfjfjfjfjffjfjfjfjfae aeeaee将其带入第一个公式可得： 21122()()cs()cos(4)xPfaafaf4、信号的函数表达式为：，其中， 为一随时间变 sin(210).5sin(230)sin(20)xtttAtdntAt化的随机过程， 为经过 390-410Hz 带通滤波器后的高斯白噪声， 为高斯白噪声，d n采样频率为 1kHz，采样时间为 2.048s。分别利用周期图谱、 ARMA、Burg 最大熵方法估计信号功率谱，其中 ARMA 方法需要讨论定阶的问题。解：由题意知采样点数一共为：10002.048=2048 个数据点。 为一随时间变化的At随机过程，由于随机过程有很多类型，如维纳过程、正态随机过程，本文采用了均值为0，方差为 1 的正态随机过程来作为演示，来代替 ，高斯白噪声采用强度为 2 的高斯At白噪声代替 ，其带通滤波后为 。其中滤波器采用的是契比雪夫数字滤波器。nt dnt可得到 x(t)如下图所示：0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2-6-4-20246原始输入信号1、周期图法matlab 中的周期图功率谱法原理是通过计算采样信号的 FFT，获得离散点的幅度，再根据幅度与功率之间的关系，转换为离散点的功率，再通过坐标变换将离散点的功率图转换为连续功率谱密度。Step1:计算采样信号 x(n)的 DFT，使用 FFT 方法来计算。如果此处将复频率处的幅度对称到物理实际频率，得到的就是单边谱，否则就是双边谱Step2:根据正余弦信号功率与幅度的关系以及直流功率与幅度的关系，将幅度转换为离散功率谱。Step3:对横纵坐标进行转换，横坐标乘以频率分辨率转换为实际连续物理频率，纵坐标除以频率分辨率转换为功率谱密度。调用 MATLAB 中自带的 matlab 中Pxx,f=periodogram(x,window,nfft,fs)函数可得计算结果如下：0 50 100 150 200 250 300 350 400 450 500f/Hz00.511.52功率/db周期图法求功率谱2、ARMA 方法参数模型估计的思想是:假定研究的过程 X(n)是一个输入序列 u(n)激励一个线性系统 H(z)的输出。有已知的 X(n)，或其自相关函数来估计 H(z)的参数。由 H(z)的参数来估计 X(n)的功率谱。不论 X(n)是确定性信号还是随机信号，u(n)与 X(n)之间总有如下输入输出关系：10()()()pqkkxnaxnbun0kh对以上两个式子两边分别取 Z 变换，并假定 b0=1，可得()BzHA其中 , , 。1()pkkAzaz1()qkkBz0()()kkh为了保证 H(z)是稳定的最小相位系统，A(z) 和 B(z)的零点都应该在单位圆内。假定 u(n)是一个方差为 的白噪声序列，由随机信号通过线性系统的理论可知，输出序列 X(n)的功2率谱为： 222*()()() jwjwjjwx jBeBePeAAARMA 阶数确定：本题目采用 AIC 准则确定 ARMA 的阶数。分别计算 p、q 从 1 到 20 阶数的计算出AIC（p,q）,如下图所示，当横坐标大概为 230 左右时，AIC(p,q)取得最小，将此时的 p,q作为带入到模型即可。0 50 100 150 200 250 300 350 4000.70.80.911.11.21.3AIC(p,q)ARMA 法谱估计结果：0 50 100 150 200 250 300 350 400 450 500f/Hz-5051015202530振幅/dBARMA法 (AIC准则 )3、Burg 最大熵法Burg 算法的具体实现步骤：步骤 1 计算预测误差功率的初始值和前、后向预测误差的初始值，并令 m = 1。210)(NnxP)(gf步骤 2 求反射系数NmnmngfK1212)()(步骤 3 计算前向预测滤波器系数),()(11iaiiam,.m步骤 4 计算预测误差功率 12)(PK步骤 5 计算滤波器输出 )()(11ngfnfmmg步骤 6 令 mm+1，并重复步骤 2 至步骤 5，直到预测误差功率 Pm 不再明显减小。最后，再利用 Levinson 递推关系式估计 AR 参数，继而得到功率谱估计。Burg 最大熵法谱估计结果如下图：0 50 100 150 200 250 300 350 400 450 500f/fs00.511.5功率谱/dBBurg法谱估计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)Q(k)=U(k)U(k)3、计算卡尔曼增益矩阵K(k)= H(k)H(k) H(k)+R(k)-1()C(k)R(k)=N(k)N(k)4、更新估计= +K(k)Y(k)-H(k) X(k) X()5、计算更新后估计协方差矩阵= I-K(k)H(k) I-K(k)H(k)+K(k)R(k)K(k)C()C()X(k+1) = (k)C(k+1) =6、重复以上步骤最终可以获得如下结果：20 40 60 80 100 120 140 160 180时间点数500100015002000250030003500电力系统负荷使用 Kalman对电力系统负荷数据进行预测负荷真实值Kalman预测值168170172174176178180182184186188190192时间点数250030003500电力系统负荷使用 Kalman对电力系统负荷数据进行预测负荷真实值Kalman预测值168170172174176178180182184186188190192时间点数-10001002003005月5日预测值与真实值误差预测值与真实值之误差本题将表中的作为观测数据，图中横坐标为 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(2t+45)+10sin(3t+30)+6sin(4t+22.5)+5sin(5t+36)式中， ， ， 分别利用小波变换、短时傅里叶变换和维格纳威利分=2f0 =30msf0=50Hz布分析故障电流信号的时频特性。解：（1）小波变换：连续小波变换的定义： *, 1(,)()()()f us tuCWTusftftds计算连续时间小波变换的 4 个步骤：1、选取一个小波，然后将其和待分析信号从起点开始的一部分进行相乘积分。 2、计算相关系数 c。3、将小波向右移，重复 1 和 2 的步骤直到分析完整个信号。4、将小波进行尺度伸缩后再重复 1，2，3 步骤，直至完成所有尺度的分析。（2）短时傅里叶变换短时傅里叶变换定义如下： ,(,)()()itf uSTFuftgftgued()1, 2it iuf tedf （3）维格纳威利分布变换维格纳威利分布定义如下： de),(WDj-*txtx在 MATLAB 中没有维格纳威利分布变换的相关函数，需要安装一个 MATLAB 版本的时频分析工具箱。调用里面的函数即可。小波变换和短时傅里叶变换 MATLAB 均自带了相关的函数。程序见附录。代码运行结果结果如下：0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1时间 t/s-20020406080幅值原始信号小波时频图0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1时间 t/s050100150200250300频率f/Hz102030405060700.2 0.3 0.4 0.5 0.6 0.7 0.8时间 /s050100150200250300350400450500频率/Hz短时傅里叶变换结果Wigner-Ville time-frequency distribution0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1时间 t/s050100150200250300频率f/Hz7.假定一电力系统谐波与间谐波信号的函数表达式如下： 12 3 4.1cos(2)cos(50).1cos(50).2cos50ynnnnn其中，采样频率为 ，相位 为独立的均匀分布 ； 为一噪声信号，04Hz14:,U信噪比取为 。分别采用三种现代信号处理方法进行谐波与间谐波频率提取与谱估计。dB解：本题目采用的频率提取的三种方法为小波变换、短时傅里叶变换和维格纳威利分布。采用周期图法、MUSIC 法、Burg 法进行谱估计。确定出谐波的频率为 50Hz 和150Hz。0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2时间 (t/秒 )-3-2-10123幅值原始信号小波时频图0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2时间 t/s100200300400500频率f/Hz0.511.522.533.50.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8时间 /s0100200300400500频率/Hz短时傅里叶变换结果Wigner-Ville time-frequency distribution0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2时间 t/s050100150200250300350400450500频率f/Hz0 50 100150200250300350400450500550600f/Hz00.10.20.30.40.5功率/db周期图法求功率谱0 50 100 150 200 250 300 350 400 450 500频率 (f/Hz)-40-200204060功率(dB)MUSIC方法0 50 100150200250300350400450500550600f/fs00.0050.010.015功率谱/dBBurg法谱估计附录代码：第四题：clc;clear;fs=1000;%采样频率T=2.048;%采样时间t=0:1/fs:T;A = normrnd(0,1,1,length(t);%方差为 1，均值为 0 的高斯分布N=wgn(1,length(t),2);%强度为 2 的高斯白噪声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(0 2.1 -7 7);% 周期图谱Pxx,f=periodogram(Sig,length(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 和 q for q=1:20%移动平均对应 ACF m=armax(z(1:length(t),p,q); AIC = aic(m); %armax(p,q)选择对应 FPE 最小，AIC 值最小模型 RecordAIC=RecordAIC;p q AIC; end end for k=1:size(RecordAIC,1) if RecordAIC(k,3)=min(RecordAIC(:,3) %选择 AIC 最小模型 pa_AIC=RecordAIC(k,1); qa_AIC=RecordAIC(k,2); break; endend mAIC=armax(z(1:length(t),pa_AIC,qa_AIC);Pxx2,f2=freqz(mAIC.c,mAIC.a,fs); P2=(abs(Pxx2).*1).2; 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 轴坐标名称% function y=bandp(x,f1,f3,fsl,fsh,rp,rs,Fs)%带通滤波%使用注意事项：通带或阻带的截止频率与采样率的选取范围是不能超过采样率的一半%即，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*f3/Fs;wsl=2*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(./负荷数据 .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);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)*TempP;endt=1:length(Z);figure;plot(t,Z,b,t,X(t),r);title(使用 Kalman 对电力系统负荷数据进行预测);xlabel(时间点数);ylabel(电力系统负荷);axis tight;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(电力系统负荷);axis tight;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 日预测值与真实值误差);axis tight;第六题：% 小波变换 clc;clear; close all;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(0 1 0 300);% 短时傅里叶变换S,F,T,P=spectrogram(x,256,250,256,Fs);figure;surf(T,F,10*log10(P),edgecolor,none); axis tight;view(0,90);xlabel(时间/s); ylabel( 频率/Hz);title(短时傅里叶变换结果);% Wigner-Ville time-frequency distribution.X=hilbert(x);tfr,t,f=tfrwv(X);figure;contour(t/Fs,f*Fs,abs(tfr);xlabel(时间 t/s);ylabel(频率 f/Hz);title(Wigner-Ville time-frequency distribution);axis(0 1 0 300)%第七题：clc;clear;close all;% 参数设置Fs = 1024; %采样频率n = 0:1/Fs:2.01;%采样时间N = length(n); % 采样点W1=0.001*cos(2*pi*n*10+unifrnd(-pi,pi)+cos(2*pi*50*n+unifrnd(-pi,pi)+0.1*cos(2*pi*n*150+unifrnd(-pi,pi)+0.002*cos(2*pi*n*50+unifrnd(-pi,pi);% 原始信号x1=awgn(W1,20); %加入噪声%原信号输出figure;plot(n,x1);xlabel(时间(t/秒),FontSize,10);ylabel(幅值,FontSize,10); axis(0 2.05 -3 3);title(原始信号 );% 小波变换wavename=cmor3-3;totalscal=256;Fc=centfrq(wavename); %小波中心频率c=2*Fc*totalscal;scals=c./(1:totalscal);f=scal2frq(scals,wavename,1/Fs); % 将尺度转换为频率coefs=cwt(x1,scals,wavename); % 求连续小波系数figure;imagesc(n,f,abs(coefs); colorbar;xlabel(时间 t/s,FontSize,14);ylabel(频率 f/Hz,FontSize,14);title(小波时频图,FontSize,16);% 短时傅里叶变换S,F,T,P=spectrogram(x1,256,250,256,Fs);figure;surf(T,F,10*log10(P),edgecolor,none); axis tight;view(0,90);xlabel(时间/s); ylabel( 频率/Hz);title(短时傅里叶变换结果);% 维格纳威利分布X=hilbert(x1);tfr,t,f=tfrwv(X);figure;contour(t/Fs,f*Fs,abs(tfr);xlabel(时间 t/s);ylabel(频率 f/Hz);title(Wigner-Villetime-frequency distribution);% 周期图谱估计Pxx,f=periodogram(x1,length(x1),Fs);%周期图法figure;plot(f,Pxx);title(周期图法求功率谱);xlabel(f/Hz); ylabel(功率/db);set(gca,XTick,0:50:600);% MUSIC 方法 谱估计nfft=1024;figure;pmusic(x1,7,1.1,nfft,Fs,32,16); grid on;xlabel(频率(f/Hz),FontSize,10);ylabel(功率(dB),FontSize,10);title(MUSIC 方法);% burg 法 谱估计Pxx,F = pburg(x1,60,length(x1),Fs);%burg法figure;plot(F,Pxx);title(Burg 法谱估计);xlabel(f/fs); %X 轴坐标名称ylabel(功率谱/dB); %Y 轴坐标名称set(gca,XTick,0:50:600);
展开阅读全文