整理实验二典型时间序列的功率谱估计.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、体会与收获
升级会员 