整理实验二典型时间序列的功率谱估计.docx

整理实验二典型时间序列的功率谱估计.docx

《整理实验二典型时间序列的功率谱估计.docx》由会员分享，可在线阅读，更多相关《整理实验二典型时间序列的功率谱估计.docx（16页珍藏版）》请在冰豆网上搜索。

整理实验二典型时间序列的功率谱估计

实验二.典型时间序列的功率谱估计

一、实验内容与目标：

了解有限长数据对谱估计的影响，重点研究周期图法和改进型周期图法的谱估计方法，并分析噪声对随机信号谱估计结果的影响。

使学生了解随机信号的功率谱分析与主要估计方法。

2、实验任务

（1）对AR

（2）模型产生的序列进行分析，并估计其数字特征。

理论知识

在《随机信号分析》课程的第五章时间序列模型中，我们对AR、MA及ARMA模型进行了分析。

对于AR

（2）模型：

其解为：

其中

和

是由初始条件确定的待定系数，

；

而根据其自相关函数，有：

；

功率谱为：

。

已知条件

设有AR

（2）模型为：

即

是高斯白噪声，均值为0，方差为4；

任务一

产生指定AR（P）模型的典型时间序列X（n）；分别画出X（n）的500、2000点和10000个观测点的波形，并估计他们的均值与方差；试根据上述理论知识推导并计算出该模型理论的均值和方差。

程序及结果

程序：

b=1;

a=[10.90.1];

noise=normrnd（0,2,1,500）;

x=filter（b,a,noise）;

subplot（211）;

plot（x）;

title（'AR

（2）随机序列500点'）;

m=mean（x）;

sigma2=var（x）;

m

sigma2

noise=normrnd（0,2,1,2000）;

x=filter（b,a,noise）;

subplot（211）;

plot（x）;

title（'AR

（2）随机序列2000点'）;

m=mean（x）;

sigma2=var（x）;

m

sigma2

noise=normrnd（0,2,1,10000）;

x=filter（b,a,noise）;

subplot（211）;

plot（x）;

title（'AR

（2）随机序列10000点'）;

m=mean（x）;

sigma2=var（x）;

m

sigma2

m=0.0404sigma2=14.3090

m=-0.0486sigma2=14.5322

m=0.0058sigma2=15.0434

理论的均值和方差

均值为0，因为是一个线性系统，w（n）为白噪声，将其看成输入，输出x（n）也为0，所以方差就等于R（0）m（x）=0δ2=r（0）=15.77

分析

随着点数的增加（500,2000,10000），波形越来越密集，随着点数的增加波形越能体现出时间序列模型的特点。

任务二

求出该AR

（2）模型理论功率谱，在

上等距选取K=500个点，画出其理论功率谱曲线；

程序及结果

程序：

fs=1000;%采样频率

w=0:

pi/2000:

pi;

G=4*（abs（1./（1+0.9*exp（-j*w）+0.1*exp（-2*j*w）））.^2）;%AR模型系统函数

G=G/max（G）;%归一化

f=w*fs/（2*pi）;

plot（f,G,'r'）

title（'理论功率谱密度曲线'）

xlabel（'f'）

ylabel（'幅值'）

任务三

估计500、2000点和10000个观测点的典型时间序列X（n）的自相关函数与功率谱，并与理论的功率谱曲线比较；

程序、结果及分析

程序：

（功率谱）

b=1;

a=[10.90.1];

noise=normrnd（0,4,1,500）;

x=filter（b,a,noise）;

window=hamming（20）;%采用hanmming窗，长度为20

noverlap=10;

Nfft=512;

fs=1000;

[Px,f]=pwelch（x,window,noverlap,Nfft,fs,'onesided'）;%估计功率谱密度

f=[-fliplr（f）f（1:

end）];%对称频率（反转）

Px=[fliplr（Px）Px（1:

end）];%对称谱

Px/max（Px）

subplot（221）;

plot（f,Px,'b'）;

holdon;

fs=1000;

w=-pi:

1/fs:

pi;

G=4*（abs（1./（1+0.9*exp（-j*w）+0.1*exp（-2*j*w）））.^2）;

G=G/max（G）;

f=w*fs/（2*pi）;

subplot（221）;

plot（f,G,'--r'）

legend（'实际功率谱曲线','理论功率谱曲线'）

title（'实际功率谱与理论功率谱曲线比较图（500点）'）

b=1;

a=[10.90.1];

noise=normrnd（0,4,1,2000）;

x=filter（b,a,noise）;

window=hamming（20）;%采用hanmming窗，长度为20

noverlap=10;

Nfft=512;

fs=1000;

[Px,f]=pwelch（x,window,noverlap,Nfft,fs,'onesided'）;%估计功率谱密度

f=[-fliplr（f）f（1:

end）];%对称频率（反转）

Px=[fliplr（Px）Px（1:

end）];%对称谱

Px/max（Px）;

subplot（222）;

plot（f,Px,'b'）;

holdon;

fs=1000;

w=-pi:

1/fs:

pi;

G=4*（abs（1./（1+0.9*exp（-j*w）+0.1*exp（-2*j*w）））.^2）;

G=G/max（G）;

f=w*fs/（2*pi）;

subplot（222）;

plot（f,G,'--r'）;

legend（'实际功率谱曲线','理论功率谱曲线'）

title（'实际功率谱与理论功率谱曲线比较图（2000点）'）

b=1;

a=[10.90.1];

noise=normrnd（0,4,1,10000）;

x=filter（b,a,noise）;

window=hamming（20）;%采用hanmming窗，长度为20

noverlap=10;

Nfft=512;

fs=1000;

[Px,f]=pwelch（x,window,noverlap,Nfft,fs,'onesided'）;%估计功率谱密度

f=[-fliplr（f）f（1:

end）];%对称频率（反转）

Px=[fliplr（Px）Px（1:

end）];%对称谱

Px/max（Px）

subplot（223）;

plot（f,Px,'b'）;

holdon;

fs=1000;

w=-pi:

1/fs:

pi;

G=4*（abs（1./（1+0.9*exp（-j*w）+0.1*exp（-2*j*w）））.^2）;

G=G/max（G）;

f=w*fs/（2*pi）;

subplot（223）;

plot（f,G,'--r'）;

legend（'实际功率谱曲线','理论功率谱曲线'）

title（'实际功率谱与理论功率谱曲线比较图（10000点）'）

由图可知：

（1）实际功率谱曲线相对于理论功率谱曲线有点误差，但其包络还是比较一致的。

（2）有限长的数据和噪声都会影响谱估计结果，数据越长，功率谱估计越准确

程序：

（自相关）

b=1;

a=[10.90.1];

noise=normrnd（0,4,1,500）;

x=filter（b,a,noise）;

R=xcorr（x,'coeff'）;

subplot（211）;

plot（R）

title（'自相关函数估计（500点）'）

noise=normrnd（0,4,1,2000）;

x=filter（b,a,noise）;

R=xcorr（x,'coeff'）;

subplot（212）;

plot（R）

title（'自相关函数估计（2000点）'）

noise=normrnd（0,4,1,10000）;

x=filter（b,a,noise）;

R=xcorr（x,'coeff'）;

subplot（213）;

plot（R）

title（'自相关函数估计（10000点）'）

（2）不同信噪比下的功率谱估计与比较。

任务

对于给定的时间序列Y（n），设计2种不同白噪声W（n），分别画出10000点的波形Y（n）+W（n），估计功率谱，比较和分析估计结果；给定时间序列

该时间序列信号离散化采样频率为1KHz。

程序：

Fs=1000;

N=1024;Nfft=10000;

n=0:

N-1;

%第一种白噪声

yn=sin（0.3*pi*n）+cos（0.6*pi*n）+0.2*randn（1,N）;

Px1=10*log10（abs（fft（yn,Nfft）.^2）/N）;

f=（0:

length（Px1）-1）*Fs/length（Px1）;

subplot（211）;plot（f,Px1）;

xlabel（'频率/Hz'）;ylabel（'功率谱/dB'）;

title（'第一种白噪声功率谱'）;

%第二种白噪声

yn=sin（0.3*pi*n）+cos（0.6*pi*n）+0.8*randn（1,N）;

Px2=10*log10（abs（fft（yn,Nfft）.^2）/N）;

f=（0:

length（Px2）-1）*Fs/length（Px2）;

subplot（212）;plot（f,Px2）;

xlabel（'频率/Hz'）;ylabel（'功率谱/dB'）;

title（'第二种白噪声功率谱'）;

（3）设计2种不同特性的线性滤波器，考察分布为N（0,4）的高斯白噪声通过不同线性滤波器后的波形变化情况，画出两种滤波器下输出前后的波形及系统的功率谱估计，并对波形变化情况加以讨论。

线性滤波器一

程序及结果

程序：

N=1000;

x=normrnd（0,1,N,1）;

b=[1];

a=[1,-0.1,-0.8];

（1）可能造成重大环境影响的建设项目，编制环境影响报告书，对产生的环境影响应进行全面评价；y=filter（b,a,x）;

（三）环境价值的定义Psd=abs（（fft（y）））.^2/N;

（二）环境影响经济损益分析的步骤subplot（2,2,1）;

plot（x）;

title（'经过滤波器前的x的波形'）;

2.间接市场评估法xlabel（'n'）;

ylabel（'x'）;

subplot（2,2,2）;

plot（y）;

（二）规划环境影响评价的技术依据和基本内容title（'经过滤波器后的y的波形'）;

四、安全预评价xlabel（'n'）;

ylabel（'y'）;

3.意愿调查评估法subplot（2,2,3）;

plot（Psd）;

[答疑编号502334050102]title（'系统功率谱'）;

规划编制单位对规划环境影响进行跟踪评价，应当采取调查问卷、现场走访、座谈会等形式征求有关单位、专家和公众的意见。

xlabel（'n'）;

ylabel（'Psd'）;

分析

用[r,p,k]=residue（[1,0],[1,a1,a2]）函数求出H（z）的极点求的其极点p1=0.9458p2=-0.8458，可得极点一个为正，一个为负，所以功率谱在w=0和π出现峰值，w=0时极点0.9458离单位圆距离比w=π时极点-0.8458离单位圆的距离小，则此时正极点作用大，w=0时|H（ejw）|最大。

即|H（ejw）|在低频时分量较大，可以看作其具有低通性质。

线性滤波器二

程序及结果

N=1000;

x=normrnd（0,1,N,1）;

b=[1];

a=[1,0.1,-0.8];

y=filter（b,a,x）;

Psd=abs（（fft（y）））.^2/N;

subplot（2,2,1）;

plot（x）;

title（'经过滤波器前的x的波形'）;

xlabel（'n'）;

ylabel（'x'）;

subplot（2,2,2）;

plot（y）;

title（'经过滤波器后的y的波形'）;

xlabel（'n'）;

ylabel（'y'）;

subplot（2,2,3）;

plot（Psd）;

title（'系统功率谱'）;

xlabel（'n'）;

ylabel（'Psd'）;

分析

用[r,p,k]=residue（[1,0],[1,a1,a2]）函数求出H（z）的极点求的其极点p1=-0.9458p2=0.8458，可得极点一个为正，一个为负，所以功率谱在w=0和π出现峰值，w=0时极点-0.9458离单位圆距离比w=π时极点0.8458离单位圆的距离小，则此时负极点作用大，w=π时|H（ejw）|最大。

即|H（ejw）|在高频时分量较大，可以看作其具有高通性质。

3、体会与收获