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

1、考虑两个谐波信号()x t 和()y t ，其中()cos()c x t A w t φ=+,()cos()c y t B w t =式中A 和cw 为正的常数,φ为均匀分布的随机变量，其概率密度函数为1,02()20,f φπφπ⎧≤≤⎪=⎨⎪⎩其他,而B 是一个具有零均值和单位方差的标准高斯随机变量，即其分布函数为 （1）求()x t 的均值()x u t 、方差2()x t σ、自相关函数()x R τ和自协方差函数()x c τ。

（2）若φ与B 为相互统计独立的随机变量，求()x t 和()x t 的互相关函数()xy R τ与互协方差函数()xy c τ。

解： （1）()x t 的均值()x u t 为： 方差2()x t σ为：自相关函数()x R τ为： 自协方差函数()x c τ为： （2）()y t 的均值为：()(())(cos())()cos()0y B B c c u t E y t E B w t E B w t ====,所以()=0E B由互相关函数的定义可知：由题意知道φ与B 为相互统计独立的随机变量，所以有 互协方差函数()xy c τ2.接收信号由下式给出：cos(),1,2,...,i c i y A i i N ωθω=++=，式中~(0,1)i N ω即i ω是零均值和单位方差的高斯噪声，c ω为载波角频率，而θ是未知的相位。

假设12,,...N ωωω相互独立，求未知相位的最大似然估计^ML θ。

解：由于12,,...N ωωω相互独立，所以1,..N y y 也相互独立并且服从高斯分布，可以得到1,..N y y 与θ的联合概率密度函数分布 由此，可以得到似然函数该似然函数对θ求偏导，并令该偏导函数为零，即可得到如下公式：因此，最大似然估计^ML θ 为上述函数的零点值。

则该函数为非线性方程，不容易求解，若忽略双倍频率2c ω ，则可简化到如下式子： 根据三角公式分解得到如下式子： 由此，可以得到如下公式所以相位的最大似然估计如下：3.离散时间的二阶AR 过程由差分方程12()(1)(2)()x n a x n a x n w n =-+-+ 描述，式中()w n 是一零均值、方差为2w σ 的白噪声。

证明()x n 的功率谱为证明：由AR 过程的功率谱公式知 其中将其带入第一个公式可得：4、信号的函数表达式为：()()()()sin(2100) 1.5sin(2300)sin(2200)x t t t A t t dn t n t πππ=++++，其中，()A t 为一随时间变化的随机过程，()dn t 为经过390-410Hz 带通滤波器后的高斯白噪声，()n t 为高斯白噪声，采样频率为1kHz ，采样时间为2.048s 。

分别利用周期图谱、ARMA 、Burg 最大熵方法估计信号功率谱，其中ARMA 方法需要讨论定阶的问题。

解：由题意知采样点数一共为：1000×2.048=2048个数据点。

()A t 为一随时间变化的随机过程，由于随机过程有很多类型，如维纳过程、正态随机过程，本文采用了均值为0，方差为1的正态随机过程来作为演示，来代替()A t ，高斯白噪声采用强度为2的高斯白噪声代替()n t ，其带通滤波后为()dn t 。

其中滤波器采用的是契比雪夫数字滤波器。

可得到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)之间总有如下输入输出关系：对以上两个式子两边分别取Z 变换，并假定b 0=1，可得 其中1()1p kk k A z a z -==+∑,1()1q kk k B z b z -==+∑,0()()k k H z h k z ∞-==∑。

为了保证H(z)是稳定的最小相位系统，A(z)和B(z)的零点都应该在单位圆内。

假定u(n)是一个方差为2σ的白噪声序列，由随机信号通过线性系统的理论可知，输出序列X(n)的功率谱为：ARMA 阶数确定：本题目采用AIC 准则确定ARMA 的阶数。

分别计算p 、q 从1到20阶数的计算出AIC （p,q ）,如下图所示，当横坐标大概为230左右时，AIC(p,q)取得最小，将此时的p,q 作为带入到模型即可。

ARMA 法谱估计结果： 3、Burg 最大熵法Burg 算法的具体实现步骤：步骤1 计算预测误差功率的初始值和前、后向预测误差的初始值，并令m = 1。

步骤2 求反射系数步骤3 计算前向预测滤波器系数 步骤4 计算预测误差功率 步骤5计算滤波器输出步骤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、重复以上步骤最终可以获得如下结果：本题将表中的作为观测数据，图中横坐标为11时刻数据，22时刻的数据，一次类推，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）短时傅里叶变换 短时傅里叶变换定义如下： （3）维格纳威利分布变换 维格纳威利分布定义如下：在MATLAB 中没有维格纳威利分布变换的相关函数，需要安装一个MATLAB 版本的时频分析工具箱。

调用里面的函数即可。

小波变换和短时傅里叶变换MATLAB 均自带了相关的函数。

程序见附录。

代码运行结果结果如下：7.假定一电力系统谐波与间谐波信号的函数表达式如下：其中，采样频率为1024Hz ，相位14φφ为独立的均匀分布[],U ππ-+；()n ξ为一噪声信号，信噪比取为20dB 。

分别采用三种现代信号处理方法进行谐波与间谐波频率提取与谱估计。

解：本题目采用的频率提取的三种方法为小波变换、短时傅里叶变换和维格纳威利分布。

采用周期图法、MUSIC法、Burg法进行谱估计。

确定出谐波的频率为50Hz和150Hz。

附录代码：第四题：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*3 00.*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和qfor q=1:20%移动平均对应ACFm=armax(z(1:length(t)),[p,q]);AIC = aic(m); %armax(p,q)选择对应FPE最小，AIC值最小模型RecordAIC=[RecordAIC;p q AIC];endendfor k=1:size(RecordAIC,1)ifRecordAIC(k,3)==min(RecordAIC(:, 3)) %选择AIC最小模型pa_AIC=RecordAIC(k,1);qa_AIC=RecordAIC(k,2);break;endendmAIC=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轴坐标名称%%functiony=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','non e'); 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,p i))+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,25 6,Fs);figure;surf(T,F,10*log10(P),'edgecolor','non e'); 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);。