Matlab复现Ungerboeck的TCM经典论文（一）部分公式推导以及M-PAM在AWGN下的信道容量仿真

1.简介

本文主要针对的UNGERBOECK的论文是Channel Coding with Multilevel/Phase Signals ，对于其中的部分公式和部分结果进行了matlab仿真，最近比较忙，博客写的比较潦草，工作量不多，就当辅助理解，抛砖引玉了。

2. 信道容量公式

论文中的公式 (3)，是互信息的变形公式，信息论中的互信息，$I(X;Y)$，其中$X$和$Y$均可以是离散和连续的，在本篇论文中，输入$a$是离散的，输出$z$是连续的。

$$p(z)=\sum _{i=0}^{N-1}Q(i)p(z|a^i)$$
所以根据互信息的定义公式可以被写成下列的形式，并简化
$$\begin{array}{rl}& \sum _{k=0}^{N-1}Q(k){\int }_{-\infty }^{\infty }p(z|a^k){\log }_2\left\{\frac{p(z|a^k)}{{\sum }_{i=0}^{N-1}Q(i)p(z|a^i)}\right\}dz\\ & =\sum _{k=0}^{N-1}Q(k){\int }_{-\infty }^{\infty }p(z|a^k)\left\{{\log }_{2}p(z|a^k)-{\log }_2\left[\sum _{i=0}^{N-1}Q(i)p(z|a^i))\right]\right\}dz\\ & =\sum _{k=0}^{N-1}Q(k){\int }_{-\infty }^{\infty }p(z|a^k)\left\{{\log }_{2}p(z|a^k)-{\log }_{2}p(z)\right\}dz\\ & =\sum _{k=0}^{N-1}{\int }_{-\infty }^{\infty }Q(k)p(z|a^k){\log }_{2}p(z|a^k)dz-\sum _{k=0}^{N-1}{\int }_{-\infty }^{\infty }Q(k)p(z|a^k){\log }_{2}p(z)dz\\ & =\sum _{k=0}^{N-1}{\int }_{-\infty }^{\infty }Q(k)p(z|a^k){\log }_{2}p(z|a^k)dz-{\int }_{-\infty }^{\infty }p(z){\log }_{2}p(z)dz\end{array}$$
其中，$Q(k)p(z|a^k)$可以看做是联合分布，所以我们有
$$\begin{array}{rl}& \sum _{k=0}^{N-1}{\int }_{-\infty }^{\infty }Q(k)p(z|a^k){\log }_{2}p(z|a^k)dz-{\int }_{-\infty }^{\infty }p(z){\log }_{2}p(z)dz\\ & =h(Z)-h(Z|A)\\ & =I(Z;A)\end{array}$$
因此我们可计算信道容量，和论文的公式 (3)一致
$$\begin{array}{rl}C& =\underset{Q(0)\dots Q(N-1)}{\max }I(Z;A)\\ & =\underset{Q(0)\dots Q(N-1)}{\max }\sum _{k=0}^{N-1}Q(k){\int }_{-\infty }^{\infty }p(z|a^k){\log }_2\left\{\frac{p(z|a^k)}{{\sum }_{i=0}^{N-1}Q(i)p(z|a^i)}\right\}dz\end{array}$$

3. $Q(k)=\frac{1}{N}$下的信道容量公式

接下来是公式(5)的推导
$$\begin{array}{rl}{C}_{Q(k)=1/N}^{\ast }& =\sum _{k=0}^{N-1}\frac{1}{N}{\int }_{-\infty }^{\infty }p(z|a^k){\log }_2\left\{\frac{p(z|a^k)}{{\sum }_{i=0}^{N-1}\frac{1}{N}p(z|a^i)}\right\}dz\\ & =\sum _{k=0}^{N-1}\frac{1}{N}{\int }_{-\infty }^{\infty }p(z|a^k)\left\{-{\log }_2\left[\sum _{i=0}^{N-1}\frac{1}{N}\frac{p(z|a^i)}{p(z|a^k)}\right]\right\}dz\\ & =\sum _{k=0}^{N-1}\frac{1}{N}{\int }_{-\infty }^{\infty }p(z|a^k)\left\{{\log }_2(N)-{\log }_2\frac{{\sum }_{i=0}^{N-1}p(z|a^i)}{p(z|a^k)}\right\}dz\\ & ={\log }_2(N)\sum _{k=0}^{N-1}{\int }_{-\infty }^{\infty }\frac{1}{N}p(z|a^k)dz-\frac{1}{N}\sum _{k=0}^{N-1}{\int }_{-\infty }^{\infty }p(z|a^k){\log }_2\frac{{\sum }_{i=0}^{N-1}p(z|a^i)}{p(z|a^k)}dz\\ & ={\log }_2(N)-\frac{1}{N}\sum _{k=0}^{N-1}{\int }_{-\infty }^{\infty }p(z|a^k){\log }_2\sum _{i=0}^{N-1}\frac{p(z|a^i)}{p(z|a^k)}dz\end{array}$$
论文中的公式(4)如下
$$p(z|a^k)=\mathrm{e}\mathrm{x}\mathrm{p}[-|z-a^k|/2{\sigma }^2]\cdot \{\begin{array}{rl}& (2\pi {\sigma }^2{)}^{-1/2}one-dimension\\ & (2\pi {\sigma }^2{)}^{-1}two-dimension\end{array}$$

将公式(4)代入，并且我们采用蒙特卡洛法进行模拟计算
$$\begin{array}{rl}& {\log }_2(N)-\frac{1}{N}\sum _{k=0}^{N-1}{\int }_{-\infty }^{\infty }p(z|a^k){\log }_2\sum _{i=0}^{N-1}\frac{p(z|a^i)}{p(z|a^k)}dz\\ & ={\log }_2(N)-\frac{1}{N}\sum _{k=0}^{N-1}E\left\{{\log }_2\sum _{i=0}^{N-1}\frac{p(z|a^i)}{p(z|a^k)}\right\}\\ & ={\log }_2(N)-\frac{1}{N}\sum _{k=0}^{N-1}E\left\{{\log }_2\sum _{i=0}^{N-1}\mathrm{e}\mathrm{x}\mathrm{p}\left[\frac{-|z-a^i{|}^2+|z-a^k{|}^2}{2{\sigma }^2}\right]\right\}\\ & ={\log }_2(N)-\frac{1}{N}\sum _{k=0}^{N-1}E\left\{{\log }_2\sum _{i=0}^{N-1}\mathrm{e}\mathrm{x}\mathrm{p}\left[\frac{-|a^k+w-a^i{|}^2+|w{|}^2}{2{\sigma }^2}\right]\right\}\end{array}$$
当$k$确定时，其中$z=a^k+w$

4. 使用蒙特卡洛法仿真M-PAM在AWGN下的信道容量

直接贴matlab代码了，这个就是论文中的Fig.2

clear
clc

% M-PAM
m_arr = [2,4,8,16,32,64];

%SNR
snr_dB=[-15:0.5:50];

%capacity for differenet SNR and M
c_plot = zeros(length(m_arr),length(snr_dB));

%iterations per M
iters=10; 

%symbols per iteration
sym_num=1000;


% for each M
for m_index = 1:length(m_arr)
    m = m_arr(m_index);
    m

    %assume E{|a|^2} =1
    d = sqrt(12/(m^2-1)); 
    constell = (1:m)*d;

    %generate constellation diagram
    constell = constell - mean(constell);
    %uniform distribution
    prop_constell=ones(1,length(constell)) * 1/ length(constell);
    power=sum(prop_constell.*constell.^2);  % =1
    for snr_index=1:length(snr_dB) 
        snr=10^(snr_dB(snr_index)/10); 
        %sigma of N
        sigma2=power/snr;  
        sigma=sqrt(sigma2);
        c_temp=zeros(1,iters);    
        for index_iters=1:iters
            %a before modualtion
            x=rand(1,sym_num);           
            pconstell=cumsum(prop_constell);
            a=zeros(1,sym_num);  
            a(x<=pconstell(1))=constell(1);     
            %  Modulation 
            for index_constell=2:length(constell)
                a(x>pconstell(index_constell-1)&x<=pconstell(index_constell))=constell(index_constell);    
            end           
            % Noise (0,sigma2)
            Noise=sigma*randn(1,sym_num); 
            % z
            z = a+Noise;
%             equation 5, the plot result is the same but slower while simulation
%             sum1=0;
%             for k =1:m
%                sum2=0;
%                for i = 1:m
%                    exp1 = exp(-(Noise).^2/(2*sigma2));
%                    inner = (Noise+constell(k)-constell(i));
%                    exp2 =exp((-(inner).^2)/(2*sigma2));
%                    sum2 = sum2+(exp2./exp1);
%                end
%                sum1 = sum1+ mean(log2(sum2));
%             end      
%             c_temp(index_iters)=log2(m) -(sum1/m);    
        
            %also equation 5 but much faster while simulation
            sum_tmp=0;
            for i = 1:m
                p1=exp(-(Noise).^2/(2*sigma2));
                p2 =exp(-(z-constell(i)).^2/(2*sigma2));
               sum_tmp = sum_tmp+ p2./p1;
            end
            c_temp(index_iters)=log2(m)-mean(log2(sum_tmp));         
        end  
        c(snr_index)=mean(c_temp);
    
    end
    %save the result
    c_plot(m_index,:) = c(:);
end

%ideal capacity
SNR = zeros(1,length(snr_dB));
targetC =zeros(1,length(snr_dB));

new_snr_dB=[-15:0.01:50];
for i =1:length(new_snr_dB)
    SNR(i) = 10^(new_snr_dB(i)/10);
    targetC(i) = 0.5*log2(1+SNR(i));
end

for i = 1:length(m_arr)
    plot(snr_dB,c_plot(i,:));
    hold on;
end
plot(new_snr_dB,targetC,'r');
grid on
legend(['Uniform ' num2str(m) '-constell']);
temp = c_plot(6,:);
xlabel('Signal to noise ratio [dB]');
ylabel('Capacity [b/dim]');
legend('2-PAM','4-PAM','8-PAM','16-PAM','32-PAM','64-PAM','1/2log2(1+SNR)')

画出来的图像应该如下所示

5. 误码率公式以及gap

在M-PAM的星座图中, 我们有$M$个星座点$p_0\cdots p_{M-1}$. 接收到的symbol $s$在星座图上，存在两种形式. 它存在于最外层的两个星座点，或者属于内部的$M-2$个星座点. 假设$d$是symbol $s$距离它的正确点$p^{\prime }$ 的距离。并且假设$D$是星座点之间的距离。
$$\begin{array}{rl}{P}_s& =\mathrm{Pr}(s\ne p^{\prime })\\ & =\frac{2}{M}\mathrm{Pr}(d<\frac{D}{2})+\frac{M-2}{M}\left[\mathrm{Pr}(d<-\frac{D}{2})+\mathrm{Pr}(d>\frac{D}{2})\right]\\ & =\frac{2}{M}\mathrm{Pr}(\frac{d}{\sqrt{N_0/2}}<\frac{D}{2\sqrt{N_0/2}})+\frac{2(M-2)}{M}\mathrm{Pr}(\frac{d}{\sqrt{N_0/2}}<\frac{D}{2\sqrt{N_0/2}})\\ & =\frac{2(M-1)}{M}Q(\sqrt{\frac{D}{2N_0}})\end{array}$$
我们又有关于$SNR$和$E_S$的公式如下
$$\begin{array}{rl}SNR& =\frac{E_s}{N_0/2}\\ E_s& =(M^2-1)D/12\end{array}$$
所以误码率如下
$$\begin{array}{r}{P}_s=\frac{2(M-1)}{M}Q(\sqrt{\frac{3SNR}{M^2-1}})\end{array}$$

由此，我们发现，在误码率$P_e=10^{-5}$下，存在一个gap损失，如下图

我们接下来就需要TCM编码来缩小这一Gap的损失

6. 总结

这是在国外选修的第一门课的其中一次作业，所以代码注释是英语的，英文不好，见谅。马上圣诞节了，无心学习，所以写的很简略，见谅见谅。