cpd_pam.m

pam 容量分析 matlab程序 原创

function C=CPD_PAM(snr_db,alpha0,alpha1,dbg)

if(nargin<1)
    snr_db = [-5:10];
    dbg=1;    
    alpha0 = [-3 -1];
    alpha1 = [1 3];
end
if(length(alpha0)~=length(alpha1))
    error('worng input alpha0 and alpha1: the length not match!');
end
alpha = [alpha0 alpha1];
M = length(alpha);
P = sum(abs(alpha.^2))./(M);
alpha = alpha./sqrt(P);
P = sum(abs(alpha.^2))./(M);
alpha0 = alpha(1:M/2);
alpha1 = alpha(M/2+1:end);

snr = 10.^(snr_db./10);
len=length(snr_db);
C=zeros(1,len);
for i=1:len
    sigma2 = P./snr(i);
    if(snr_db(i)>14)
        C(i)=1-quadv(@(x)intfun(x,alpha0,alpha1,sigma2),-2,2,1e-7)./(M*sqrt(2*pi*sigma2));
    elseif(snr_db(i)>0)
        C(i)=1-quadv(@(x)intfun(x,alpha0,alpha1,sigma2),-4,4,1e-7)./(M*sqrt(2*pi*sigma2));
    else
        C(i)=1-quadv(@(x)intfun(x,alpha0,alpha1,sigma2),-5,5,1e-7)./(M*sqrt(2*pi*sigma2));
    end
end
if dbg
    plot(snr_db,C);
end




function z = intfun(x,alpha0,alpha1,sigma2)
t0 = sum(exp(-(x-alpha0).^2/(2*sigma2)));
t1 = sum(exp(-(x-alpha1).^2/(2*sigma2)));
z = t0.*log2(1+t1./t0)+t1.*log2(1+t0./t1);