随机信号与系统第五章习题部分答案

第五章 习题
5-1 设某信号为
1000||()t x t e -=
（1）试求x (t )的傅里叶变换X (j ω)，并绘制X (j ω)曲线；
（2）假设分别以采样频率为f s =5000Hz 和f s =1000Hz 对该信号进行采样，得到一组采样序列x k ，说明采样频率对序列x k 频率特性X (e j Ω)的影响。

解：（1）1000||62
2000
()()10j t t j t X j x t e dt e e dt ωωωω
∞
∞
----∞
-∞
=
==
+⎰
⎰. X (j ω)的曲线如下图所示：
（2）设采样周期为T ，则采样输出为()()()()k k k x x t t kT x kT t kT δδ∞∞
=-∞
=-∞=-=-∑∑.由
时域相乘等于频域卷积，有
1122()()*[()]()*[
()]22j k k X e X j t kT X j k
T T
π
π
δδππ∞
∞
Ω
=-∞
=-∞
=Ω-=ΩΩ-∑∑F 121212()()()2k k X j k d X j jk T T T T T
πππωδωωπ∞∞∞-∞=-∞=-∞=⋅=Ω--=Ω-∑∑⎰. 即序列x k 频率特性X (e j Ω)是原信号频谱X (j ω)以
2T
π
为周期进行延拓而成的，而采样频率1122s f T T
ππ=
=⋅，所以采样频率越高，序列x k 频率特性的各周期越分散，越不容易发生频谱混叠。

5-2 假设平稳随机过程x (t )和y (t )满足下列离散差分方程
11;k k k k k k k x ax e y ay x v ---=-=+
式中，|a|<1；e k ，v k ～N (0，σ 2)分布，且二者互不相关。

试求随机序列y k 的功率谱。

解：对1k k k x ax e --=进行离散时间傅里叶变换（DTFT ），且记DTFT(x k )=X (e j Ω)，DTFT(e k )=E (e j Ω)，则有
j j j ()(1)()X e ae E e ΩΩΩ--=
式中，Ω=ωT s ，称为数字频率（rad ），ω为实际频率（rad/s ），T s 为采样周期（s ）。

如果把上式看作以e k 为输入、以x k 为输出的线性滤波器，则有
j def j j j ()1
()()()1X e H e H E e ae
ΩΩ
ΩΩ
-===Ω- 式中，H (Ω)称为线性滤波器的频率传递函数。

根据线性系统的功率传递函数，若记输入序列e k 的功率谱为N (Ω)，则输出序列x k 的功率谱S x (Ω)可表示为
2
2
j 2()()|()|ˆ|1|x S N H ae
ΩσΩΩΩ-==
-
将AR(1)模型参数估计值a ˆ和噪声功率估计值σˆ2代入上式中相应的项，即可得到AR(1)序列x k 的AR 谱估计：
2j 2
ˆˆ()ˆ|1|
x
S ae
Ω
σΩ-=-
由11;
k k k k k k k x ax e y ay x v ---=-=+可知：11k k k k k y ay ax e v ----=+
对11k k k k k y ay ax e v ----=+进行离散时间傅里叶变换（DTFT ），且记DTFT(y k )=Y (e j Ω)，DTFT(v k )=V (e j Ω)，则有：()()()()()1j j j j j j Y e
ae ae
X e
E e V e Ω
-Ω
-Ω
Ω
Ω
Ω
--=+
j def ()
1
j ()()j j j ()()1Y e
G e
G X e V e ae
Ω
Ω
ΩΩΩ===Ω-+-
所以：
()()222
2
j 4j 2()()|()||()|
ˆˆ|1||1|y x S N G S V G ae
ae
ΩΩΩΩΩΩσσ--==Ω+Ω⎡⎤⎣⎦=
+
--
5-3 已知某一线性系统的单位脉冲响应函数为
,0()t e t h t -⎧≥=⎨
⎩0,
其它
假定输入x (t )是一零均值的高斯白噪声，其功率谱为S x (f )=N 0，试求该线性系统输出响应
y (t )的功率谱和协方差函数。

解：系统的频率响应函数为
(1)0
1
()()1j t
t j t
t j H j h t e
dt e e
dt e dt j ωωωωω
∞
∞
∞
----+-∞
====
+⎰⎰⎰
则 2
2
111
()()*()111H j H j H j j j ωωωωωω
==
⋅=+-+ 而输出响应y (t )的功率谱2
()()()y x S H j S ωωω=，所以该线性系统输出响应y (t )的功率谱为 022
1
()()11y x
N S S ωωωω=
=++. 可见y (t )的功率谱不是常数，所以输出响应y (t )不再是白噪声。

y (t )的自相关函数为
1
002
1
()[()]212
j y y N N R S e d e τ
ωττωωπ
ω∞
---∞==
=+⎰F
∴ 直流分量功率 2
0()lim 02
y y N R e ττμ-→∞=∞==
∴ 输出响应y (t )的协方差函数2
00()()022
y y y N N C R e e ττττμ--=-=-=.
5-4 请分别用Levinson 递推算法和Burg 算法估计信号
123()cos(2π)cos(2π)cos(2π)()x t f t f t f t e t =+++
的功率谱。

式中，f 1=150Hz ，f 2=200Hz ，f 3=210Hz ；e (t )是方差为0.1的白噪声过程。

解：MATLAB 程序如下：
f1=150;f2=200;f3=210; % 信号频率f1，f2，f3； Fs=1000; N=155;
% 采样频率Fs ，序列点数N ；
k=0:1/Fs:N/Fs;
x=cos(2*pi*f1*k)+cos(2*pi*f2*k)+cos(2*pi*f3*k)+randn(size(k)); % x k 序列
R=zeros(1,N+1);
% 相关函数初始化，R0=R(1)
% 估计 x k 的相关函数
for m=1:N+1 RXm=0;
for k=1:N+2-m
RX=x(k+m-1)*x(k);
RXm=RXm+RX;
end
R(m)=RXm/N;
end
a=zeros(N+1,N+1); sigma2=zeros(1,N+1); % 参数a k和估计量方差 k2
FPE=zeros(1,N+1); % 最终预测误差准则FPE
% 计算一阶AR模型的未知参数
sigma2(1)=R(1); a(1,1)=-R(2)/R(1);
sigma2(2)=(1-(abs(a(1,1)))^2)*sigma2(1);
FPE=sigma2(2)*(N+2)/N; % 一阶AR模型的最终预测误差FPE % Levinson递推算法
for k=2:N
RXk=0;
for m=1:k-1
aRX=a(k-1,m)*R(k-m+1);
RXk=RXk+aRX;
end
a(k,k)=-(R(k+1)+RXk)/sigma2(k);
for m=1:k-1
a(k,m)=a(k-1,m)+a(k,k)*a(k-1,k-m);
end
sigma2(k+1)=(1-(abs(a(k,k))^2))*sigma2(k);
FPE(k)=sigma2(k+1)*(N+k+1)/(N-k+1); % k阶AR模型的最终预测误差FPE
end
% 确定AR模型阶次
min=FPE(1);
for k=2:N
if FPE(k)<min
min=FPE(k);
p=k;
end
end
disp(‘输出模型阶次 p');
disp(p);
disp(‘输出AR(p)模型参数 a');
for k=1:p
disp(a(p,k));
end
disp(‘估计量协方差 sigma2');
disp(sigma2(p+1));
% AR谱估计
H=0; W=0:0.01:pi;
for k=1:p
H=H+a(p,k).*exp(-j*k*W);
end
Sx=sigma2(p+1)./((abs(1+H)).^2);
f=W*Fs/(2*pi); % 将角频率转化为频率
figure（1）；
plot(f,10*log10(Sx));grid % 打印对数功率谱
Burg算法程序：
f1=150;f2=200;f3=210; % 信号频率f1，f2，f3；
Fs=1000; N=155; % 采样频率Fs，序列点数N；
k=0:1/Fs:N/Fs;
x=cos(2*pi*f1*k)+cos(2*pi*f2*k)+cos(2*pi*f3*k)+randn(size(k)); % xk 序列
R=zeros(1,N+1); % 相关函数初始化，R0=R(1) % 估计 xk 的相关函数
for m=1:N+1
RXm=0;
for k=1:N+2-m
RX=x(k+m-1)*x(k);
RXm=RXm+RX;
end
R(m)=RXm/N;
end
a=zeros(N+1,N+1); sigma2=zeros(1,N+1); % 参数ak和估计量方差?k2
FPE=zeros(1,N+1); % 最终预测误差准则FPE
% 计算一阶AR模型的未知参数
sigma2(1)=R(1); a(1,1)=-R(2)/R(1);
sigma2(2)=(1-(abs(a(1,1)))^2)*sigma2(1);
FPE=sigma2(2)*(N+2)/N; % 一阶AR模型的最终预测误差FPE
% Levinson递推算法
for k=2:N
RXk=0;
for m=1:k-1
aRX=a(k-1,m)*R(k-m+1);
RXk=RXk+aRX;
end
a(k,k)=-(R(k+1)+RXk)/sigma2(k);
for m=1:k-1
a(k,m)=a(k-1,m)+a(k,k)*a(k-1,k-m);
end
sigma2(k+1)=(1-(abs(a(k,k))^2))*sigma2(k);
FPE(k)=sigma2(k+1)*(N+k+1)/(N-k+1); % k阶AR模型的最终预测误差FPE
end
% 确定AR模型阶次
min=FPE(1);
for k=2:N
if FPE(k)<min
min=FPE(k);
p=k;
end
end
order=p;
nfft=1024;
pburg(x,order,nfft,Fs);
disp('输出模型阶次 p');
disp(p);
disp('输出AR(p)模型参数 a');
for k=1:p
disp(a(p,k));
end
disp('估计量协方差 sigma2');
disp(sigma2(p+1));
% AR谱估计
H=0; W=0:0.01:pi;
for k=1:p
H=H+a(p,k).*exp(-j*k*W);
end
Sx=sigma2(p+1)./((abs(1+H)).^2);
f=W*Fs/(2*pi); % 将角频率转化为频率figure(1);
plot(f,10*log10(Sx));grid % 打印对数功率谱
5-5 分别考虑表p5-1和表p5-2所给出的时间序列值。

（1） 已知表p5-1所给出序列的数学模型为
121.200.55k k k k x x x e --=-+
式中，e k 是均值为零、方差为0.5的高斯白噪声。

试分别按DFT 和AR 谱估计方法，估计该时间序列的功率谱，并比较二者的异同之处。

（2）已知表p5-2所给出序列的数学模型为
10.72k k k x e e -=-
式中，e k 是均值为零、方差为0.7的高斯白噪声。

试分别按DFT 和AR 谱估计方法，估计该时间序列的功率谱，并比较二者的异同之处。

表p5-1（习题5-5）
表p5-2（习题5-5）
解：（1）DFT法估计该序列的功率谱程序如下：
x = [4.200 5.800 6.900 7.620 5.570 3.340 2.000 1.700 2.020 2.710 3.630 5.180 7.110 8.260 7.960 6.780 5.070 5.040 6.020 7.610 10.320];
Nfft = 21;
Fs=100;
XN = fft(x,Nfft);
Pxx = 10*log10(abs(XN.^2)/(Nfft+1));
f =Fs*(0:length(Pxx)-1)/length(Pxx);
plot(f,Pxx);
xlabel('频率'),ylabel('功率/dB'),title('功率谱');grid;
得到的功率谱图像为
由数学模型，知该序列为MA(1)序列，AR法估计该序列的功率谱程序如下：
Fs=100; N=20;
k=0:1/Fs:N/Fs;
x=[4.200 5.800 6.900 7.620 5.570 3.340 2.000 1.700 2.020 2.710 3.630 5.180 7.110 8.260 7.960 6.780 5.070 5.040 6.020 7.610 10.320];
R=zeros(1,3);
for m=1:3
RXm=0;
for k=1:N+2-m
RX=x(k+m-1)*x(k);
RXm=RXm+RX;
end
R(m)=RXm/(N+1-m);
end
a=zeros(2,2); sigma2=zeros(1,3);
sigma2(1)=R(1); a(1,1)=-R(2)/R(1);
sigma2(2)=(1-(abs(a(1,1)))^2)*sigma2(1);
a(2,2)=-(R(3)+a(1,1)*R(2))/sigma2(2);
a(2,1)=a(1,1)+a(2,2)*a(1,1);
sigma2(3)=(1-(abs(a(2,2))^2))*sigma2(2);
H=0; W=0:0.01:pi;
for k=1:2
H=H+a(2,k).*exp(-j*k*W);
end
Sx=sigma2(3)./((abs(1+H)).^2);
f=W*Fs/(2*pi);
figure(1);
plot(f,10*log10(Sx));grid
得到的功率谱图像为
与DFT法相比，图像更为平滑。

（2）DFT法估计该序列的功率谱程序如下：
x = [10.5 10.1 8.8 9.9 11.3 12.2 11.3 9.8 9.7 9.5 10 8.9 8.2 10.2 8.8 8.4 9.6 10.2 10.6 11.1 4.7];
Nfft = 21;
Fs=100;
XN = fft(x,Nfft);
Pxx = 10*log10(abs(XN.^2)/(Nfft+1));
f =Fs*(0:length(Pxx)-1)/length(Pxx);
plot(f,Pxx);
xlabel('频率'),ylabel('功率/dB'),title('功率谱');grid;
得到的功率谱图像为
由数学模型，知该序列为MA(1)序列，MA法估计该序列的功率谱程序如下：
Fs=100; N=20;
x=[10.5 10.1 8.8 9.9 11.3 12.2 11.3 9.8 9.7 9.5 10 8.9 8.2 10.2 8.8 8.4 9.6 10.2 10.6 11.1 4.7];
R=zeros(1,2);
for m=1:2
RXm=0;
for k=2:N+1
RX=x(k+1-m)*x(k);
RXm=RXm+RX;
end
R(m)=RXm/N;
end
b=-(2*sqrt((R(1))^2-(R(2))^2)-2*R(1))/(2*R(2));
sigma2=R(1)+0.5*b*R(2);
W=0:0.01:pi;
H=b.*exp(-j*W);
Sx=sigma2.*((abs(1+H)).^2);
f=W*Fs/(2*pi);
figure(1);
plot(f,10*log10(Sx));grid
得到的功率谱图像为
与DFT 法相比，图像更为平滑。

5-6 设观测方程为
6sin(0.4π)sin(0.43π),(0,1,,127)k k x k k e k =++=
其中，e k 是均值为零、方差为1的高斯白噪声。

试用最小二乘法估计随机序列x k 的AR 模型参数（模型阶次分别取4和6）；并分别用DFT 和AR 谱估计方法，估计正弦信号分量的频率及其统计结果（均值和方差）。

解：程序如下（利用burg 算法实现）
clear all
N =512;
f1 = 0.20;
f2 = 0.215;
n= 1:N;
x(n) = 6*sin(2*pi*n*f1) + sin(2*pi*n*f2)+randn(size(n));
subplot(2,1,1);
plot(n,x(n));xlabel('n');ylabel('x(n)');
title('两个正弦信号与白噪声叠加的时域波形');
p=input('Input a Number > ')
ef=zeros(p,N);
eb=zeros(p,N);
de=zeros(p,p);
ef(1,:)=x(n);
eb(1,:)=x(n);
cov(1)=x*x'/N;
k(1)=0;
for m=2:p+1
mol=0;
den=0;
for n=m:N
mol= mol+ (-2)*ef(m-1,n)*eb(m-1,n-1);
den=den+(ef(m-1,n))^2+(eb(m-1,n-1))^2;
end
k(m)=mol./den;
de(m,m)=k(m);
h(1)=cov(1)*(1-de(1,1)^2);
for n=m:N
ef(m,n)=ef(m-1,n)+k(m).*eb(m-1,n-1);
eb(m,n)=eb(m-1,n-1)+k(m).*ef(m-1,n);
end
end
k=k(1,2:p+1);
de=de(2:p+1,2:p+1);
for m=2:p
for i=1:m-1
de(m,i)=de(m-1,i)+k(m)*de(m-1,m-i);
h(m)=h(m-1)*(1-k(m)^2);
end
end
z=de(p,:);
s=[1,z];
n=(0:511)/512;
Hw=fft(s,512);
subplot(2,1,2);
plot(n,h(p)./abs(Hw).^2);
xlabel('频率（Hz）');
ylabel('10log(PSD)');
title('基于burg算法的功率谱估计');
运行结果：
P=4
P=6。