Image:Channel capacity for complex constellations.svg

Plot of the performance of several modulation techniques against the SNR. The calculations were made using an Octave source (but works on Matlab as well), the plot was made using Gnuplot. Here are two helper functions:

function p = HW1_pYget(y, X)
    M = length(X);
    p = 1/M * 1/pi * sum(exp(-abs(y - X).^2));

and

function X = HW1_pskX(SNR, M)
    m = 0:(M-1);
    X = sqrt(SNR)*exp(j*2*pi*m/M);

and here is the code that makes the calculations. It saves the data in a file called channel.dat:

% Channel Capacity for complex constellations
clear all;
% N is the number of samples we're going to be simulating
% To make the plots smoother, jack up this number
N = 3000;

% The x-axis
SNR = logspace(-.3, 2);
EBNO = logspace(-.16, 2);

%Initialize the noise
Z = ((randn(1,N)) + j*(randn(1,N)))/sqrt(2);
qam = [-1 -1/3 1/3 1];
qam16 = [];
for i=1:4
    qam16 = [qam16 (qam(i) + j*qam)];
end
bits = [1 2 3 4 4];

for i = 1:5
    if(i==1)
        the_xs = HW1_pskX(1, 2);
        r = 1;
    elseif(i==2)
        the_xs = HW1_pskX(1, 4);
        r = 2;
    elseif(i==3)
        the_xs = HW1_pskX(1, 8);
        r = 3;
    elseif(i==4)
        the_xs = HW1_pskX(1, 16);
        r = 4;
    elseif(i==5)
        the_xs = qam16;
        r = 4;
end

    for s = 1:length(SNR);
        snr = SNR(s);
        the_x = the_xs*sqrt(snr);
        M = length(the_x);
        X = the_x(randint(1,N,M)+1);
        
        Y = X + Z;
        lP = 0;
        for y = Y
            p = HW1_pYget(y, the_x);
            if(p > 1 || p < 0)
                fprintf('Bad p=%d, y=%d, i=%d\n', p,y,i);
            end
            lP = lP + log2(p);
        end
        HY = -lP/length(Y);
        C(i, s) = HY - log2(pi*exp(1));
    end
    
    for s = 1:length(EBNO);
        snr = (EBNO(s)*r);
        the_x = the_xs*sqrt(snr);
        M = length(the_x);
        X = the_x(randint(1,N,M)+1);
        
        Y = X + Z;
        lP = 0;
        for y = Y
            p = HW1_pYget(y, the_x);
            if(p > 1 || p < 0)
                fprintf('Bad p=%d, y=%d, i=%d\n', p,y,i);
            end
            lP = lP + log2(p);
        end
        HY = -lP/length(Y);
        C2(i, s) = HY - log2(pi*exp(1));
    end
end

%save data to external file
F = [log10(SNR)*10; C(1:5,:)];
F = F';
save -ascii 'channel.dat' F;

The file channel.dat is used in the following Gnuplot code to get the plot:

# set the output
set terminal svg enhanced fname "Times" fsize 18
set output "channel.svg"

# axis properties
set xrange [-5:20]
set yrange [0:4.5]
set xlabel "SNR [dB]"
set ylabel "Channel capacity"
set grid xtics ytics
set key 1.5,4.1

plot "channel.dat" using 1:2 title "BPSK" with lines linewidth 2,\
     "channel.dat" using 1:3 title "QPSK" with lines linewidth 2,\
     "channel.dat" using 1:4 title "8 PSK" with lines linewidth 2,\
     "channel.dat" using 1:5 title "16 PSK" with lines linewidth 2,\
     "channel.dat" using 1:6 title "16 QAM" with lines linewidth 2

that created channel.svg. Finally I have post-processed it with Inkscape and re-named it.

self-made

17th May 2007

see below