通信系统仿真速成第5天：PAM系统在AWGN信道下的互信息

信息论也好，通信原理也好，有个互信息的概念。

互信息(Mutual Information)是信息论里一种有用的信息度量。

它可以看成是一个随机变量中包含的关于另一个随机变量的信息量，

或者说是一个随机变量由于已知另一个随机变量而减少的不肯定性。

百度百科：互信息

数字通信原理 Digital Communication (作者：Proakis)

还有一个东西叫AWGN下的信道容量。当初上课的时候，我确实没有想过，书上这种图是怎么画出来的……

后来写了些代码画了一些图。这里放上求PAM的互信息的代码。

基本想法就是，不断地发，然后统计误符号率，估算系统传递矩阵，根据公式求出I（X;Y）。

百度百科：PAM（脉冲幅度调制）

%% PAM-m的AIR% 内容：希望作一条AIR v.s Eb/N0 的curve% 作者：qcy% 版本：v1.5% 说明：希望多次循环，求I的平均%%clear;close all;clc%%d = 1; % 固定星座点的距离mod_order = 4;x_table = (1:mod_order).' * d;x_table = x_table - mean(x_table);Nsym = 1e5; % 一次循环中，仿多少个符号% rand('seed',1);iMsg = randi(mod_order,Nsym,1);x_data = x_table(iMsg);% scatterplot(x);% grid on;% title('PAM4');N_iter = 2e1; % 多次循环，求I的平均%%Eb_N0_dB_array = (2:.5:12).';I_array = zeros(size(Eb_N0_dB_array));for k_Eb_N0_dB_array = 1:length(Eb_N0_dB_array)        Eb_N0_dB = Eb_N0_dB_array(k_Eb_N0_dB_array);    px = ones(length(x_table),1) * 1/ length(x_table); % 等概率分布        for k_iter = 1:N_iter        x_noise = add_awgn_pam_1(x_data,x_table,Eb_N0_dB,px);        msg_decision = pam_decision_1( x_noise,x_table );        system_property = calc_system_property_1( iMsg , msg_decision , mod_order );        I_array(k_Eb_N0_dB_array) = I_array(k_Eb_N0_dB_array) + system_property.IXY;    end        I_array(k_Eb_N0_dB_array) = I_array(k_Eb_N0_dB_array) / N_iter;end%%figure(1);clf;plot(Eb_N0_dB_array,I_array);grid on;axis([Eb_N0_dB_array(1) Eb_N0_dB_array(end) 1 log2(mod_order)*1.1]);legend(['Uniform PAM-' num2str(mod_order)]);xlabel('E_b/N_0 [dB]');ylabel('I(X;Y) [bit/symbol]');title(['Uniform PAM-' num2str(mod_order)]);

function [ x_noise,noise ] = add_awgn_pam_1( x_data,x_table,Eb_N0_dB,px )%function [ x_noise,noise ] = add_awgn_pam_1( x_data,x_table,Eb_N0_dB,px )%   向1-D信号加入awgn%   输入参数%       x_data：基带信息 int型msg%       x_table：星座图映射表%       Eb_N0_dB：用dB表示的比特信噪比%   返回参数%       x_noise：带噪基带信号%       noise：噪声mod_order = length(x_table);if nargin < 3 % 如果没有指定分布，就默认先验等概   px = ones(mod_order,1) * 1/ mod_order;endNsym = length(x_data);Es = sum((abs(x_table).^2).*px);Eb = Es / log2(mod_order); % 【问】 每个比特携带的能量% 【注意】好像是这样。。。（你们有兴趣可以再检查一下）% 是除以-_- 不是乘以% e.g. 2个比特携带Es的能量，那么1个比特携带Eb = Es/2的能量;Eb_N0_linear = 10 ^ (Eb_N0_dB/10);N0 = Eb / Eb_N0_linear;% N0_2 = N0/2;sigma2  = N0/2; % sigma^2;sigma = sqrt(sigma2);noise = sigma * randn(Nsym,1);x_noise = x_data + noise;end

function [ msg_decision ] = pam_decision_1( x_input,x_table )%function [ msg_decision ] = pam_decision_1( x_input,x_table )%   1-D信号判决函数。%   输入参数%       x_input：带判决%       x_table：星座映射表%   返回参数%       msg_decision：判决成星座表中的下标idx%%   作者：qcy%   版本：v1.0msg_decision = zeros(size(x_input));for m = 1:length(x_input)    x_temp = x_input(m);    dist_min = inf;    k_min = 0;    for n = 1 : length(x_table)        dist = abs(x_temp - x_table(n));        if dist < dist_min            dist_min = dist;            k_min = n;        end    end    msg_decision(m) = k_min;endend

function [ system_property ] = calc_system_property_1( iMsg,msg_decision,mod_order )%function [ system_property ] = calc_system_property_1( iMsg,msg_decision，mod_order )% 计算任意系统特性（信源熵、传输概率矩阵）% [注意]% 此时的信源熵，是根据iMsg中，各信源符号出现的频率代替概率来计算的。% 并不是根据信源的分布律来算的。% 此时算出来的先验概率，其实也只是频率。% 转移矩阵的各个元素，其实也只是频率。% 依据：频率是可以趋近于概率的。%%   输入参数%       iMsg：原始的发送符号（1、2、3、4、……）%       msg_decision：接收判决后的符号（1、2、3、4、……）%       mod_order：调制阶数%   返回参数%       system_property：能零误码达到的最大信息速率（bit/symbol）%       system_property.TP：转移概率矩阵%       system_property.Px：信源先验概率%       system_property.IXY：系统平均互信息%%   作者：qcy%   版本：v1.0Nsym = length(iMsg);% 初始化structuresystem_property.TP = zeros(mod_order,mod_order); % 转移概率矩阵system_property.IXY = 0; % 互信息system_property.Hx = 0; % 信源熵for m = 1:mod_order    for n = 1:mod_order        send_x_idx = (iMsg==m); % 发m，那一位为1，否则为0        decision_y_idx = (msg_decision==n); % 判n，那一位为1，否则为0        temp = (send_x_idx + decision_y_idx); % 加起来 都是 2 ，那一位才代表 发1，判1        f_num = length(find(temp==2));        system_property.TP(m,n) = (f_num / (nnz(iMsg == m)+eps)); % 发m，判n的概率 -_-?    endend% system_property.Psystem_property.Px = zeros(mod_order,1);for k = 1:mod_order    system_property.Px(k) = nnz(iMsg == k) / Nsym;end% system_property.HxHx = -sum(system_property.Px .* log2(system_property.Px+eps));system_property.Hx = Hx;I_part1 = Hx;I_part2 = 0;mod_order = length(system_property.Px);for m = 1:mod_order % i    for n = 1:mod_order % j        pxi = system_property.Px(m);        pyjxi = system_property.TP(m,n);                xxx1 = 0; % log里的分子        for k = 1:mod_order            pxk = system_property.Px(k);            pyjxk = system_property.TP(m,k);            xxx1 = xxx1 + pxk * pyjxk;        end                xxx2 = pyjxi * pxi + eps;% log里的分母                mylog = log2((xxx1+eps)/xxx2); % log                I_part2 = I_part2 + pxi*pyjxi*mylog;    endendI_part2 = I_part2 * (-1);I = I_part1 + I_part2;system_property.IXY = I;end

麻烦就麻烦在：

1. 好像要简单用一下Bayes公式（好像我是在Z99火车上写的。

之前的博文也提到了，我是去广州看张学友演唱会。。。

周围的人都把我看着。看得出来搞这东西好像有点另类）

2. 三重求和有点麻烦（好像是3重……）

纵轴应该就是信道容量了吧。