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

1、考虑两个谐波信号x(t)和y(t)，其中x(t)?Acos(wct??),y(t)?Bcos(wct)式中A和wc为正的常数,?为均匀分布的随机变量，其概率密度函数为

?1?,0???2?, f(?)??2??0,其他?而B是一个具有零均值和单位方差的标准高斯随机变量，即其分布函数为

1fB(b)?exp(?b2/2),???b??

2?（1）求x(t)的均值ux(t)、方差?2x(t)、自相关函数Rx(?)和自协方差函数cx(?)。

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

解： （1）

x(t)的均值ux(t)为：

1ux?E(x(t))?E(Acos(wct??))?2??2?01Acos(wct??)d??Asin(wct??) ?0

2?02?方差?2x(t)为：

A2A2A2A2 ??E?(x(t))?E(Acos(wct??))?E(1?cos(2wct?2?))??E(cos(2wct?2?))?2222自相关函数Rx(?)为：

2x222Rx(?)?E(Acos(wct??)Acos(wct+wc???))?A2E(cos(wct??)cos(wct+wc???)) A2A2A2A2?E(cos(2wct+2wc??2?)?cos(wc?))?cos(wc?)?E(cos(2wct+2wc??2?))?cos(wc?)2222自协方差函数cx(?)为：

A2cx(?)?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(wct??)Bcos(wct?wc?))

由题意知道?与B为相互统计独立的随机变量，所以有

Rxy(?)?E(Acos(wct??)E(Bcos(wct?wc?))?AE(cos(wct??)E(B)cos(wct?wc?))?A?0?0?cos(wct?wc?)?0

互协方差函数cxy(?)

cxy(?)?Rxy(?)?0

2.接收信号由下式给出：yi?Acos(?ci??)??i,i?1,2,...,N，式中?i~N(0,1)即?i是零均

?c为载波角频率，值和单位方差的高斯噪声，而?是未知的相位。假设?1,?2,...?N相互独立，求未知相位的最大似然估计?ML。

解：由于?1,?2,...?N相互独立，所以y1,..yN也相互独立并且服从高斯分布，可以得到y1,..yN与?的联合概率密度函数分布

f(y1,..yN|?)?1(2?)N?^e?[yi?Acos(?ci??)]22i?1N

由此，可以得到似然函数

N1NL??ln(2?)??[yi?Acos(?ci??)]

22i?12该似然函数对?求偏导，并令该偏导函数为零，即可得到如下公式：

?LN??[yi?Acos(?ci??)]sin(?ci??)?0 ??i?1因此，最大似然估计?ML为上述函数的零点值。 则

^?Acos(?ci??ML)sin(?ci??ML)??yisin(?ci??ML)

i?1i?1N^^N^该函数为非线性方程，不容易求解，若忽略双倍频率2?c，则可简化到如下式子：

?yisin(?ci??)?0

i?1N^根据三角公式分解得到如下式子：

sin?ML?yicos(?ci)??cos?ML?yisin(?ci)

i?1i?1^N^N由此，可以得到如下公式

tan?ML??^i?1N?ysin?iicicN

?ycos(?i)i?1所以相位的最大似然估计如下：

?ML??arctan(^?ysin(?i)icN?ycos(?i)ici?1i?1N)

3.离散时间的二阶AR过程由差分方程x(n)?a1x(n?1)?a2x(n?2)?w(n)描述，式中w(n)2是一零均值、方差为?w的白噪声。证明x(n)的功率谱为

2?w21?a12?a2?2a1(1?a2)cos(2?f)?2a2cos(4?f)Px(f)?

证明：

由AR过程的功率谱公式知

Px(f)?2???j4?f21?a1e?j2?f?a2e

其中

1?a1e?j2?f?a2ej4?f2?(1?a1e?j2?f?a2ej4?f)(1?a1e?j2?f?a2ej4?f)2?1?a12?a2?a2(e?j4?f?e4?f)?a1a2ej2?f?a1a2e?j2?f?a1(e?j2?f?ej2?f) 2?1?a12?a2?2a1(1?a2)cos(2?f)?2a2cos(4?f)将其带入第一个公式可得：

Px(f)?2??21?a12?a2?2a1(1?a2)cos(2?f)?2a2cos(4?f)

x?t??sin(2?100t)?1.5sin(2?300t)?A?t?sin(2?200t)?dn?t??n?t?，4、信号的函数表达式为：

其中，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)之间总有如下输入输出关系：

x(n)???akx(n?k)??bku(n?k)

k?1k?0pqx(n)??h(k)u(n?k)

k?0?对以上两个式子两边分别取Z变换，并假定b0=1，可得

H(z)?pq?kB(z) A(z)?其中A(z)?1??akz,B(z)?1??bkz?k,H(z)??h(k)z?k。

k?1k?1k?0为了保证H(z)是稳定的最小相位系统，A(z)和B(z)的零点都应该在单位圆内。假定u(n)是一个方差为?2的白噪声序列，由随机信号通过线性系统的理论可知，输出序列X(n)的功率谱为：

Px(ejw)?2?2B*(ejw)B(ejw)A*(ejw)A(ejw)??2B(ejw)A(e)jw2

ARMA阶数确定：

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

ARMA法谱估计结果：

3、Burg最大熵法

Burg算法的具体实现步骤：

步骤1 计算预测误差功率的初始值和前、后向预测误差的初始值，并令m = 1。

1P0?N2

f0(n)?g0(n)?x(n)

n?1?x(n)N步骤2 求反射系数

?Km?Nn?m?1?fNm?1?(n)gm?1(n?1)122[fm?1(n)?gm?1(n?1)]?2n?m?1

步骤3 计算前向预测滤波器系数

?am(i)?am?1(i)?Kmam,...,m?1 ?1(m?i),i?1am(m)?Km

步骤4 计算预测误差功率

Pm?(1?Km)Pm?12

步骤5计算滤波器输出

fm(n)?fm?1(n)?Kmgm?1(n?1)

?gm(n)?Kmfm?1(n)?gm?1(n?1)

步骤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、预估计

?=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)=C(k)×H(k)'×[H(k)×C(k)×H(k)'+R(k)]-1

R(k)=N(k)×N(k)'

4、更新估计

???=X(k)+K(k)×[Y(k)-H(k)×X(k)] X(k)5、计算更新后估计协方差矩阵

??= [I-K(k)×H(k)]×C(k)×[I-K(k)×H(k)]'+K(k)×R(k)×K(k)' C(k)?X(k+1) =X(k) ?C(k+1) =C(k)

6、重复以上步骤

最终可以获得如下结果：

本题将表中的作为观测数据，图中横坐标为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 2ωt+45° +10sin 3ωt+30° +6sin 4ωt+22.5° +5sin 5ωt+36° 式中，ω=2πf0，τ=30ms，f0=50Hz分别利用小波变换、短时傅里叶变换和维格纳威利分布分析故障电流信号的时频特性。 解：

（1）小波变换：

连续小波变换的定义：

CWTf(u,s)??f(t),?u,s(t)???????f(t)1s?*(t?u)dt s计算连续时间小波变换的4个步骤：

1、选取一个小波，然后将其和待分析信号从起点开始的一部分进行相乘积分。 2、计算相关系数c。

3、将小波向右移，重复1和2的步骤直到分析完整个信号。

4、将小波进行尺度伸缩后再重复1，2，3步骤，直至完成所有尺度的分析。 （2）短时傅里叶变换

短时傅里叶变换定义如下：

STFTf(u,?)??f(t),gu,?(t)???????f(t)g(t?u)e?i?tdt

??STFTf(u,?)??????f(t)g(t?u)e?i?tdt?12?????(?)g??(???)eiu(???)d? f（3）维格纳威利分布变换

维格纳威利分布定义如下：

???????WDx(t,?)??x?t??x*?t??e-j??d????2??2?

在MATLAB中没有维格纳威利分布变换的相关函数，需要安装一个MATLAB版本的时频

分析工具箱。调用里面的函数即可。小波变换和短时傅里叶变换MATLAB均自带了相关的函数。程序见附录。代码运行结果结果如下：

7.假定一电力系统谐波与间谐波信号的函数表达式如下：

y?n??0.001cos(2??10n??1)?cos(2??50n??2)?0.1cos(2??150n??3)?0.002cos?2??50n??4????n?

其中，采样频率为1024Hz，相位?1??4为独立的均匀分布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*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; end end

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=[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); for i=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('电力系统负荷');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);

本文来源：https://www.bwwdw.com/article/gljw.html