File talk:Channel Capacity for complex constellations.jpg

source

If you post the Matlab source as well, I can create a better SVG version Alessio Damato 17:32, 15 May 2007 (UTC)

Matlab Source code

I haven't quite figured out all of the wikipedia markup. So to view this code, just edit the entry and copy and paste.

% 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

figure(1);clf; hold on;
subplot(211);
colors = 'rgbmc';
for i=1:5
    plot(log10(SNR)*10, C(i,:),colors(i),'LineWidth',4); hold on;
end
legend('BPSK', 'QPSK', '8 PSK', '16 PSK', '16 QAM','Location','NorthWest');
xlabel('SNR db');
ylabel('Capacity (bits per channel use)');


% This second plot is for the capacity plots vs EB/N0 but I'm pretty
% sure there is a bug in it, so don't use it.
subplot(212); hold on;
for i=1:5
    plot(log10(EBNO)*10, C2(i,:),colors(i),'LineWidth',4); hold on;
end
legend('BPSK', 'QPSK', '8 PSK', '16 PSK', '16 QAM','Location','NorthWest');
xlabel('E_b / N_0 db');
ylabel('Capacity (bits per channel use)');
title('Channel Capacity for complex constellations')
% figure(2); clf; 
% the_xp = HW1_pskX(snr, 16);
% plot(real(the_xp), imag(the_xp));


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

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

Speedplane 04:54, 17 May 2007 (UTC)

if you want to post source code, just put in within <pre>...</pre> (as I did for your code). Thanks for the source, I'll use to soon. Alessio Damato 08:43, 17 May 2007 (UTC)