随机信号分析编程作业

随机信号分析
编程作业
姓名：
学号：
学院：计算机与信息学院班级：通信工程14-1班
1.23编写一个产生均值为1、方差为4的高斯分布随机数程序，求其最大值、最小值、均值和方差，并与理论值比较。

解：分析:本题可用累加近似法产生标准正太分布随机数，首先产生12个相互独立的均匀分布随机数，计算这十二个数的和后减去六后就可以得到N(0，1)分布的随机数。

以下是代码及结果：
从运行结果可以看出，
产生的1024个随机
数，其均值为
0.9557，方差为
3.8582与理论均值1
方差值4较为接近，故
此程序比较理想的产生
了均值为1方差值为4
的高斯随机数。

2.26编写一个产生协方差函数为C(τ)=4e−2|τ|的平稳高斯过程的程序，产生若
干样本函数，估计所产生的时间自相关函数和功率谱密度，并统计自相关函数和功率谱密度，最后将结果和理论值比较。

解：本题可根据教材例题5.13差分方程得到自相关函数为b2
1−a2
a|m|的随机序列，代码如下：
N=10000;
Ts=0.001;
sigma=2;
beta=2;
a=exp(-beta*Ts);
b=sigma*sqrt(1-a*a);
w=normrnd(0,1,[1,N]);
x=zeros(1,N);
x(1)=sigma*w(1);
for i=2:N
x(i)=a*x(i-1)+b*w(i);
end;
Rxx=xcorr(x)/N;
m=[-N+1:N-1];
Rxx0=(sigma^2)*exp(-beta*abs(m*Ts));
plot(m*Ts,Rxx0,'b.',m*Ts,Rxx,'r');title(‘理论与实测自相关函数’);
用matlab运行后绘制的图如下：
例题3.6.4仿真一个平均功率为1的白噪声通过带通系统，白噪声为高斯分布，带通系统的俩个截至频率分别为3kHz 和4kHz，求输出的自相关函数和功率谱密度。

%准备工作N = 500;
xt = random('norm',0,1,1,N);
ht = fir1(101,[0.3 0.4]);
HW = fft(ht,2*N);
%仿真
Rxx = xcorr(xt,'biased');
Sxx = abs(fft(xt,2*N).^2)/(2*N);
HW2 = abs (HW).^2;
Syy = Sxx.*HW2;
Ryy = fftshift(ifft(Syy));
%画曲线
w = (1:N)/N;
t = (-N : N-1)/N * (N/20000);
subplot(4,1,1); plot(w,abs(Sxx(1:N)));
subplot(4,1,2); plot(w,abs(HW2(1:N)));
subplot(4,1,3); plot(w,abs(Syy(1:N)));
subplot(4,1,4); plot(Ryy);
M = 100;
N = 500;
xt= random('norm',0,1,M,N);
ht = fir1(101,[0.3 0.4]);
HW =fft (ht,2*N);
Sxx = abs(fft(xt,2*N,2).^2)/(2*N);
Sxxav = mean(Sxx);
HW2 =abs(HW).^2;
Syy = Sxxav .* HW2;
Ryy = fftshift(ifft(Syy));
运行结果见下图：
教材习题6.16编写MATLAB程序，模拟产生功率谱为
S(ω)=16
（ω+ω0）2
+64
+16
（ω−ω0）
2
+64
的高斯带通随机信号，其中ω0=400π，绘制带通信号相关函数与功率谱。

解：由傅里叶变换性质知
σ2e−β|τ|cosω0↔σ2β
（ω+ω0）2
+β2
+σ2β
（ω−ω0）
2
+β2
因此先产生俩个自相关函数同为R(τ)=σ2e−β|τ|的独立平稳过程i(t)和q(t),然后由x(t)=i(t)cos(w0t)-q(t)sinw0t就可得出要求的带通信号。

代码如下：
NFFT = 1024;
fs = 1000;
Ts = 1/fs;
B= 0.5 * fs;
df = fs/NFFT;
f = -B:df:B-df;
sigma = sqrt(2);
beta = 8;
a = exp(-beta*Ts);
b = sigma * sqrt(1 - a*a);
f0 = 200;
N = 10000;
wi = normrnd (0,1,[1,N]);
wq = normrnd (0,1,[1,N]);
xi = zeros(1,N);
xq = zeros(1,N);
xi(1) = sigma*wi(1);
xq(1) = sigma * wq(1);
for i = 2:N
xi(i)=a*xi(i-1) + b*wi(i); xq(i)=a*xq(i-1) + b*wq(i); end
t = 0:Ts:(N-1)*Ts;
x0=xi.*cos(2*pi*f0*t) -
xq.*sin(2*pi*f0*t);
x = x0(N-1000+1:N);
t = t(N-1000+1:N);
Rxx = xcorr(x)/1000;
t0=[-1000+1:1000-1];
Rxx0=(sigma^2)*exp(-
beta*abs(t0*Ts));
subplot(141);plot(t,x);title('样本x');
subplot(142);plot(t0,Rxx);title( '实测自相关函数');
subplot(143);periodogram(x,[],10 00,fs);
subplot(144);plot(t0,Rxx0,'r');t itle('理论自相关函数');
以下是运行结果：。