M-PAM信道容量仿真：

Matlab程序代码：

snr_indB = -10:1:30;   %定义信噪比单位是db

snr = 10 .^ (0.1 * snr_indB);%信噪比，单位无

 

EbN0 = [];

C_mpam = [];

logpy = [];

alphabet = [];

trials = 50000;

 

% b = 3;

% M = 2 ^ b;

M=2;

while M<=16

for k = 1 : 1 : M;

alphabet(k) = 2*k - 1 - M;

end

average = sum(alphabet .^ 2) / M; % average energy Es

 

for i = 1 : 1 : length(snr_indB);

sigma2 = average / snr(i);%sigma2指N0

 

for j = 1 : 1 : trials;

x(j) = alphabet(floor(rand * M) + 1);

end

 

y = x + sqrt(sigma2) * randn(1, trials);

for j = 1 : 1 : trials;

pyx = 0; % compute p(y) from p(y|x)

for k = 1 : 1 : M;

pyx = pyx + exp(-(y(j) - alphabet(k)).^2 / (2*sigma2));

end

logpy(j) = log2(1/M / sqrt(2*pi*sigma2) * pyx);%??

end

% I(X; Y) = H(Y) – H(Y|X) = H(Y) – H(N)

 

C_mpam(i) = -mean(logpy) - 0.5*log2(2*pi * exp(1) * sigma2);

end

EbN0 = 10 * log10(snr ./ C_mpam / 2);

plot(snr_indB, C_mpam, 'k-');

%plot(EbN0, C_mpam, 'k-');

hold on;

M=M*2;

end

 

 

C = []; % Shannon capacity

for i = 1 : 1 : length(snr_indB);

C(i) = 0.5 * log2(1 + snr(i));

EbN0(i) = 10 * log10(snr(i) / C(i) / 2);

end

plot(snr_indB, C, 'b-');

%plot(EbN0, C, ‘b-‘);

legend('M-PAM bound', 'Shannon bound',-1);

xlabel('SNR (dB)');

%xlabel('Eb/N0 (dB)');

ylabel('Capacity (bits/symbol)');

grid on;

text(25,1.2,'2-PAM');

text(25,2.2,'4-PAM');

text(25,3.2,'8-PAM');

text(25,3.7,'16-PAM');

 

M-QAM信道容量仿真：

Matlab程序代码：

 

snr_indB = -10:1:35;

snr = 10 .^ (0.1 * snr_indB);  % = Es / N0

EbN0 = [];

C_qam = [];

logpy = [];

alphabet = [];

trials = 10000;

%b = 6;

for b=2:1:6;

M = 2 ^ b;

 

if b < 4

    for k = 1 : 1 : M;

        alphabet(k) = complex(cos(2*pi*k/M), sin(2*pi*k/M));

    end

    average = 1; %sum(alphbet .^ 2) / M;

    %d = sqrt(snr ./ average);

elseif b == 4      % for 16QAM

    alphabet(1) = complex(3, 3);

    alphabet(2) = complex(1, 3);

    alphabet(3) = complex(-1, 3);

    alphabet(4) = complex(-3, 3);

    alphabet(5) = complex(3, 1);

    alphabet(6) = complex(1, 1);

    alphabet(7) = complex(-1, 1);

    alphabet(8) = complex(-3, 1);

    alphabet(9) = complex(3, -1);

    alphabet(10) = complex(1, -1);

    alphabet(11) = complex(-1, -1);

    alphabet(12) = complex(-3, -1);

    alphabet(13) = complex(3, -3);

    alphabet(14) = complex(1, -3);

    alphabet(15) = complex(-1, -3);

    alphabet(16) = complex(-3, -3);

elseif b == 5       % for 32QAM

    alphabet(1) = complex(3, 5);

    alphabet(2) = complex(1, 5);

    alphabet(3) = complex(-1, 5);

    alphabet(4) = complex(-3, 5);

    alphabet(5) = complex(5, 3);

    alphabet(6) = complex(3, 3);

    alphabet(7) = complex(1, 3);

    alphabet(8) = complex(-1, 3);

    alphabet(9) = complex(-3, 3);

    alphabet(10) = complex(-5, 3);

    alphabet(11) = complex(5, 1);

    alphabet(12) = complex(3, 1);

    alphabet(13) = complex(1, 1);

    alphabet(14) = complex(-1, 1);

    alphabet(15) = complex(-3, 1);

    alphabet(16) = complex(-5, 1);

    alphabet(17) = complex(5, -1);

    alphabet(18) = complex(3, -1);

    alphabet(19) = complex(1, -1);

    alphabet(20) = complex(-1, -1);

    alphabet(21) = complex(-3, -1);

    alphabet(22) = complex(-5, -1);

    alphabet(23) = complex(5, -3);

    alphabet(24) = complex(3, -3);

    alphabet(25) = complex(1, -3);

    alphabet(26) = complex(-1, -3);

    alphabet(27) = complex(-3, -3);

    alphabet(28) = complex(-5, -3);

    alphabet(29) = complex(3, -5);

    alphabet(30) = complex(1, -5);

    alphabet(31) = complex(-1, -5);

    alphabet(32) = complex(-3, -5);

else    % for 64QAM,b = 6

    alphabet(1) = complex(3, 3);

    alphabet(2) = complex(1, 3);

    alphabet(3) = complex(-1, 3);

    alphabet(4) = complex(-3, 3);

    alphabet(5) = complex(3, 1);

    alphabet(6) = complex(1, 1);

    alphabet(7) = complex(-1, 1);

    alphabet(8) = complex(-3, 1);

    alphabet(9) = complex(3, -1);

    alphabet(10) = complex(1, -1);

    alphabet(11) = complex(-1, -1);

    alphabet(12) = complex(-3, -1);

    alphabet(13) = complex(3, -3);

    alphabet(14) = complex(1, -3);

    alphabet(15) = complex(-1, -3);

    alphabet(16) = complex(-3, -3);

 

    for k = 1 : 1 : 16;

        alphabet(k) = alphabet(k) + complex(4, 4);

    end

    for k = 1 : 1 : 16;

        alphabet(k + 16) = alphabet(k) + complex(-8, 0);

    end

    for k = 1 : 1 : 16;

        alphabet(k + 32) = alphabet(k) + complex(-8, -8);

    end

    for k = 1 : 1 : 16;

        alphabet(k + 48) = alphabet(k) + complex(0, -8);

    end

end

average = sum(abs(alphabet).^2) / M;   % average Es

 

for i = 1 : 1 : length(snr_indB);

    sigma2 = average / (2 * snr(i));

   

    for j = 1 : 1 : trials;

        x(j) = alphabet(floor(rand * M) + 1);

    end

    y = x + sqrt(sigma2)*complex(randn(1, trials), randn(1, trials));

    for j = 1 : 1 : trials;

        pyx = 0;

        for k = 1 : 1 : M;

            pyx = pyx + exp(-(abs(y(j) - alphabet(k)).^2) / (2*sigma2));

        end

        logpy(j) = log2(1/M / (2*pi*sigma2) * pyx);

    end

    C_qam(i) = -mean(logpy) - log2(2*pi * exp(1) * sigma2);

end

 

EbN0 = 10 * log10(snr ./ C_qam);

%plot(snr_indB, C_qam, 'k-');

plot(EbN0, C_qam, 'k-');

hold on;

 

C = [];

for i = 1 : 1 : length(snr_indB);

    C(i) = log2(1 + snr(i));

    EbN0(i) = 10 * log10(snr(i) / C(i));

end

%plot(snr_indB, C, 'b-');

plot(EbN0, C, 'b-');

 

legend('M-QAM bound', 'Shannon bound');

%xlabel('SNR (dB)');

xlabel('Eb/N0 (dB)');

ylabel('Capacity C (bits/symbol)');

grid on;

text(25,2.5,'QPSK' )

text(25,3.5,'8PSK')

text(25,4.5,'16QAM')

text(25,5.5,'32QAM')

text(25,6.5,'64QAM')

end