哈尔滨工业大学威海随机信号分析实验一报告.docx
《哈尔滨工业大学威海随机信号分析实验一报告.docx》由会员分享,可在线阅读,更多相关《哈尔滨工业大学威海随机信号分析实验一报告.docx(30页珍藏版)》请在冰豆网上搜索。
哈尔滨工业大学威海随机信号分析实验一报告
《随机信号分析》实验报告
班级:
学号:
姓名:
《随机信号分析》实验报告
实验一
1、实验目的:
熟悉并练习使用随机信号Matlab的函数
二、实验内容:
1、熟悉并练习使用下列Matlab的函数,给出各个函数的功能说明和内部参数的意义,并给出至少一个使用例子和运行结果:
1)rand()11)unifpdf()
2)randn()12)unifcdf()
3)normrnd()13)raylpdf()
4)mean()14)raylcdf()
5)var()15)exppdf()
6)xcorr()16)expcdf()
7)periodogram()17)chol()
8)fft()18)ksdensity()
9)normpdf()19)hist()
10)normcdf()20)int()
用法、功能、程序如下:
1)randn(m,n)
功能:
返回一个从标准正态分布中得到的伪随机标量。
>>r=randn(5)%由标准正态分布随机数组成的5×5矩阵。
r=
-1.0689-0.75490.31920.6277-1.2141
-0.80951.37030.31291.0933-1.1135
-2.9443-1.7115-0.86491.1093-0.0068
1.4384-0.1022-0.0301-0.86371.5326
0.3252-0.2414-0.16490.0774-0.7697
2)rand(m,n)
功能:
返回一个从开区间(0,1)上的标准均匀分布得到的伪随机标量。
r=rand(5)%生成一个由介于0和1之间的均匀分布的随机数组成的5×5矩阵
>>r=
0.54690.95720.91570.84910.3922
0.95750.48540.79220.93400.6555
0.96490.80030.95950.67870.1712
0.15760.14190.65570.75770.7060
0.97060.42180.03570.74310.0318
3)normrnd(mu,sigma,m,n)
功能:
以均值μ和标准差σ为参数的正态分布随机数mxn
>>normrnd(0,1,3,4)%生成均值μ=0,σ=1的3x4正态分布随机数
ans=
0.27610.3919-0.74110.0125
-0.2612-1.2507-0.5078-3.0292
0.4434-0.9480-0.3206-0.4570
4)mean(A,dim)
功能:
数组的平均值
mean(A,dim)dim=1,返回列平均数,默认为1
dim=2,返回列平均数
dim>2,返回A
A=[011;232;132;422]%M=mean(A)沿A的大小不等于1的第一个数组维度返回均值。
M=mean(A)
>>A=
011
232
132
422
M=
1.75002.25001.7500
mean(A,2)
ans=
0.6667
2.3333
2.0000
2.6667
5)var(A,dim)
功能:
方差
var(A,dim)dim=0,按N-1算,默认是0
dim=1,按N算
a=
123
>>var(a,1)%按N=3算
ans=
0.6667
>>var(a,0)%按N=2算
ans=
1
6)xcorr(x,y)
功能:
互相关
%返回两个离散时间序列的互相关函数,返回长度为2*N-1互相关序列,其中x,y均为N
%------------------------------------------------------------------
x=[1,2i,3];y=[4,5,6];%离散序列的互相关函数
[c1,lags]=xcorr(x,y);
c1=mat2str(c1,2),lags
>>c1=
[6-8.9e-16i5+12i22+10i15+8i12+8.9e-16i]
lags=
-2-1012
%------------------------------------------
t=[0:
1:
100];%余弦序列自相关
x=cos(t);
[a,b]=xcorr(x);
stem(b,a)
图像如下:
7)[PXX,W]=periodogram(x)
功能:
x是输入。
周期图法(矩形窗截取)估计的功率谱密度pxx,如果输入x是实数值,返回的是单边PSD,w=(0,pi);如果输入x是复数值,w=(0,2pi),返回的是双边PSD。
%获取一个输入信号由一个(pi/4)/弧度/样本离散正弦周期图(0,1)白噪声的功率谱密度。
n=0:
319;
x=cos(pi/4*n)+randn(size(n));
[pxx,w]=periodogram(x);
plot(w,10*log10(pxx));
periodogram(x);
图像如下:
8)fft(x,n)
功能:
快速傅里叶变换,x是输入,返回N点DFT
(a)
x1=[13536839];
fft_x1=fft(x1,8);
x_axis=[0:
1:
7];
figure
(1);
subplot(3,1,1);
stem(x_axis,x1,'.');
xlabel('n');title('时间序列');
subplot(3,1,2);
stem(x_axis,abs(fft_x1),'x');gridon;
xlabel('频率k');ylabel('幅度');
title('幅度谱');
subplot(3,1,3);
stem(x_axis,angle(fft_x1));gridon;
xlabel('频率k');ylabel('相位');
title('相位谱');
图像如下:
(b)
clear;
Fs=1000;%采样频率
T=1/Fs;%采样时间
L=1000;%信号长度
t=(0:
L-1)*T;%时间向量
%50Hz正弦信号与120Hz正弦信号之和
x=0.7*sin(2*pi*50*t)+sin(2*pi*120*t);
y=x+2*randn(size(t));%有用信号加噪声
figure
(1)
plot(Fs*t(1:
50),y(1:
50))
title('零均值随机噪声信号')
xlabel('时间(毫秒)')
NFFT=2^nextpow2(L);%比L更大的2的指数的幂次
Y=fft(y,NFFT)/L;%对y做1024点的DFT,除以L对幅值作限制(因为做NFFT,幅值有叠加)
f=(Fs/2)*linspace(0,1,NFFT/2+1);%'linspace(0,1,NFFT/2+1)'在(0,1)内返回(NFFT/2+1)个线性间隔点
%绘制单边振幅谱
figure
(2)
plot(f,2*abs(Y(1:
NFFT/2+1)))%将负频映射到正频,绘制单边谱
title('单边幅度谱图')
xlabel('频率(Hz)')
ylabel('|Y(f)|')
%------------------------------------------
nextpow2(2的更高次幂的指数)
P=nextpow2(A)返回对A中每个元素满足:
(2^p)>=|A|
的最小的2的幂的指数。
您可以使用nextpow2填充传递到fft的信号。
当信号长度并非2次幂时,这样做可以加快FFT的运算速度。
%------------------------------------------
图像如下:
9)normpdf(x,mu,sigma)默认mu=0,sigma=1
功能:
返回计算正态分布概率密度函数在x点的值,x可以是矢量
normpdf(2,0,1)%计算标准正态分布的pdf在x=2时的值
ans=
0.0540
10)normcdf(x,mu,sigma)默认mu=0,sigma=1
功能:
返回正态分布累积分布函数在x点的值,x可以是矢量
%P=normcdf(x)返回值x的标准正态分布的参数μ=0,σ=1标准正态累积分布函数
p=normcdf([-11]);%计算x=[-11]两点的cdf的值
p
(2)-p
(1)
ans=
0.6827
11)unifpdf(x,a,b)
功能:
返回计算[a,b]区间连续均匀分布概率密度函数在x处的值
x=0.1:
0.1:
0.6;%求对应的概率密度
y=unifpdf(x)
y=
111111
12)unifcdf(x,a,b)
功能:
计算[a,b]上均匀分布累积分布函数F(x),默认[0,1]。
>>probability=unifcdf(0.75)
probability=
0.7500
13)raylpdf(x,b)
功能:
返回计算瑞利概率密度函数在参数b下的在x处的值
x=[0:
0.01:
2];%sigma=0.5的瑞利分布
p=raylpdf(x,0.5);
figure;
plot(x,p);
title('瑞利分布(b=0.5)')
图像如下:
14)raylcdf(x,b)
功能:
瑞利分布参数为b时累积分布函数
x=0:
0.1:
3;%sigma=1的瑞利累积分布
p=raylcdf(x,1);
plot(x,p)
title('瑞利分布cdf(b=1)')
图像如下:
15)exppdf(x,mu)
功能:
返回计算指数概率密度函数在参数为mu时在x处的值
y=exppdf(5,5)%mu=5对应的指数概率密度值
y=
0.0736
16)expcdf(x,mu)
功能:
返回计算指数累积分布函数在参数为mu时在x处的值
x=0:
0.1:
30;
p=expcdf(x,2);
plot(x,p);
title('指数分布cdf(\mu=2)')
图像如下:
17)chol()
功能:
Cholesky矩阵分解
%-----------------------------------------------------------------------------------------
R=chol(A)
从矩阵A的对角线和上三角生成一个上三角矩阵R,满足R'*R=A。
下三角被认为是上三角的(复共轭)转置,矩阵A必须是正定的,否则,MATLAB显示错误信息。
L=chol(A,'lower')
从矩阵A的对角线和下三角生成一个下三角矩阵L,满足L*L'=A。
当A为稀疏矩阵,chol是非常快的。
矩阵A必须是正定的,否则,MATLAB显示错误信息。
[R,p]=chol(A)
对于正定矩阵A,从矩阵A的对角线和上三角生成一个上三角矩阵R,满足R'*R=A并且p=0。
对于A不是正定的,则p是正整数,MATLAB不产生任何信息。
当A是全矩阵时,R是阶数为q=p-1的上三角矩阵并有R'*R=A(1:
q,1:
q)。
book.iLoveM
当A是稀疏矩阵时,R是大小为q*n的上三角矩阵。
[L,p]=chol(A,'lower')
对于正定矩阵A,从矩阵A的对角线和下三角生成一个下三角矩阵L,满足L*L'=A并且p=0。
对于A不是正定的,则p是正整数,MATLAB不产生任何信息。
当A是全矩阵时,R是阶数为q=p-1的下三角矩阵并有L*L'=A(1:
q,1:
q)。
当A是稀疏矩阵时,R是大小为q*n的下三角矩阵。
[R,p,S]=chol(A)
当A是稀疏矩阵时,返回一个置换矩阵S。
当p=0时,R是一个上三角矩阵并有R'*R=S'*A*S。
当p不为0时,R大小为q*n的上三角矩阵。
[R,p,s]=chol(A,'vector')
当p=0时,返回置换信息作为矩阵s并有A(s,s)=R'*R。
可以用'matrix'代替'vector'来包含默认情况。
[L,p,s]=chol(A,'lower','vector')
当p=0时,仅仅用矩阵A的对角线和下三角并返回一个下三角矩阵L和置换矩阵S有A(s,s)=L*L'。
可以用'matrix'代替'vector'来包含置换矩阵。
%---------------------------------------------------------------------------------------------
A=gallery('moler',5)
A=
1-1-1-1-1
-12000
-10311
-10142
-10125
C=chol(A)
ans=
1-1-1-1-1
01-1-1-1
001-1-1
0001-1
00001
18)[f,xi]=ksdensity(x)
功能:
核心平滑密度估计,x是输入。
计算样本向量x的概率密度估计,返回在xi点的概率密度f,此时我们使plot(xi,f)就可以绘制出概率密度曲线。
该函数,首先统计样本x在各个区间的概率,再自动选择xi,计算对应的xi点的概率密度f=ksdensity(x,xi)。
rngdefault%forreproducibility
x=[randn(30,1);5+randn(30,1)];%由正态随机数估计概率密度函数曲线
[f,xi]=ksdensity(x);
figure
plot(xi,f);
title('x的pdf估计')
图像如下:
19)hist()
功能:
hist是histogram的缩写,直方图。
1、N=hist(Y)
功能:
将向量Y的元素平均到十个等间隔的容器中,并且返回每个容器的元素
个数;如果Y是一个矩阵,hist指令逐列元素操作
2、N=hist(Y,M)
功能:
使用M个箱子
3、N=hist(Y,X)
功能:
X是向量,以X中的元素为区间中心可获得一系列的区间,执行命令可获得Y在这些区间中的分布情况。
>>x=[02925873194358100129510];
hist(x)
%直方图的横轴在x最小和最大值之间分成10等分,每个矩形的高度表示在本单元数
%第一个矩形的高度为0<=x<=1的值的个数,第二个矩形高度为1图像如下:
20)int()
功能:
求定积分和不定积分(换句话说,可用来求cdf、pdf表达式)
(a)
>>symsx%求表达式的积分
int(-3*x^2+8*x)
ans=
-x^2*(x-4)
(b)
>>symsx%把这个表达式从0到1求定积分
int(x*log(1+x),0,1)
ans=
1/4
2、产生高斯随机变量
(1)产生数学期望为0,方差为1的高斯随机变量;
(2)产生数学期望为5,方差为10的高斯随机变量;
(3)利用计算机求上述随机变量的100个样本的数学期望和方差,并与理论值比较;
程序如下:
(1)X=randn
X=
0.0326
或
X=normrnd(0,1)
X=
0.5525
(2)Y=normrnd(5,sqrt(10))
Y=
8.4804
或
Y=5+sqrt(10)*randn
Y=
9.8832
(3)
x=randn(1,100);
y=normrnd(5,sqrt(10),[1100]);
mx1=mean(x)%计算均值
mx2=mean(y)
Dx1=var(x)%计算方差
Dx2=var(y)
xmx=abs(mx1-0)%与理论值进行比较,xmi(i=x,y),xDi(i=x,y)为计算值与理论值的差值的绝对值
xmy=abs(mx2-5)
xDx1=abs(Dx1-1)
xDx2=abs(Dx2-10)
mx1=
-0.0690
mx2=
4.7735
Dx1=
0.9781
Dx2=
8.9237
xmx=
0.0690
xmy=
0.2265
xDx1=
0.0219
xDx2=
1.0763
%-------------------------------------
%用symsum或int计算均值的方法
RPX的分布列为:
P{X=k}=(1/2)^k%k=1,2,3,...,n,...
求EX
%-------------------------------------
%由均值的计算公式,程序如下:
symsk;
symsum(k*(1/2)^k,k,1,inf);%1,inf为求和上下限(inf为无穷大)
3、产生
分布的随机变量
(1)产生自由度为2,数学期望为2,方差为4的具有中心
分布的随机变量;
(2)产生自由度为2,数学期望为4,方差为12的具有非中心
分布的随机变量;
(3)利用计算机求上述随机变量的100个样本的数学期望和方差,并与理论值比较;
程序如下:
(1)
x1=normrnd(2,2)
x2=normrnd(2,2)
y=(x1).^2+(x2).^2
x1=
2.6010
x2=
1.2539
y=
8.3372
(2)
x1=normrnd(2,sqrt(12))%由非中心卡方分布的均值与方差公式
x2=normrnd(2,sqrt(12))
y=(x1).^2+(x2).^2
x1=
4.8249
x2=
4.7674
y=
46.0083
(3)
(a)
x1=normrnd(2,2,1,100);
x2=normrnd(2,2,1,100);
y=(x1).^2+(x2).^2;
y_m=mean(y)
y_sigma=var(y)
y_m=
17.9989
y_sigma=
262.3949
(b)
x1=normrnd(2,sqrt(12));%由非中心卡方分布的均值与方差公式
x2=normrnd(2,sqrt(12));
y=(x1).^2+(x2).^2;
y_m=mean(y)
y_sigma=var(y)
y_m=
99.9632
y_sigma=
0
4、利用Matlab现有pdf和cdf函数,画出均值为零、方差为4的高斯随机变量的概率密度
曲线和概率分布曲线。
程序如下:
x=-10:
0.1:
10;
y_pdf=normpdf(x,0,2);
y_cdf=normcdf(x,0,2);
subplot(2,1,1)
plot(x,y_pdf,'b');
title('高斯随机变量的概率密度曲线');
subplot(2,1,2)
plot(x,y_cdf,'r');
title('高斯随机变量的累积分布曲线');
图像如下:
5、产生长度为1000数学期望为5,方差为10的高斯随机序列,并根据该序列值画出其概率密度曲线。
(不使用pdf函数)
程序如下:
x=normrnd(5,sqrt(10),[11000]);
[f,xi]=ksdensity(x);
plot(xi,f);
title('\mu=5,\sigma=sqrt(10)的高斯随机分布的pdf')
图像如下:
6、利用Matlab求随机变量的统计特性
MATLAB程序如下:
symsxyA;
f=A*exp(-x-y);
C=int(int(f,x,0,inf),y,0,inf);%计算分布函数的最大值,最大值C=1,由此可C=A=1
A=solve(C-1,A)
P1=int(int(f,y,0,1-x),x,0,1);%计算P{0fx=int(f,y,0,inf)%计算f(x)
fy=int(f,x,0,inf)%计算f(y)
symsuv;
f=C*exp(-u-v);
F=int(int(f,u,0,y),v,0,x);%求分布函数表达式
A=
1
fx=
A*exp(-x)
fy=
A*exp(-y)
6、1)
MATLAB程序如下:
symsxyA;
f=A*exp(-2*x-y);
C=int(int(f,x,0,inf),y,0,inf);%计算分布函数的最大值,C=1
A=solve(C-1,A)
P1=int(int(f,y,1,inf),x,2,inf);%计算P{X>2,Y>1}
fx=int(f,y,0,inf)%计算fx
fy=int(f,x,0,inf)%计算fy
symsuv;
f=C*exp(-2*u-v);
F=int(int(f,u,0,y),v,0,x);%求分布函数表达式
A=
2
fx=
A*exp(-2*x)
fy=
(A*exp(-y))/2
3、实验结论
通过实验了解了与随机信号分析有关的Matlab的函数,明确了各函数的功能和内部参数的意义,学会了使用Matlab产生高斯随机变量、产生
分布随机变量,利用Matlab现有pdf和cdf函数画出概率密度曲线和概率分布曲线。