MATLAB 随机过程基本理论

一、平稳随机过程

1、严平稳随机过程

clc
clearn=0:1000;
x=randn(1,1001);
subplot(211),plot(n,x);
xlabel('n');ylabel('x(n)');
grid on
subplot(212),hist(x,50)
grid on

2、

3、随机相位正弦信号

3.1

%随机相位正弦信号clc
clear
t=0:0.01:30;
w=pi/2;
x1=cos(w*t+2*pi*rand*ones(1,3001));
x2=cos(w*t+2*pi*rand*ones(1,3001));
x3=cos(w*t+2*pi*ones(1,3001));
plot(t,x1,'-.',t,x2,'--',t,x3,'-');
xlabel('n');ylabel('x(n)');
grid on
legend('样本1','样本2','样本3')

样本1和样本2的初始相位是随机生成的，样本3的初始相位为定值：2*pi。

注意一个信号只有一个初始相位，而不是说一个点一个初始相位

3.2 随机相位正弦信号均值

clc
clear
t=0:0.01:30;
N=10000;
w=2*pi;
y=zeros(1,3001);
for ii=1:Nx=cos(w*t+2*pi*rand*ones(1,3001));y=y+x;
end
y_mean=y/N;%集平均plot(t,y_mean,'-');
xlabel('n');ylabel('E[X(t)]');
grid on

可以看出1、振幅和角频率相同而初始相位不相同的正弦信号叠加还是同频率的正弦信号，只是初始相位和幅值发生了变化2、随机相位正弦信号的均值为0。

3.3 随机相位正弦信号自相关函数

clc
clear
tao=-5:0.01:5;%时间差
t=1; %选取一个时间点
N=10000;%仿真次数
w=2*pi;%角频率
y=0;
for ii=1:Nfai=rand;%随机相位x=cos(w*t+w*fai).*cos(w*t+w*tao+w*fai);%自相关函数，改变时间差y=y+x;
end
y_mean=y/N;%求平均值plot(tao,y_mean,'-');
xlabel('tao');ylabel('Rx(tao)');
grid on

更改时间点

clc
clear
tao=-5:0.01:5;%时间差
t=1.5; %选取一个时间点
N=10000;%仿真次数
w=2*pi;%角频率
y=0;
for ii=1:Nfai=rand;%随机相位x=cos(w*t+w*fai).*cos(w*t+w*tao+w*fai);%自相关函数，改变时间差y=y+x;
end
y_mean=y/N;%求平均值，循环10000次，求出一个点的一个值plot(tao,y_mean,'-');
xlabel('tao');ylabel('Rx(tao)');
grid on

时间点更改为1.5后，自相关函数没有改变，也就证明了E[cos(wt+fai)*cos(w(t+tao)+fai)]的图形与t无关，仅仅与tao有关。

4、宽平稳过程的均方遍历性

从这可以看出时间均值和集均值的差异。时间均值是一个随机函数各个值平均，类似于一个人五门课的平均分。而集平均，是指多个随机信号的平均值，类似于全班同学的平均分。

4.1 时间平均

clear all
clc
t=0:0.01:1000; %选取一个时间点
w=2*pi;%角频率
x=cos(w*t+2*pi*rand*ones(1,100001));%初始相位一直是定值
plot(t(1:500),x(1:500))
grid
Ex=mean(x') %时间平均值

可以看出时间均值和集均值都是约等于0。

4.2 自相关函数的时间平均和集平均

clc
clear
corr=[];
w=2*pi;%角频率
a=rand;
for tao=-5:0.01:5t=0:0.01:1000;x=cos(w*t+w*a*ones(1,100001)).*cos(w*t+w*tao+w*a);%自相关函数，改变时间差corr_mean=mean(x');corr=cat(1,corr,corr_mean);%将值连接起来到一个矩阵里，求出对应tao的值
end
tao=-5:0.01:5
plot(tao,corr,'-');
xlabel('tao');ylabel('corr');
grid on

自相关的时间平均和集平均相等。

所以随机相位正弦信号具有各态历经性。

二、随机信号谱分析

通过傅里叶分析来研究频域特性就是功率谱分析

1、连续周期傅里叶级数

程序实现

clc
clear
t=-2:0.0001:2;
omega=2*pi;
y=square(2*pi*t,50);%周期为1的防波信号
n_max=[1,5,20,100];
for k=1:4fk=zeros(1,length(t));for n=1:2:n_max(k)bn=4/(pi*n);%系数fk=fk+bn*sin(n*omega*t);%公式endsubplot(2,2,k)plot(t,y,t,fk,'linewidth',1.5);xlabel('t(sec)');ylabel('部分和波形');String=['最大谐波数=',num2str(n_max(k))];axis([-2 2 -3 3]);grid on;title(String);disp([String,'时，在信号跳变点附近过冲幅度（ %）']);%注意连接情况f_max=(max(fk)-max(y))*100%谐波最大值减去方波最大值
end

当谐波数分别为1、5、20、100时，在信号跳变点附近的过冲幅度分别为27.320%、18.8357%、17.9814%、17.9013%，这就是著名的Gibbs现象。该现象在一般情形下都会发生，不仅仅限于周期性方波。

将具有不连续点的周期函数（如矩形脉冲)进行傅立叶级数展开后，选取有限项进行合成。当选取的项数越多，在所合成的波形中出现的峰起越靠近原信号的不连续点。当选取的项数很大时，该峰起值趋于一个常数，大约等于总跳变值的9%。这种现象称为吉布斯效应。

2、离散周期，连续非周期，离散非周期

3、功率谱密度

随机过程的功率谱密度与确定信号频谱密度对应，不包含任何相位信息，这是由随机信号的随机性决定的。

自相关函数从时域描述随机信号，而功率谱密度从频域描述随机信号

三、随机信号自相关函数估计

3.1直接估计法

%64个标准正态分布的自相关函数有偏和无偏估计
clc
clear
N=64;
x=randn(1,N);
Rx1=xcorr(x,'unbiased');
Rx2=xcorr(x,'biased');
m=(-N+1):(N-1);
plot(m,Rx1,'r--');
hold on
plot(m,Rx2);
xlabel('m');ylabel('Rx');
axis([-N+1 N-1 -1 1.5]);
grid on;
legend('无偏估计','有偏估计');

可以看出当m比较小时，两者几乎一致；m变大时，差异明显；关于m=0对称。

2、快速傅里叶变换

clear all
clear
x=1:5;
y=reshape(1:25,5,5)
plot(x,y)

plot画图方法，按列开始画的。矩阵画法。

（1）正态分布的自相关函数

%4096点个标准正态分布的自相关函数xcorr和FFT方法比较
clc
clear
N=4096;
x=randn(1,N);tic
Rx1=xcorr(x,'biased');
time1=toc
m=(-N+1):(N-1);%个苏
subplot(121)
plot(m,Rx1);
xlabel('m');ylabel('Rx1');
axis([-N+1 N-1 -1 1.5]);
grid on;
title('xcorr命令法');tic
Xk=fft(x,2*N);
Rx2=ifft(abs(Xk).^2/N);
time2=toc
m=(-N):(N-1);%注意变化后的个数
subplot(122)
plot(m,fftshift(Rx2));%注意 fftshift的使用
xlabel('m');ylabel('Rx2');

时间FFT方法要花费少。

FFT方法比xcorr命令方法多一个点。

（2）余弦的自相关函数

clc
cleart=1:0.001:20;
N=length(t);
x=cos(2*pi*t);tic
Rx1=xcorr(x,'biased');
time1=toc
m=(-N+1):(N-1);
subplot(121)
plot(m,Rx1);
xlabel('m');ylabel('Rx1');
%axis([-N+1 N-1 -1 1.5]);
grid on;
title('xcorr命令法');tic
Xk=fft(x,2*N);
Rx2=ifft(abs(Xk).^2/N);
time2=toc
m=(-N):(N-1);
subplot(122)
plot(m,fftshift(Rx2));
xlabel('m');ylabel('Rx2');
%axis([-N+1 N-1 -1 1.5]);
grid on;
title('FFT法');

（3）随机信号与高斯白噪声信号叠加

%随机相位正弦信号与高斯白噪声信号的叠加
clc
clearN=1024;%点数，一个信号采样点
n=32;%信号数
t=0:N-1;
m=-N:N-1;
x1n=random('norm',0,1,N,n);%N行，n列的数，高斯白噪声
% x1k=fft(x1n,2*N);%2n行，32列
% R1x=ifft((abs(x1k).^2)/N);%2n行，32列
A=random('unif',0,1,1,32)*2*pi;%32个[0,2*pi]随机相位
for k=1:nx2n(:,k)=cos(2*pi*4*t(:)/N+A(k));%n个相位，随机相位正弦信号，除以N是为了使周期个数变小，只有四个周期。xn(:,k)=x1n(:,k)+x2n(:,k);%叠加信号
endsubplot(321)
plot(t,x1n(:,1));
xlabel('m');ylabel('Rn(m)');
axis([0 N -5 5]);
title('一个高斯白噪声')
grid on;subplot(322)
plot(t,x1n);
xlabel('m');ylabel('Rn(m)');
axis([0 N -5 5]);
title([num2str(n),'个高斯白噪声'])
grid on;subplot(323)
plot(t,x2n(:,1));
xlabel('m');ylabel('Rn(m)');
axis([0 N -2 2]);
title('一个随机相位正弦信号')
grid on;subplot(324)
plot(t,x2n);
xlabel('m');ylabel('Rn(m)');
axis([0 N -2 2]);
title([num2str(n),'个随机相位正弦信号'])
grid on;subplot(325)
plot(t,xn(:,1));
xlabel('m');ylabel('Rn(m)');
axis([0 N -5 5]);
title('一个叠加信号')
grid on;subplot(326)
plot(t,xn);
xlabel('m');ylabel('Rn(m)');
axis([0 N -5 5]);
title([num2str(n),'个叠加信号'])
grid on;

可以看出来，时域展示叠加信号时，混乱，找不出规律

%随机相位正弦信号与高斯白噪声信号的叠加
clc
clearN=1024;%点数，一个信号采样点
n=32;%信号数
t=0:N-1;
m=-N:N-1;
x1n=random('norm',0,1,N,n);%N行，n列的数，高斯白噪声
x1k=fft(x1n,2*N);%2n行，32列
R1x=ifft((abs(x1k).^2)/N);%2n行，32列
A=random('unif',0,1,1,32)*2*pi;%32个[0,2*pi]随机相位
for k=1:nx2n(:,k)=cos(2*pi*4*t(:)/N+A(k));%n个相位，随机相位正弦信号xn(:,k)=x1n(:,k)+x2n(:,k);%叠加信号
end
x2k=fft(x2n,2*N);%2n行，32列
R2x=ifft((abs(x2k).^2)/N);%2n行，32列
xk=fft(xn,2*N);%2n行，32列
Rx=ifft((abs(xk).^2)/N);%2n行，32列subplot(321)
plot(m,fftshift(R1x(:,1)));
xlabel('m');ylabel('Rn(m)');
axis([-N N-1 -1 1.5]);
title('一个高斯白噪声自相关')
grid on;subplot(322)
plot(m,fftshift(R1x));
xlabel('m');ylabel('Rn(m)');
axis([-N N-1 -1 1.5]);
title([num2str(n),'个高斯白噪声自相关'])
grid on;subplot(323)
plot(m,fftshift(R2x(:,1)));
xlabel('m');ylabel('Rn(m)');
axis([-N N-1 -1 1.5]);
title('一个随机相位正弦信号自相关')
grid on;subplot(324)
plot(m,fftshift(R2x));
xlabel('m');ylabel('Rn(m)');
axis([-N N-1 -1 1.5]);
title([num2str(n),'个随机相位正弦信号自相关'])
grid on;subplot(325)
plot(m,fftshift(Rx(:,1)));
xlabel('m');ylabel('Rn(m)');
axis([-N N-1 -1 2]);
title('一个叠加信号自相关')
grid on;subplot(326)
plot(m,fftshift(Rx));
xlabel('m');ylabel('Rn(m)');
axis([-N N-1 -1 2]);
title([num2str(n),'个叠加信号自相关'])
grid on;

高斯白噪声仅仅在m=0时取得最大值，其他点很小。随机相位正弦信号的自相关基本符合正弦信号外形。合成信号时的自相关是两个分量自相关的相加。

四、平稳信号的谱估计

1、直接法——周期图法

（1）FFT方法

直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。

当采样频率fs.max大于信号中最高频率fmax的2倍时(fs.max>2fmax)，采样之后的数字信号完整地保留了原始信号中的信息，一般实际应用中保证采样频率为信号最高频率的2.56～4倍；采样定理又称奈奎斯特定理。

clc
clearN=1024;%采样点数
fs=1000;%采样频率，满足采样定理
t=(0:N-1)/fs;%时间，平分了1024份。（0~1.24s）
fai=random('unif',0,1,1,2)*2*pi;%两个相位
A=[1,2];%两个振幅
f=[30,70];%两个频率，可以计算出信号周期
x=A(1)*cos(2*pi*f(1)*t+fai(1))+A(2)*cos(2*pi*f(2)*t+fai(2))+randn(size(t));%该信号
Sx=(abs(fft(x)).^2)/N;%功率谱密度公式f=(0:N/2-1)*fs/N;%频率范围（一半）
subplot(311),plot(f,10*log10(Sx(1:N/2)));%化为dB表示
xlabel('f/Hz');ylabel('Sx(f)/Hz');
grid onf1=(0:N-1)*fs/N;%频率范围（正方向）
subplot(312),plot(f1,10*log10(Sx));%化为dB表示%频率和幅值不对应
xlabel('f/Hz');ylabel('Sx(f)/Hz');
grid onf2=(-N/2:(N-2)/2)*fs/N;%频率范围（正方向）
subplot(313),plot(f2,fftshift(10*log10(Sx)));%化为dB表示%频率和幅值对应
xlabel('f/Hz');ylabel('Sx(f)/Hz');
grid on

close all
clc
clear
N=1024;
fs=1000;
t=(0:N-1)/fs;
fai=random('unif',0,1,1,2)*2*pi;
A=[1,2];
f=[30;70];
x=A*cos(2*pi*f*t+(fai)'*t)+randn(size(t));%利用矩阵乘法，简化了信号叠加的写法
Sx=(abs(fft(x)).^2)/N;
f=(0:N/2-1)*fs/N;
plot(f,10*log10(Sx(1:N/2)));
xlabel('f/Hz');ylabel('Sx(f)/Hz');
grid on

直接用傅里叶变换公式推导的话得到的是双边谱，既-fs——fs,在matlab中用FFT得到的是单边谱，既0——2fs,因为在FFT之后我们得到的是一组复数，需要做abs。乘以2是因为我们的分析频率在0——fs,具体为什么除以点数N，在索引中有详细解释。

（2） periodogram方法

周期图法估计功率谱MATLAB命令为periodogram.

Pxx=periodogram(X)。返回随机序列组成的向量X的功率谱密度，X用与其长度相等的矩形窗。对于实信号，默认返回单边功率谱，复信号，返回双边谱

clear all
clearN=1024;%采样点数
fs=1000;%采样频率，满足采样定理
t=(0:N-1)/fs;%时间，平分了1024份。（0~1.24s）
fai=random('unif',0,1,1,2)*2*pi;%两个相位
A=[1,2];%两个振幅
f=[30,70];%两个频率，可以计算出信号周期
x=A(1)*cos(2*pi*f(1)*t+fai(1))+A(2)*cos(2*pi*f(2)*t+fai(2))+randn(size(t));%该信号
Px=periodogram(x);
f=(0:N/2)*fs/N;%频率范围（一半）
plot(f,10*log10(Px));
xlabel('f/Hz');ylabel('Px(f)/Hz');
title('periodogram法')
grid on

可以看到Px只有513个数据，返回的是单边谱。

clear all
clearN=1024;%采样点数
fs=1000;%采样频率，满足采样定理
t=(0:N-1)/fs;%时间，平分了1024份。（0~1.24s）
fai=random('unif',0,1,1,2)*2*pi;%两个相位
A=[1,2];%两个振幅
f=[30,70];%两个频率，可以计算出信号周期
x=A(1)*cos(2*pi*f(1)*t+fai(1))+A(2)*cos(2*pi*f(2)*t+fai(2))+randn(size(t));%该信号Px=periodogram(x);
f=(0:N/2)*fs/N;%频率范围（一半）
subplot(211),plot(f,10*log10(Px));
xlabel('f/Hz');ylabel('Px(f)/Hz');
title('periodogram法功率谱估计')
grid onSx=(abs(fft(x)).^2)/N;%功率谱密度公式
f=(0:N/2-1)*fs/N;%频率范围（一半）
subplot(212),plot(f,10*log10(Sx(1:N/2)));%化为dB表示
xlabel('f/Hz');ylabel('Sx(f)/Hz');
title('FFT法功率谱估计')
grid on

两种方法计算出来的数值有差异。具体差异可以参见：http://cn.mathworks.com/help/signal/ug/power-spectral-density-estimates-using-fft.html;jsessionid=3e406929b1cfa28e576c6818960d和http://www.ilovematlab.cn/thread-173424-2-1.html

（3）

clear all
clear
%多次频率谱估计平均
N=1024;%采样点数
fs=1000;%采样频率，满足采样定理
M=20;%估计次数
t=(0:N-1)/fs;%时间
fai=random('unif',0,1,1,M)*2*pi;%M个相位
x1=random('norm',0,1,N,M);
A=[1,2];%两个振幅
f=[30,70];%两个频率，可以计算出信号周期for ii=1:Mx(:,ii)=A(1)*cos(2*pi*f(1)*t(:)+fai(ii))+A(2)*cos(2*pi*f(2)*t(:)+fai(ii))+x1(:,ii);%该信号Px(:,ii)=periodogram(x(:,ii));
end
f=(0:N/2)*fs/N;%频率范围（一半）
subplot(121),plot(f,10*log10(Px));
xlabel('f/Hz');ylabel('Px(f)(dB/Hz)');
grid on
title('periodogram法功率谱估计')EPx=mean(Px,2)
subplot(122),plot(f,10*log10(EPx));
xlabel('f/Hz');ylabel('EPx(f)/Hz');
title('periodogram法功率谱估计平均值')
grid on

可以看出提高采样点数，可以提高频率分辨率；利用平均，可以降低功率谱方差。

2、间接法——自相关函数法

clear all
clearN=1024;%采样点数
fs=1000;%采样频率，满足采样定理
t=(0:N-1)/fs;%时间，平分了1024份。（0~1.24s）
fai=random('unif',0,1,1,2)*2*pi;%两个相位
A=[1,2];%两个振幅
f=[30,70];%两个频率，可以计算出信号周期
x=A(1)*cos(2*pi*f(1)*t+fai(1))+A(2)*cos(2*pi*f(2)*t+fai(2))+randn(size(t));%该信号Rx=xcorr(x,'biased');自相关函数
Sx=abs(fft(Rx));
f=(0:N/2)*fs/N;%频率范围（一半）
plot(f,10*log10(Sx(1:N/2+1)));
xlabel('f/Hz');ylabel('Sx(f)(dB/Hz)');
title('自相关法求功率谱(N=1024)')
grid on

可以看到随着N点数的增加，功率谱分辨率得到了提升，但是起伏加剧，也就是分辨率和方差是一对矛盾。

3、改进的功率谱估计——分段平均法和加窗平均法

（1）分段平均（分段加矩形窗）

clear all
clearN=4096;%采样点数
fs=1000;%采样频率，满足采样定理
t=(0:N-1)/fs;%时间，平分了1024份。（0~1.24s）
fai=random('unif',0,1,1,2)*2*pi;%两个相位
A=[1,2];%两个振幅
f=[30,70];%两个频率，可以计算出信号周期
x=A(1)*cos(2*pi*f(1)*t+fai(1))+A(2)*cos(2*pi*f(2)*t+fai(2))+randn(size(t));%该信号
Nseg=1024;
win=rectwin(1024)';%加矩形窗
Sx1=abs(fft(win.*x(1:1024)).^2/norm(win)^2);%2范数，平方求和开根号还是 1024
Sx2=abs(fft(win.*x(1025:2048)).^2/norm(win)^2);%2范数，平方求和开根号还是 1024
Sx3=abs(fft(win.*x(2049:3072)).^2/norm(win)^2);%2范数，平方求和开根号还是 1024
Sx4=abs(fft(win.*x(3073:4096)).^2/norm(win)^2);%2范数，平方求和开根号还是 1024
f=(0:Nseg/2-1)*fs/Nseg;%频率范围（一半）(0~500Hz)
Sx=[Sx1;Sx2;Sx3;Sx4];
for ii=1:4subplot(211),plot(f,10*log10(Sx(ii,1:Nseg/2)));%4段估计hold on
end
legend('Sx1','Sx2','Sx3','Sx4')
xlabel('f/Hz');ylabel('Sx(f)(dB/Hz)');
title('分段估计功率谱')
grid on ESx=(10*log10(Sx1)+10*log10(Sx2)+10*log10(Sx3)+10*log10(Sx4))/4;%求平均
subplot(212),plot(f,ESx(1:Nseg/2));%平均值
xlabel('f/Hz');ylabel('Sx(f)(dB/Hz)');
title('分段平均估计功率谱')
grid on

（2）加汉宁窗

将上面程序矩形窗改为汉宁窗

win=hanning(1024)';%汉宁窗

可以看出，分辨率和方差性能都得到提升。

4、Welch方法

clear all
clearN=4096;%采样点数
fs=1000;%采样频率，满足采样定理
t=(0:N-1)/fs;%时间，平分了1024份。（0~1.24s）
fai=random('unif',0,1,1,2)*2*pi;%两个相位
A=[1,2];%两个振幅
f=[30,70];%两个频率，可以计算出信号周期
x=A(1)*cos(2*pi*f(1)*t+fai(1))+A(2)*cos(2*pi*f(2)*t+fai(2))+randn(size(t));%含噪信号Nseg=1024;%分段长度
win=hanning(1024)';%1024点汉宁窗
noverlap=Nseg/2;%重叠点数
f=(0:Nseg/2)*fs/Nseg;%频率坐标
Sx=pwelch(x,win,256,Nseg,fs,'onesided')*fs/2;%Welch法单边功率谱估计plot(f,10*log10(Sx));%平均值
xlabel('f/Hz');ylabel('Sx(f)(dB/Hz)');
title('Welch法估计功率谱')
grid on

5、歌曲实例

（1）

clear all
clear
[song,fs]=audioread('李健 - 贝加尔湖畔 (自白版).mp3',[1290001,1910000]); % loads
N=length(song(:,1));%采样点数
fs=16000;%采样频率，满足采样定理
for ii=1:2x=song(:,ii);%该信号Sx(:,ii)=(abs(fft(x).^2/N));
end
f=(0:N/2-1)*fs/N;%频率
subplot(211),plot(1:N,song(:,1),1:N,song(:,2))
xlabel('n');ylabel('f(x)');
title('原始双声道数据')
legend('1声道','2声道')
grid on
subplot(212),plot(f,10*log10(Sx(1:N/2,1)),f,10*log10(Sx(1:N/2,2)));
xlabel('f/Hz');ylabel('Px(f)(dB/Hz)');
title('两个声道声道频谱图')
legend('1声道','2声道')
grid onfigure
subplot(221),plot(f,10*log10(Sx(1:N/2,1)));
xlabel('f/Hz');ylabel('Px(f)(dB/Hz)');
title('1声道声道频谱图')
grid on
%Yule-Walker谱估计对象
Hs=spectrum.yulear(16);
subplot(222)
psd(Hs,song(:,1),'fs',16000,'nfft',1024)subplot(223),plot(f,10*log10(Sx(1:N/2,2)));
xlabel('f/Hz');ylabel('Px(f)(dB/Hz)');
title('2声道声道频谱图')
grid on
%Yule-Walker谱估计对象
Hs=spectrum.yulear(16);
subplot(224)
psd(Hs,song(:,2),'fs',16000,'nfft',8096)

两个声道的时域和频域图。可以看出两个声道声音几乎一致。

（2）　　

clear all
clear
clear all
clear
[song1,fs]=audioread('李健 - 贝加尔湖畔 (自白版).mp3',[1290001,1910000]); % loads
% sound(song,fs)
[song2,fs]=audioread('李健 - 贝加尔湖畔 (自白版).mp3',[1910001,2530000]); % loads
% sound(song,fs)
song(:,1)=song1(:,1);
song(:,2)=song2(:,1);
sound(song,fs)

左右耳机出来不一样的音乐。　　

6、典型随机过程

（1）

白噪声：均值为0

平稳白噪声：均值为0，功率谱为定值（所有频率都有，并且值一样），自相关函数只在0处有一个冲击

%平稳高斯白噪声clear all
clear
x=randn(5000,1);
N=length(x);
fs=200;
subplot(221)
plot(x);
xlabel('n'),ylabel('x(n)');
title('平稳高斯白噪声');
grid onsubplot(222)
y=histogram(x,20)
xlabel('x'),ylabel('n/个');
grid onRx1=xcorr(x,'biased');%自相关函数
m=(-N+1):(N-1);%
subplot(223)
plot(m,Rx1);
xlabel('m');ylabel('Rx1');
axis([-N+1 N-1 -1 1.5]);
grid on;
title('平稳高斯白噪声自相关函数');Sx=(abs(fft(x)).^2)/N;%功率谱密度公式
f=(0:N/2-1)*fs/N;%频率范围（一半）
subplot(224),plot(f,Sx(1:N/2));
xlabel('f/Hz');ylabel('Sx(f)/Hz');
grid on

（2）　

clc
clear all
a=0.95;
sigma=2;
N=200;
u=randn(N,1);
x(1)=sigma*u(1)/sqrt(1-a^2);
for ii=2:Nx(ii)=a*x(ii-1)+sigma*u(ii);
end
plot(x);
axis([0,200,-15,15])
xlabel('n')
ylabel('x(n)')
grid on
figure
R=xcorr(x,'coeff');%归一化
plot(R,':r');
hold on
M=1:399;
Rx=sigma^2/(1-a^2)*a.^abs(M-200);
Rx=Rx/((sigma^2)/(1-a^2));%归一化
plot(Rx,'-B')
grid on
legend('随机序列的自相关','理论上的自相关')

转载于:https://www.cnblogs.com/ruo-li-suo-yi/p/7801235.html

MATLAB 随机过程基本理论相关推荐

求解平稳分布matlab,随机过程课程设计--应用马尔科夫链的平稳分布预测空调市场的占有率.doc...
课程名称: <随机过程> 课程设计(论文) 题目: 应用马尔科夫链的平稳分布预测市场占有率学院: 理学院专业: 数学与应用数学班级: 15-1 学生姓名: 邹光睿学 ...
汽车动力性仿真matlab程序,汽车理论课程设计：基于Matlab的汽车动力性的仿真
汽车理论课程设计:基于Matlab的汽车动力性的仿真 2009 届届汽车工程系汽车工程系汽汽车车理理论论课课程程设设计计题题目目汽车动力性的仿真学学院院机电工程学 ...
汽车理论matlab编程,汽车理论课后作业matlab编程详解(带注释)[试题学习]
<汽车理论课后作业matlab编程详解(带注释)[试题学习]>由会员分享,可在线阅读,更多相关<汽车理论课后作业matlab编程详解(带注释)[试题学习](11页珍藏版)>请在 ...
汽车理论matlab编程,汽车理论课后作业matlab编程详解带注释[10页]
<汽车理论课后作业matlab编程详解带注释[10页]>由会员分享,可在线阅读,更多相关<汽车理论课后作业matlab编程详解带注释[10页](11页珍藏版)>请在读根文库上搜 ...
汽车理论课后习题matlab程序,汽车理论课后习题程序
free]<汽车理论>课后部分习题程序[/free] 一.确定一轻型货车的动力性能. 1)绘制汽车驱动力与行驶阻力平衡图: 2)求汽车最高车速与最大爬坡度: 3)绘制汽车行驶加速度倒数曲线 ...
汽车理论matlab编程,汽车理论1.3和2.7matlab编程答案
汽车理论1.3和2.7matlab编程答案汽车理论中期作业1孙野 200812681.3(1)绘制汽车驱动力与行驶阻力平衡图选用 5 挡变速器进行整车性能计算发动机转速与汽车行驶速度之间的关系: 0 ...
汽车理论课后习题matlab程序,汽车理论课后作业matlab编程详解(带注释).doc
汽车理论课后作业matlab编程详解(带注释) 1.3matlab程序: (1)%驱动力-行驶阻力平衡图%货车相关参数. m=3880;g=9.8;nmin=600;nmax=4000;G=m*g;i ...
matlab 摄动波浪理论,基于MATLAB的三维海浪模型数值仿真_齐宁.pdf
ISSN1009-3044 E-mail:eduf@ 第9卷第25期 (2013年09月) ComputerKnowledgeandTechnology电脑知识与技术 ComputerKnowledg ...
最优控制理论与应用matlab,最优控制：理论、方法与应用
评分☆☆☆☆☆ 最优控制是现代控制理论的重要分支,目前已广泛应用于工业生产.经济管理以及国防军事等领域.<最优控制:理论.方法与应用>系统地介绍了最优控制理论内容,包括变分法.极小值原理. ...

MATLAB 随机过程基本理论

MATLAB 随机过程基本理论相关推荐

最新文章

热门文章