典型时间序列模型分析(可打印修改) (2)

1. 产生样本函数，并画出波形

2. 题目中的AR 过程相当于一个零均值正态白噪声通过线性系统后的输出，可以按照上面的方法进行描述。

clear all;

b=[1]; a=[1 0.3 0.5]; % 由描述的差分方程，得到系统传递函数

h=impz(b,a,20); % 得到系统的单位冲激函数，在20 点处已经可以认为值是0

randn('state',0);

w=normrnd(0,2,1,500); % 产生题设的白噪声随机序列，标准差为2

x=filter(b,a,w); % 通过线形系统，得到输出就是题目中要求的2 阶AR 过程

plot(x,'r');

ylabel('x(n)');

title('邹先雄——产生的AR 随机序列');

grid on;

得到的输出序列波形为：

2. 估计均值和方差

可以首先计算出理论输出的均值和方差，得到

x

m

，对于方差可以先求出理论自相

关输出，然后取零点的值。

并且，，带入有

在最大值处输出的功率，也就是方差，为

对实际数据进行估计，均值为mean(x)=-0.0703，而方差为var(x)=5.2795，两者合理论值吻合得比较好。

程序及运行结果图如下，其中y_mean表示均值，y_var表示方差。

3.画出理论的功率谱密度曲线

理论的功率谱为，

用下面的语句产生：

delta=2*pi/1000;

w_min=-pi;

w_max=pi;

Fs=1000;

w=w_min:delta:w_max; % 得到数字域上的频率取样点，范围是[-pi,pi]

Gx=4*(abs(1./(1+0.3*exp(-i*w)+0.5*exp(-2*i*w))).^2); % 计算出理论值

Gx=Gx/max(Gx); % 归一化处理

f=w*Fs/(2*pi); % 转化到模拟域上的频率

plot(f,Gx);

title('邹先雄——理论功率谱密度曲线');

grid on;

得到的图形为：

可以看出，这个系统是带通系统。

4. 估计自相关函数和功率谱密度

用实际数据估计自相关函数和功率谱的方法前面已经讨论过，在这里仅给出最后的仿真图形。

Mlag=20; % 定义最大自相关长度

Rx=xcorr(x,Mlag,'coeff');

m=-Mlag:Mlag;

stem(m,Rx,'r.');

title('邹先雄——自相关函数');

最终的值为

可以看出，它和上面的理论输出值吻合程度很好。实际的功率谱密度可以用类似于上面的方法进行估计，

window=hamming(20); % 采用hanmming 窗，长度为20

noverlap=10; % 重叠的点数

Nfft=512; % 做FFT 的点数

Fs=1000; % 采样频率，为1000Hz

b=[1]; a=[1 0.3 0.5]; % 由描述的差分方程，得到系统传递函数

h=impz(b,a,20); % 得到系统的单位冲激函数，在20 点处已经可以认为值是0

randn('state',0);

w=normrnd(0,2,1,500); % 产生题设的白噪声随机序列，标准差为2

x=filter(b,a,w); % 通过线形系统，得到输出就是题目中要求的2 阶AR 过程

[Px,f]=pwelch(x,window,noverlap,Nfft,Fs, 'onesided'); % 估计功率谱密度

f=[-fliplr(f) f(1:end)]; % 构造一个对称的频率，范围是[-Fs/2, Fs/2]

Py=[-fliplr(Px) Px(1:end)]; % 对称的功率谱

plot(f,10*log10(Py),'b');

title('邹先雄——实际的功率谱密度曲线');

估计出来的功率谱密度为，

将两幅图画在一起，可以看到拟合的情况比较好（两者相位刚好相反，但是基本波形相似）：

代码如下：

clear all;

delta=2*pi/1000;

w_min=-pi;

w_max=pi;

Fs=1000;

w=w_min:delta:w_max; % 得到数字域上的频率取样点，范围是[-pi,pi]

Gx=4*(abs(1./(1+0.3*exp(-i*w)+0.5*exp(-2*i*w))).^2); % 计算出理论值

Gx=Gx/max(Gx); % 归一化处理

f=w*Fs/(2*pi); % 转化到模拟域上的频率结束

plot(f,Gx,'r');

hold on;

title('邹先雄——理论和实际的功率谱密度曲线拟合');

window=hamming(20); % 采用hanmming 窗，长度为20

noverlap=10; % 重叠的点数

Nfft=512; % 做FFT 的点数

Fs=1000; % 采样频率，为1000Hz

b=[1]; a=[1 0.3 0.5]; % 由描述的差分方程，得到系统传递函数

h=impz(b,a,20); % 得到系统的单位冲激函数，在20 点处已经可以认为值是0

randn('state',0);

w=normrnd(0,2,1,500); % 产生题设的白噪声随机序列，标准差为2

x=filter(b,a,w); % 通过线形系统，得到输出就是题目中要求的2 阶AR 过程