整理实验二典型时间序列的功率谱估计.docx
- 文档编号:8551479
- 上传时间:2023-01-31
- 格式:DOCX
- 页数:16
- 大小:329.75KB
整理实验二典型时间序列的功率谱估计.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、体会与收获
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 整理 实验 典型 时间 序列 功率 估计
![提示](https://static.bdocx.com/images/bang_tan.gif)