Harbin Institute of Technology
实验报告
课程名称： 随机信号分析
院 系： 电子与信息工程学院
班 级：
姓 名：
学 号：
实验时间：
哈尔滨工业大学
实验一
一.产生均匀分布的随机数，检验均值和方差。
1.程序
a=random('unif',1,5,1,2000);
b1=mean(a) %为样本均值
b2=var(a) %为样本方差
figure
subplot(211),plot(a);
subplot(212),hist(a);
2.结果
实验分析：
实验得到的平均数为：2.997；方差为：1.2940，理论的平均数应该是：3，方差为：1.33.与仿真结果相近，误差主要原因是产生的随机变量是有限个，随着随机数数量的增加，平均数和方差会更加接近理论值。
二.产生均值为0，方差为1的高斯随机数。
1.程序
b=-4:0.1:4;
a=random('normal',0,1,1,2000);
b1=mean(a)
b2=var(a)
figure
subplot(211),plot(a);
subplot(212),hist(a,b);
2.结果
实验分析：
实验得到的平均数为：0.0366；方差为：1.0013.仍与理论值有偏差，误差主要原因是产生的随机变量是有限个，随着随机数数量的增加，平均数和方差会更加接近理论值。
三．指数分布
1.程序
function [aver,sigema]=exp
N=20000;
G1=random('Normal',0,1,1,N);
G2=random('Normal',0,1,1,N);
Exp=G1.*G1+G2.*G2;
aver=mean(Exp);
sigema=var(Exp);
g=0:0.1:20;
subplot(211);plot(Exp);
subplot(212);hist(Exp,g);
end
实验结果如下：
实验分析：
MATLAB计算所得的均值和方差分别为，理论计算得均值和方差分别为，误差很小。
四．互相独立的，在（0,1）内均匀分布的随机数作和来得到高斯分布。
1.程序
X1=random('unif',0,1,1,1024);
X2=random('unif',0,1,1,1024);
X3=random('unif',0,1,1,1024);
X4=random('unif',0,1,1,1024);
X5=random('unif',0,1,1,1024);
X6=random('unif',0,1,1,1024);
Y1=X1+X2; %两组随机数作和
Y2=X1+X2+X3+X4+X5+X6; %六组随机数作和
Y3=random('normal',0,1,1,1024); %一个高斯数列与之作对比
subplot(311);hist(Y1,0:0.05:2);
subplot(312);hist(Y2,0:0.05:6);
subplot(313);hist(Y3,-3:0.05:3);
结果
实验分析：
仿真显示当两个独立的均匀分布随机变量相加时得到的分布情况接近于一个三角形，与理论计算相同。如果产生的随机变量是无穷多个即可看出两个均匀分布的随机变量相加得到的是一个标准三角形。
仿真显示当多个独立的均匀分布随机变量相加时得到的分布情况接近于高斯分布，如果均匀分布的随机数组数增加，每组的随机变量的个数增加，就会更加接近高斯分布。
五．用N（0,1）高斯分布随机数分别仿真瑞利分布随机变量和4个和7个自由度的 分布随机变量，并画出直方图。
1.程序
N=3000;
G1=random('Normal',0,1,1,N);
G2=random('Normal',0,1,1,N);
G3=random('Normal',0,1,1,N);
G4=random('Normal',0,1,1,N);
G5=random('Normal',0,1,1,N);
G6=random('Normal',0,1,1,N);
G7=random('Normal',0,1,1,N);
R=sqrt(G1.*G1+G2.*G2);
X2=G1.*G1+G2.*G2+G3.*G3+G4.*G4;
X3=G1.*G1+G2.*G2+G3.*G3+G4.*G4+G5.*G5+G6.*G6+G7.*G7;
subplot(411);hist(G1,-5:0.05:5);
subplot(412);hist(R,0:0.05:4);
subplot(413);hist(X2,0:0.1:20);
subplot(414);hist(X3,0:0.1:30);
mean1=mean(G1);
var1=var(G1);
mean2=mean(R);
var2=var(R);
mean3=mean(X2);
var3=var(X2);
mean4=mean(X3);
var4=var(X3);
2.结果
实验分析：
可以看出分布随着参加运算的高斯分布数量的增加，分布会更加趋近于高斯分布。
实验二
一．仿真一个均值为0，方差为1的高斯随机过程样本，并用直接法求自相关函数及功率谱。
自相关函数
N1=256;
xn1=random('norm',0,1,1,N1);
Rx1=xcorr(xn1,'biased');
m1=-N1+1:N1-1;
figure(1)
subplot(211),plot(m1,Rx1);
axis([-N1+1 N1-1 -0.5 1.5]);
N2=1024;
xn2=random('norm',0,1,1,N2);
Rx2=xcorr(xn2,'biased');
m2=-N2+1:N2-1;
subplot(212),plot(m2,Rx2);
axis([-N2+1 N2-1 -0.5 1.5]);
功率谱程序
N1=256;
fs=1000;
xn1=random('norm',0,1,1,N1);
Sk1=periodogram(xn1);
figure(1)
subplot(211),plot(10*log10(Sk1));
N2=1024;
xn2=random('norm',0,1,1,N2);
Sk2=periodogram(xn2);
subplot(212),plot(10*log10(Sk2));
3.修改均值为1，即仿真一个均值为1，方差为1的高斯随机过程样本
N1=256;
xn1=random('norm',1,1,1,N1);
Rx1=xcorr(xn1,'biased');
m1=-N1+1:N1-1;
figure(3)
subplot(211),plot(m1,Rx1);
axis([-N1+1 N1-1 -0.5 2.0]);
N2=1024;
xn2=random('norm',1,1,1,N2);
Rx2=xcorr(xn2,'biased');
m2=-N2+1:N2-1;
subplot(212),plot(m2,Rx2);
axis([-N2+1 N2-1 -0.5 2.0]);
功率谱程序
N1=256;
xn1=random('norm',5,1,1,N1);
Sk1=periodogram(xn1);
figure(3)
subplot(211),plot(10*log10(Sk1));
N2=1024;
xn2=random('norm',5,1,1,N2);
Sk2=periodogram(xn2);
subplot(212),plot(10*log10(Sk2));
二．条件同一，用FFT估计法求自相关函数
N1=256;
xn1=random('norm',0,1,1,N1);
Xk=fft(xn1,2*N1);
Rx=ifft((abs(Xk).^2)/N1);
figure(2)
m=-N1:N1-1;
subplot(211),plot(m,fftshift(Rx));
axis([-N1 N1-1 -0.5 1.5]);
N2=1024;
xn2=random('norm',0,1,1,N2);
Xk=fft(xn2,2*N2);
Rx=ifft((abs(Xk).^2)/N2);
m=-N2:N2-1;
subplot(212),plot(m,fftshift(Rx));
axis([-N2 N2-1 -0.5 1.5]);
三．均匀分布的自相关函数和功率谱密度
程序
N1=256;
xn1=random('unif',1,5,1,N1);
Rx1=xcorr(xn1,'biased');
m1=-N1+1:N1-1;
figure(1)
subplot(211),plot(m1,Rx1);
axis([-N1+1 N1-1 -0.5 15]);
N2=1024;
xn2=random('unif',1,5,1,N2);
Rx2=xcorr(xn2,'biased');
m2=-N2+1:N2-1;
subplot(212),plot(m2,Rx2);
axis([-N2+1 N2-1 -0.5 15]);
功率谱密度
程序
Sk1=abs(fft(Rx1));
Sk2=abs(fft(Rx2));
figure(2)
subplot(211),plot(10*log10(Sk1));
subplot(212),plot(10*log10(Sk2))； %运行时需要和上面的程序一起运行。
实验三
一．相位均匀分布的随机余弦信号和高斯白噪声之和的自相关函数估计
实验产生的随机信号是，其中为[0,]内均匀分布的随机变量，是数学期望为0、方差为1的高斯白噪声。
实验程序如下：
Function 1
N=1024; %每个样本函数的数据长度
t=0:N-1; %时间变量
m=-N:N-1; %自相关函数的时间变量
x1n=random('norm',0,1,N,8); %产生个高斯随机数，作为8个样本序列
x2n=zeros(N,8);
xn=x2n;
X1k=fft(x1n,2*N); %对x1n做2N点FFT
R1x=ifft((abs(X1k).^2)/N); %反变换，得到白噪声的自相关函数
A=random('unif',0,1,1,8)*2*pi; %产生8个[0,]内均匀随机数，作为初始相位
for k=1:8
x2n(:,k)=cos(2*pi*4*t(:)/N+A(k)); %产生余弦信号，随机相位为A(k)
xn(:,k)=x1n(:,k)+x2n(:,k); %噪声与信号之和
end
X2k=fft(x2n,2*N);
R2x=ifft((abs(X2k).^2)/N); %得到随机相位余弦信号的自相关函数
Xk=fft(xn,2*N);
Rx=ifft((abs(Xk).^2)/N); %得到噪声与信号之和的自相关函数
subplot(311);plot(m,fftshift(R1x)); %将数据位移后画图
axis([-N N-1 -0.5 1.5]); %图的坐标范围
title('高斯白噪声自相关函数');
subplot(312);plot(m,fftshift(R2x));
axis([-N N-1 -0.5 1.5]);
title('随机相位余弦自相关函数');
subplot(313);plot(m,fftshift(Rx));
axis([-N N-1 -0.5 1.5]);
title('二者相加自相关函数');
end
实验分析：
图示的自相关函数虽然是时间自相关函曲线，没有做统计平均，但是我们认为这三个过程都是各态历经过程，所以我们可以用其中的一个样本来估计整个过程的情况。可以看出白噪声是非周期平稳序列，随机相位余弦信号则是周期平稳序列。
二．相位均匀分布的两个随机余弦信号与高斯白噪声之和的功率谱密度估计
Function 2
N=1024; %序列长度
fs=2000; %采样频率
t=(0:N-1)/fs; %时间序列
fai=random('unif',0,1,1,2)*2*pi; %产生2个[0,]内均匀随机数
xn=cos(2*pi*40*t+fai(1))+3*cos(2*pi*120*t+fai(2))+rand(1,N);
%产生含噪声的随机变量
Sx=abs(fft(xn).^2)/N; %估计功率谱
f=(0:N-1)*fs/N; %频率轴坐标
plot(f,10*log10(Sx(1:N))); %用dB/Hz做功率谱单位，画图
axis([0 2000 -40 40]); %图的坐标范围
xlabel('f/Hz');
ylabel('dB/Hz');
End
结果
实验分析
从图中可以看出，功率谱大致在40Hz和120Hz处有两个谱峰，与两个余弦的频率相同，符合理论值。
三．信号和噪声同时作用于带通滤波器
1. 信号在带通滤波器通带内
实验产生的随机信号是，其中，和为分布在[0,]内均匀分布的随机变量，是数学期望为0、方差为1的高斯白噪声。所用的带通滤波器是100阶、截止频率为0.05（50Hz）和0.5（500Hz）的带通滤波器。
实验程序如下：
function 3
N=1024; %样本长度
fs=2000; %抽样频率
t=(0:N-1)/fs; %时间序列
fai=random('unif',0,1,1,2)*2*pi; %产生2个[0,]内均匀随机数
x=8*cos(2*pi*200*t+fai(1))+3*cos(2*pi*300*t+fai(2))+rand(1,N);
%产生含噪声的随机序列
h=fir1(100,[0.05,0.5]); %产生100阶，截止频率为0.1和0.4的带通滤波器
H=fft(h,2*N); %2N点滤波器频率响应（系统传输函数）
Rx=xcorr(x,'biased'); %直接估计含噪声的随机序列的自相关函数
Sx=abs(fft(x,2*N).^2)/(2*N); %估计含噪声的随机序列的功率谱密度
HW=abs(H).^2; %系统的功率传输函数
Sy=Sx.*HW; %输出信号的功率谱
Ry=fftshift(ifft(Sy)); %用IFFT求输出信号的自相关函数
f=(0:N-1)*fs/N;
w=(1:N)/N;
t1=(-N:N-1)/N*(N/20000);
subplot(411);plot(f,10*log10(Sx(1:N))); %画出输入信号功率谱密度
title('输入信号功率谱密度');
subplot(412);plot(w,abs(HW(1:N))); %画出系统功率传输函数
title('系统功率传输函数');
subplot(413);plot(f,10*log10(Sy(1:N))); %画出输出信号功率谱密度
title('输出信号功率谱密度');
subplot(414);plot(t1,Ry); %画出输出信号自相关函数
title('输出信号自相关函数');
实验结果如下：
其他条件不变只是将滤波器的阶数变为1000阶
实验分析：
从图中可以看到，当信号在带通滤波器通带内的时候，输出信号功率谱密度在200Hz和300Hz处有谱峰，基本符合理论计算。同时在滤波器的阶数增加时可以看到在带外的功率谱密度下降的更多，也就是说阶数的增加会增强滤波器的性能。
2.信号在带通滤波器通带外
实验产生的随机信号是，其中，和为分布在[0,]内均匀分布的随机变量，是数学期望为0、方差为1的高斯白噪声。所用的带通滤波器是1000阶、截止频率为0.2和0.4的带通滤波器。
实验程序如下：
function 4
N=1024; %样本长度
fs=2000; %抽样频率
t=(0:N-1)/fs; %时间序列
fai=random('unif',0,1,1,2)*2*pi; %产生2个[0,]内均匀随机数
x=8*cos(2*pi*500*t+fai(1))+3*cos(2*pi*700*t+fai(2))+rand(1,N);
%产生含噪声的随机序列
h=fir1(1000,[0.2,0.4]); %产生100阶，截止频率为0.1和0.4的带通滤波器
H=fft(h,2*N); %2N点滤波器频率响应（系统传输函数）
Rx=xcorr(x,'biased'); %直接发估计含噪声的随机序列的自相关函数
Sx=abs(fft(x,2*N).^2)/(2*N); %周期法估计含噪声的随机序列的功率谱密度
HW=abs(H).^2; %系统的功率传输函数
Sy=Sx.*HW; %输出信号的功率谱
Ry=fftshift(ifft(Sy)); %用IFFT求输出信号的自相关函数
f=(0:N-1)*fs/N;
w=(1:N)/N;
t1=(-N:N-1)/N*(N/20000);
subplot(411);plot(f,10*log10(Sx(1:N))); %画出输入信号功率谱密度
title('输入信号功率谱密度');
subplot(412);plot(w,abs(HW(1:N))); %画出系统功率传输函数
title('系统功率传输函数');
subplot(413);plot(f,10*log10(Sy(1:N))); %画出输出信号功率谱密度
title('输出信号功率谱密度');
subplot(414);plot(t1,Ry); %画出输出信号自相关函数
title('输出信号自相关函数');
实验结果如下：
实验误差分析：
由第一个实验可得，用1000阶滤波器可以有很好的效果。从图中可以看出，当信号在带通滤波器通带外的时候，输出信号功率谱密度在500Hz和700Hz处谱峰有很大的衰减。