第15章 离散信道模型 MATLAB 仿真笔记

本篇文章是《通信系统仿真原理与无线应用》这本书的学习笔记，仅作个人学习使用。

对于有记忆信道，最常用的模型是离散时间的有限状态马尔可夫模型。马尔可夫模型得到广泛应用有几方面的原因。

这类模型易干分析处理，在统计方面的文献中已建立了完备的理论，并且它们已经被成功地应用千解决各种类型的重要通信问题。马尔可夫模型一直被用来建立英文文本之类的离散信源模型，（在一段典型的英文文章中，出现的一列英语字母可以用一个马尔可夫序列来建模）。

15.3.1 两状态模型

为了准备讨论马尔可夫通道模型，先来考虑一个衰落信道。在这个衰落信道中，有一部分时间内，接收信号的强度在可接受的性能阖值之上，而在深度衰落时会处千阖值之下。如果我们只对“在阖值之上”还是“在阖值之下” 的情况感兴趣，那么我们可以将信道建模为处千如图15-5所示的两个状态之一。在图15-5 中，有一个好状态 $g$, 其系统性能为可接受的状态（$P_E<10^{-3}$ ），另一个是坏状态 $b$ , 其接收信号很弱，以至千差错概率高得不可接受（即 $P_E>10^{-3}$）。

因此，在如图15-5所示的简单两状态模型中，状态可用如下集合表示

S = { g , b } (5) S=\{g , b\}\tag5 S={g,b}(5)

随着时间的推进，信道可能会从好状态转移到坏状态的， 也可能从坏状态转移到好状态， 转移的速率和在某一个状态的驻留时间都将取决于衰落过程的时间相关性。

令 $S_{t+1}$ 和 $S_t$ ，分别表示信道在时刻 t+1 和时刻 t 时的状态。（ 注意 t + 1 表示比 t 大一个时间增量的时刻，而不是比 t 大 1s 的时刻，所以 t 可看作离散时间序数。）使用这个记号，我们定义四个转移概率为

a g g ( t ) = Pr ⁡ { S t + 1 = g ∣ S t = g } a g b ( t ) = Pr ⁡ { S t + 1 = b ∣ S t = g } a b g ( t ) = Pr ⁡ { S t + 1 = g ∣ S t = b } a b b ( t ) = Pr ⁡ { S t + 1 = b ∣ S t = b } \begin{array}{c} a_{g g}(t)=\operatorname{Pr}\left\{S_{t+1}=g \mid S_{t}=g\right\} \\ a_{g b}(t)=\operatorname{Pr}\left\{S_{t+1}=b \mid S_{t}=g\right\} \\ a_{b g}(t)=\operatorname{Pr}\left\{S_{t+1}=g \mid S_{t}=b\right\} \\ a_{b b}(t)=\operatorname{Pr}\left\{S_{t+1}=b \mid S_{t}=b\right\} \end{array} agg​(t)=Pr{St+1​=g∣St​=g}agb​(t)=Pr{St+1​=b∣St​=g}abg​(t)=Pr{St+1​=g∣St​=b}abb​(t)=Pr{St+1​=b∣St​=b}​

这些转移概率可以用如下状态转移矩阵表示

$$A(t)=\left[\begin{array}{ll}{a}_{gg}(t)& a_{gb}(t)\\ a_{bg}(t)& a_{bb}(t)\end{array}\right]$$

注意，由千我们考虑的是两状态模型，给定状态 $S_t=g$ 必定有 $S_{t+1}=g$ （信道保持在好状态）或 $S_{t+1}=b$（信道转移到坏状态），因此，状态转移矩阵每一行的和都必为1 。两状态模型的状态转移图如图15-6所示。

在下面的部分，我们将假定信道模型是静态的，因此状态转移矩阵 $A(t)$ 是固定的，即 $A(t)=A$。然而， 如果仿真运行中初始状态概率分布 ${\Pi }_0$ 不同于稳态分布，那就需要一段时间来使概率分布演化到稳态值 ${\Pi }_{ss}$。我们关心的是模型处于给定状态下的概率，定义 ${\Pi }_t$，为 $t$ 时刻的状态概率分布，具体而言，有
$${\Pi }_t=\left[\begin{array}{ll}{\pi }_{t,g}& {\pi }_{t,b}\end{array}\right]$$

式中，${\pi }_{t,g}$ 和 ${\pi }_{t,b}$ 分别表示在 $t$ 时刻信道处于好状态或坏状态的概率。通过定义状态转移矩阵，在时刻 $t+1$ 的状态分布为
$${\Pi }_{t+1}={\Pi }_{t}A\text{\hspace{1em}(9)}$$

类似地，在时刻 $t+2$ 时的状态分布是

$${\Pi }_{t+2}={\Pi }_{t+1}A=\left({\Pi }_{t}A\right)A={\Pi }_t{A}^2$$

因此， 一般地有

$${\Pi }_{t+k}={\Pi }_t{A}^k\text{\hspace{1em}(11)}$$

式中， $A^k$ 代表 $k$ 步转移矩阵。

大多数情况下（但不是所有情况下），随时间的推进，马尔可夫过程都会演化到一个稳态概率分布。假定对干足够大的时间 $t$ 和任意的 $k$, 马尔可夫过程收敛到一个稳态分布 ${\Pi }_{t+k}={\Pi }_t$,则对任意 $k$, 有下式

$${\Pi }_{ss}={\Pi }_{ss}{A}^k=\left[\begin{array}{ll}{\pi }_g& {\pi }_b\end{array}\right]\text{\hspace{1em}(12)}$$

例15-1 作为一个简单的例子，令
$$A=\left[\begin{array}{ll}0.98& 0.02\\ 0.05& 0.95\end{array}\right]$$
初始状态分布为
$${\Pi }_0=\left[\begin{array}{ll}0.50& 0.50\end{array}\right]$$
下面的MATLAB程序说明了从 ${\Pi }_0$ 到 ${\Pi }_{ss}$ 的收敛过程：

% File: c15_MMtransient.m
% Software given here is to accompany the textbook: W.H. Tranter, 
% K.S. Shanmugan, T.S. Rappaport, and K.S. Kosbar, Principles of 
% Communication Systems Simulation with Wireless Applications, 
% Prentice Hall PTR, 2004.
%
N = 100;
pie = zeros(N,2);
A = [0.98 0.02; 0.05 0.95];
pie(1,:) = [0.50 0.50];
for k=2:N
   pie(k,:) = pie(k-1,:)*A;
end
kk = 1:N;
plot(kk,pie(:,1),'k-',kk,pie(:,2),'k:')
xlabel('Iteration')
ylabel('Probability')
text1 = ['The steady-state probabilities are ',...
      num2str(pie(N,1)),' and ' ,num2str(pie(N,2)),'.'];
legend('State 1','State 2')
disp(text1)
disp(' ')
disp('The value of A^N is'); A^N
% End of script file.

执行程序得到如下结果

>> c15_MMtransient
The steady-state probabilities are 0.71412 and 0.28588.

The value of A^N is

ans =

0.7145    0.2855
0.7138    0.2862

上述代码较为简单，基本没有什么好解释的，唯一需要注意这里使用了矩阵pie存储马尔可夫模型的状态$\left[\begin{array}{ll}{\pi }_{t,g}& {\pi }_{t,b}\end{array}\right]$。最后使用了语句plot(kk,pie(:,1),'k-',kk,pie(:,2),'k:')来绘制逐渐收敛的信道状态。

15.3.4 平稳性

因为在建模倌道模拟部分时一般会假定平稳性，因此在离散信道的马尔可夫模型中常常也假定平稳性。平稳性意味着信道模型的参数，即概率${\Pi }_{t,i}$、$a_{ij}(t)$和$b_i(e_k)$ 跟时间无关。在这种情况下，我们有

$$\mathrm{Pr}\left[S_{t+1}=i\right]={\Pi }_i=\sum _{k=1}^N\mathrm{Pr}\left[S_t=k\right]\mathrm{Pr}\left[S_{t+1}=i\mid S_t=k\right]=\sum _{k=1}^N{\pi }_k{a}_{ki},i=1,\cdots ,N$$
式中 $a_{ki}$ 是状态转移矩阵第 $i$ 列的元素

$$A=\left[\begin{array}{cccccc}{a}_{11}& a_{12}& \cdots & a_{1i}& \cdots & a_{1N}\\ a_{21}& a_{22}& \cdots & a_{2i}& \cdots & a_{2N}\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮\\ a_{N1}& a_{N2}& \cdots & a_{Ni}& \cdots & a_{NN}\end{array}\right]$$

而 ${\Pi }_k$ 是状态概率分布向量的元素
$$\Pi =\left[\begin{array}{llll}{\pi }_1& {\pi }_2& \cdots & {\pi }_N\end{array}\right]$$

跟两状态模型的情况一样， n 步转移矩阵由 $A^n$ 给出。正如例15-1 对两状态模型所作的阐述，式 (15-34) 及其约束

$$\sum _{i=1}^N{\pi }_i=1$$

意味状态概率 ${\pi }_i$ 可由转移概率 $a_{ij}$ 唯一确定。因此马尔可夫模型完全由状态转移概率矩阵 $A$ 和 $B$ 确定，其中矩阵 $B$ 是前面所讨论的输人－输出符号转移（即差错）概率：

$$\mathbf{B}=\left[\begin{array}{ccc}{b}_{11}& \cdots & b_{1N}\\ ⋮& \ddots & ⋮\\ b_{M1}& \cdots & b_{MN}\end{array}\right]$$

15.3.5 马尔可夫模型仿真

一旦推导出了离散马尔可夫模型，对模型进行仿真就相对容易。

给定模型处在状态 $k$, 转移概率由状态转移矩阵 $A$ 的第 $k$ 行定义，矩阵 $A$ 的第 $k$ 行可用向量 $A_k$ 给
出，其中

$$A_k=\left[\begin{array}{lllllll}{a}_{k1}& a_{k2}& \cdots & a_{ki}& a_{k,i+1}& \cdots & a_{kN}\end{array}\right]$$

$A_k$ 的累积概率记作 ${\beta }_k$， 由下式给出

$${\beta }_k=\sum _{j=1}^k{a}_{kj}$$

注意， MATLAB 的库函数cumsum将向量 $A_k$ 映射到具有元素 ${\beta }_k$ 的向量。定义从状态 $k$ 转移到状
态 $i$ 的概率为 $a_{ki}$ ，其值由下式给出：

$$a_{ki}={\beta }_i-{\beta }_{i-1}$$

这个概率如图15-8(a) 中的阴影面积所示。

在确定了新的状态之后，下一步是决定在新状态是否有差错发生。为了完成这项工作，我们将第二个随机变量 $U_2$ 与阈值 $\gamma$ 比较，如图15-8(b) 所示。正确判决的概率是 P r { U 2 < γ } Pr\{ U_2<\gamma\} Pr{U2​<γ} ，而错误判决的概率是 P r { U 2 > γ } Pr\{ U_2>\gamma\} Pr{U2​>γ}，如图15-8(b) 的阴影区域所示。由于当前状态是 $i$, $\gamma$ 的值由状态观察矩阵 $B$ 第 1 行第 $i$ 列的元素 $b_{1,i}$ 给出。

实现仿真过程的MATLAB代码如下：

%A是状态转移矩阵
u1=rand(1);
cum_sum=[0 cumsum(A(state,:))];
for i=i:total_states
    if u1>cum_sum(i)&u1<cum_sum(i+1);
        state=i;
    end
end

state_seq(t)=state;
%B是符号差错概率矩阵
u2=rand(1);
if u2>B(1,state)
    out(t)=1;
end

例15-3 在这个例子中， 产生长为 $N=20000$ 的二进制序列来表示信道上的误码，用一个三状态马尔可夫模型来产生这个差错序列。如果序列中给定元素是二进制 ”1" ，则在给定元素对应的位置记录一个信道误码，二进制 “0” 则表示对应的位置没有发生差错。

我们假定三个状态中的任意一个都可能发生差错。我们具体地作以下假定：

Pr ⁡ { E ∣ S 1 } = 0.0010 Pr ⁡ { E ∣ S 2 } = 0.0500 Pr ⁡ { E ∣ S 3 } = 0.0100 \begin{array}{l} \operatorname{Pr}\left\{E \mid S_{1}\right\}=0.0010 \\ \operatorname{Pr}\left\{E \mid S_{2}\right\}=0.0500\\ \operatorname{Pr}\left\{E \mid S_{3}\right\}=0.0100 \end{array} Pr{E∣S1​}=0.0010Pr{E∣S2​}=0.0500Pr{E∣S3​}=0.0100​

它对应了差错概率矩阵，其中 $B$ 的第 1 列即为状态 1

$$B=\left[\begin{array}{lll}0.9990& 0.9500& 0.9900\\ 0.0010& 0.0500& 0.0100\end{array}\right]$$

假设马尔可夫模型的状态转移矩阵为
$$A=\left[\begin{array}{lll}0.80& 0.10& 0.10\\ 0.20& 0.60& 0.20\\ 0.02& 0.08& 0.90\end{array}\right]$$

我们假设模型最初处千状态 1 。

% File: c15_errvector.m
% Software given here is to accompany the textbook: W.H. Tranter, 
% K.S. Shanmugan, T.S. Rappaport, and K.S. Kosbar, Principles of 
% Communication Systems Simulation with Wireless Applications, 
% Prentice Hall PTR, 2004.
%
disp(' ')
disp('Default values are:')
N = 20000											% default N			
A = [0.8 0.1 0.1; 0.2 0.6 0.2; 0.02 0.08 0.90]		% default A
B = [0.999 0.95 0.99; 0.001 0.05 0.01]				% default B
disp(' ')
disp('Accept default values?')
dtype = input('Enter y for yes or n for no  > ','s');
if dtype == 'n'
   N = input('   Enter N, the number of points to be generated > ');
   A = input('   Enter A, the state transition matrix > ');
   B = input('   Enter B, the error distribution matrx > ');
end
state = 1;							            % initial state
total_states = size(A,1);
out = zeros(1,N);					            % initialize error vector
state_seq = zeros(1,N);			                % initialize state sequence
h = waitbar(0,'Calculating Error Vector');
%

u2 = rand(1);	 					            % get random number
if u2>B(1,state)					            % test for error
   out(1) = 1;						            % record error
end
state_seq(1) = state;			                % record state
for t=2:N
   u1 = rand(1);					            % get random number
   cum_sum = [0 cumsum(A(state,:))];
   for i=1:total_states			                % loop to determine new state
      if u1>=cum_sum(i) & u1<cum_sum(i+1);
         state = i;				                % assign new state
      end
   end	
   state_seq(t) = state;		                % new record state
   u2 = rand(1);					            % get random number
   if u2>B(1,state)
      out(t) = 1;					            % record error
   end
   waitbar(t/N)
end
close(h)
% End of script file.

h = waitbar(0,'Calculating Error Vector');

生成一个进度条窗口，0 表示刚开始时候，进度条上显示英文 Calculating Error Vector。

这行代码创建了一个进度条（progress bar），用于显示代码运行的进度。它的第一个参数是当前进度，这里是0，第二个参数是进度条显示的文本。在这里，文本是“Calculating Error Vector”。在代码中的循环中，通过调用 waitbar 函数并将当前进度作为参数传递给它来更新进度条的进度。当循环结束时，通过调用 close 函数关闭进度条。这样可以方便地显示代码的运行进度，并且用户可以随时中止程序的运行。

在 waitbar 函数中，第一个参数是一个标量，表示当前的进度百分比，取值范围是 0 到 1。在这里，初始值是 0，因为代码还没有开始运行。随着代码的执行，循环中的语句逐步执行，waitbar 函数也会被调用多次，每次更新进度条的进度。waitbar(t/N)在每次调用 waitbar 函数时，第一个参数都会更新为当前的进度百分比。最终，当代码执行完成后，进度条的进度会更新为 100%，close(h)。因此，0 的用法是为了初始化进度条，表示代码还没有开始执行，进度为 0。

附录 A 给出了用来产生差错序列的 MATLAB 程序c15_errvector.m 。程序产生两个输出，一个是差错向量 out，它给出了 N 个传输序列中的差错位置，另一个是状态序列向量 state_seq, 它给出了状态转移的细节。模型处千给定状态下的概率可以由状态序列向量计算出来，仿真差错概率可以从差错向量计算出来。

% File: c15_hmmtest.m
% Software given here is to accompany the textbook: W.H. Tranter, 
% K.S. Shanmugan, T.S. Rappaport, and K.S. Kosbar, Principles of 
% Communication Systems Simulation with Wireless Applications, 
% Prentice Hall PTR, 2004.
%
pe = sum(out)/N;%pe表示误码率
state_sum = zeros(1,total_states);
for k=1:N
   if state_seq(k)==1
      state_sum(1)=state_sum(1)+1;
   end   
   if state_seq(k)==2
      state_sum(2)=state_sum(2)+1;
   end
   if state_seq(k)==3
      state_sum(3)=state_sum(3)+1;
   end
end   
a = ['The probability of State 1 is ',num2str(state_sum(1)/N),'.'];
b = ['The probability of State 2 is ',num2str(state_sum(2)/N),'.'];
c = ['The probability of State 3 is ',num2str(state_sum(3)/N),'.'];
d = ['The error probability is ',num2str(pe),'.'];
disp('Simulation results:')
disp(a)			% display probability of state 1
disp(b)			% display probability of state 2
disp(c)			% display probability of state 3
disp(d)			% display error probability
% End script file.   

总结一下，c15_errvector.m 实现了对差错向量和状态序列向量的仿真，c15_hmmtest t.m 实现了对差错概率以及每种状态概率的计算，如图。