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

实验二. 典型时间序列的功率谱估计 一、实验内容与目标： 了解有限长数据对谱估计的影响，重点研究周期图法和改进型周期图法的谱估计方法，并分析噪声对随机信号谱估计结果的影响。使学生了解随机信号的功率谱分析与主要估计方法。 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=[1 0.9 0.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.0404 sigma2=14.3090 m=-0.0486 sigma2=14.5322 m=0.0058 sigma2=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=[1 0.9 0.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'); hold on; 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=[1 0.9 0.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'); hold on; 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=[1 0.9 0.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'); hold on; 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=[1 0.9 0.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]; 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.9458 p2=-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.9458 p2=0.8458，可得极点一个为正，一个为负，所以功率谱在w=0和π出现峰值，w=0时极点-0.9458离单位圆距离比w=π时极点0.8458离单位圆的距离小，则此时负极点作用大，w=π时|H(ejw)|最大。即|H(ejw)|在高频时分量较大，可以看作其具有高通性质。 3、 体会与收获 _1234567893.unknown _1234567895.unknown _1234567897.unknown _1234567899.unknown _1234567900.unknown _1234567898.unknown _1234567896.unknown _1234567894.unknown _1234567891.unknown _1234567892.unknown _1234567890.unknown