随机过程的核心概念与Matlab实现

摘要

本文系统讲解随机过程的核心理论与Matlab实现，涵盖随机变量分布、蒙特卡罗仿真、信息熵计算及平稳过程特性。通过高斯、瑞利分布的生成代码、蒙特卡罗积分估计、窄带信号仿真等案例，结合功率谱分析与自相关函数推导，演示随机过程建模与统计特性验证方法。

关键词：随机过程、Matlab仿真、蒙特卡罗方法、信息熵、平稳过程

1. 随机变量与分布

1.1 概率分布函数与密度函数

随机变量的统计特性通过**累积分布函数（CDF）和概率密度函数（PDF）**描述：

CDF： $\leq x)$
PDF： $\frac{dF(x)}{dx}$

1.2 常见分布生成与Matlab实现

下表列出常见分布的Matlab生成方法：

分布名称	Matlab函数	关键代码片段
均匀分布	`rand()`	`x = rand(1, 1000);`
高斯分布	`randn()`	`x = mu + sigma*randn(1, N);`
瑞利分布	基于高斯变量生成	见代码示例

代码示例：生成瑞利分布随机变量

function s = rayleigh(sigma2, m, n)  x = sqrt(sigma2/2) * randn(m, n);  y = sqrt(sigma2/2) * randn(m, n);  s = sqrt(x.^2 + y.^2);  
end

1.3 二项分布与0-1分布生成

二项分布可通过多个独立0-1分布变量的和生成：

function s = rand01(p, m, n)  x = rand(m, n);  s = (sign(x - p + eps) + 1)/2; % 将均匀变量映射为0-1分布  
end  function s = binom_dist(p, N, m)  y = rand01(p, N, m);  s = sum(y); % 求和生成二项分布  
end

2. 蒙特卡罗仿真方法

2.1 算法原理与积分估计

蒙特卡罗方法通过随机采样近似积分：
$\int_a^b f(x)g(x)dx \approx \frac{1}{N} \sum_{i=1}^N g(x_i)$

示例代码：估计积分 $\int_0^1 x^2 dx$

N = 1e6;  
x = rand(1, N);  
y = mean(x.^2); % 结果接近理论值1/3

2.2 应用示例：通信系统误码率分析

通过蒙特卡罗仿真模拟信号传输过程，统计误码率（BER），代码框架如下：

bits = randi([0,1], 1, N);  
tx_signal = 2*bits - 1; % BPSK调制  
rx_signal = tx_signal + sigma*randn(1, N); % 加高斯噪声  
decoded_bits = (rx_signal > 0);  
BER = sum(bits ~= decoded_bits)/N;

3. 信息论基础

3.1 熵与信息量计算

离散随机变量熵：
$-\sum p(x_i)\log p(x_i)$

代码示例：高斯分布熵的近似计算

x = 1 + randn(1, 1e5);  
px = 1/sqrt(2*pi) * exp(-(x-1).^2 / 2);  
H = -mean(log2(px)); % 结果接近理论值

4. 平稳随机过程

4.1 严平稳与宽平稳特性

严平稳：概率分布与时间起点无关。
宽平稳：均值恒定，自相关函数仅依赖时间差 $\tau$ 。

4.2 自相关函数与功率谱密度

宽平稳过程满足维纳-辛钦定理：
$P(\omega) = \int_{-\infty}^\infty R(\tau)e^{-j\omega\tau}d\tau$

5. 窄带随机过程分析

5.1 等效基带信号生成

窄带过程表达式：
$X_c(t)\cos(2\pi f_c t) - X_s(t)\sin(2\pi f_c t)$

示例：
高斯白噪声是一个功率谱密度为常数，且各时刻满足高斯分布的随机过程。设高斯白噪声的双边功率谱密度为 $N_{0}/2$ ，现将高斯白噪声经过一个中心频率为 $f_{c}$ ，带宽为 $B$ 的带通滤波器，得到的输出信号即为窄带随机过程。
(1) 用Matlab产生高斯平稳窄带随机过程, 设 $f_{c}=10$ ， $B = 1 Hz$ ， $N_{0}=1$
(2)画出窄带高斯过程的等效基带信号；
(3)求出窄带高斯过程的功率，等效基带信号的实部功率、虚部功率。

解：

%例：窄带高斯过程，文件zdpw.m
clear all; close all;
NO = 1; %单边功率谱密度
fc = 10; %中心频率
B = 1; %带宽dt = 0.01;
T = 100;
t = 0:dt:T - dt;
%产生功率为NO*B的高斯白噪声
P = NO*B;
st = sqrt(P)*randn(1,length(t)); 
%将上述白噪声经过窄带带通系统
[f,sf]=T2F(t,st); %高斯信号频谱
figure(1)
plot(f,abs(sf)); %高斯信号的幅频特性[tt,gt]=bpf(f,sf,fc - B/2,fc + B/2);%高斯信号经过带通系统glt = hilbert(real(gt)); %窄带信号的解析信号，调用hibert函数
glt = glt.*exp(-j*2*pi*fc*tt);%得到解析信号[ff,glt]=T2F(tt,glt);
figure(2)
plot(ff,abs(glt));
xlabel('频率(Hz)');ylabel('窄带高斯过程样本的幅频特性')figure(3)
subplot(411);
plot(tt,real(gt));
title('窄带高斯过程样本')subplot(412)
plot(tt,real(glt).*cos(2*pi*fc*tt)-imag(glt).*sin(2*pi*fc*tt))
title('由等效基带重构的窄带高斯过程样本')subplot(413)
plot(tt,real(glt));
title('窄带高斯过程样本的同相分量')subplot(414)
plot(tt,imag(glt));
xlabel('时间(秒)');title('窄带高斯过程样本的正交分量')%求窄带高斯信号功率；注：由于样本的功率近似等于随机过程的功率，因此可能出现一些偏差
P_gt = sum(real(gt).^2)/T;
P_glt_real = sum(real(glt).^2)/T;
P_glt_imag = sum(imag(glt).^2)/T;
%验证窄带高斯过程的同相分量、正交分量的正交性
a = real(glt)*(imag(glt))'/T;function [t,st]= bpf(f,sf,B1,B2)
% This function filter an input at frequency domain by an ideal
% bandpass filter
% Inputs:
%   f:frequency samples
%   sf: input data spectrum samples
%   B1:bandpass's lower frequency 
%   B2:bandpass's higher frequency
% Outputs:
%   st:output data's time samples 
%   t:frequency samplesdf = f(2)-f(1);T = 1/df;hf = zeros(1,length(f));bf= [floor(B1/df):floor(B2/df)];bf1 = floor(length(f)/2)+bf;bf2 = floor(length(f)/2)-bf;hf(bf1)=1/sqrt(2*(B2 - B1));hf(bf2)=1/sqrt(2*(B2 - B1));yf = hf.*sf.*exp(-j*2*pi*f*0.1*T);[t,st]= F2T(f,yf);
endfunction [t,st]=F2T(f,sf)% This function calculate the time signal using ifffunction for the input% signal's spectrumdf=f(2)-f(1);Fmx=(f(end)-f(1)+df);dt=1/Fmx;N = length(sf);T=dt*N;t=-T/2:dt:T/2-dt;%t=0:dt:T-dt;sf = fftshift(sf);st= Fmx* ifft(sf);
endfunction [f,sf]= T2F(t,st)% This is a function using the FFT function to calculate a signal's Fourier% Translation% Input is the time and the signal vectors,the length of time must greater% than 2%Output is the frequency and the signal spectrumdt=t(2)-t(1);T=t(end);df=1/T;N= length(st);f = (-N/2:N/2 - 1) * df;sf = fft(st);sf = T/N* fftshift(sf);
end

核心代码：窄带高斯过程

function zdpw()  N0 = 1; fc = 10; B = 1;  st = sqrt(N0*B) * randn(1, N);  [f, sf] = T2F(t, st); % 傅里叶变换  [tt, gt] = bpf(f, sf, fc-B/2, fc+B/2); % 带通滤波  glt = hilbert(real(gt)) .* exp(-1j*2*pi*fc*t);  
end

5.2 功率计算与正交性验证

P_total = sum(real(gt).^2)/T;  
P_real = sum(real(glt).^2)/T;  
P_imag = sum(imag(glt).^2)/T;

总结

随机过程理论结合Matlab实现，可高效解决通信系统建模、信号分析与信息量化问题。通过代码示例与理论推导，本文展示了从分布生成到复杂系统仿真的全流程方法，为工程实践提供直接参考。