function dsmud_E44
% E4 盲多用戶檢測算法在AWGN信道下的性能
% RLS-LMF算法


clear; 
% clear global;
global  x offset ny nsy pskout px pnseries;
% 為提高算法修改效率，不重復(fù)計算調(diào)制解調(diào)過程
% 若需改變數(shù)據(jù)，請操作 clear all;

seedx=5;
seedoffset = 3;
seed = 8;

ds_parm = struct('pn',{[]},'pnoffset',{[]},'fd',{[]},...
                  'fc',{[]},'fs',{[]},'N',{[]},'mod',{[]},'rdcdma',{[]});

global pnseries;
pn=[1,0,0,0,0,1,0,0,0]; % PN generation polynominal
pnseries = [pn,pncode(pn,0,2^length(pn)-1,0)];
rand('state',seedoffset);
if length(pskout) == 0
    offset=round(32*rand);
end;

% modulation
Fs = 8;
Fd =1;
Fc = 2;
N = 32;

global h;
% ----------------UWB Multipath Fading Channel---------------------
cm =1;
arsv =2/Fs; %ns
Nch =17; % channel realization of num_channels
% cms = sprintf('cm%d_imr.mat',cm);load(cms);
% % now reduce continuous-time result to a discrete-time result
% ts = 0.01; % ns, 因為uwb_sv_cnvrt_ct所得最小解析度為100Ghz,即0.01ns,取最小解析度
% [hN,Nts] = uwb_sv_cnvrt_ct( h_ct, t_ct, np, num_channels, ts );
% ts = ts /Nts; % nsec
% [h,rsv]=shrinkMP(hN(:,Nch)',ts,arsv);

% ----------------Narrowband Multipath Fading Channel---------------------
% h = [0.8,zeros(1,7),0.8,zeros(1,15),0.6,zeros(1,15),0.4,zeros(1,7),0.3]; 

% ----------------NO Multipath Fading Channel (AWGN)---------------------
h = [1];

h = shrinkMP(h,arsv,2);
h = inflateMP(h,2,arsv);

snr = 10;% AWGN channel
tmp = find(abs(h)>0.1*max(abs(h)));
mph = tmp(1);
h = h(mph:end);

if h(1) < 0
    h = -h;     % the stongest path is innegative
end;
mph =1;

usrnum =6;
offset_step = 34;  % users' pn offset steps, must > N

% Original Data
bitnum = 4000;

ds_parm.pn = pn;
ds_parm.pnoffset = offset;
ds_parm.fd = Fd;
ds_parm.fs = Fs;
ds_parm.fc = Fc;
ds_parm.N = N;;
ds_parm.mod = 'psk';
ds_parm.rdcdma = 'D';  % 在采用短碼'D'時，由于用戶間相關(guān)性不變，造成
                        % 某些用戶干擾很大，有些用戶無多用戶干擾的情況
                        % 必須注意排除這種因素對結(jié)果的影響

s1 = sprintf('\n---CDMA Multiuser Detection------------------------------');
s2 = sprintf('--- CDMA AWGN Channel SNR=%ddB',snr);
s4 = sprintf('---User number=%d and Asynchronous ',usrnum);
s5 = sprintf('---pnoffset=%d, N=%d',offset,N);
disp(s1);disp(s2);disp(s4);disp(s5);

% 如果全局變量中已經(jīng)有解調(diào)結(jié)果，則不再進行調(diào)制解調(diào)過程，直接進入自適應(yīng)算法！
if length(pskout) == 0   

rand('state',seedx); 
x = round(rand(usrnum,bitnum));    

rand('state',seed);
timeoffset = zeros(1,usrnum-1);%floor(N*Fs*rand(1,usrnum-1));
% modulation-------------------------------------------
sy = dspsk_MOD(x(1,:),ds_parm);  
y = sy;
for usr=1:usrnum -1
    ds_parm.pnoffset = ds_parm.pnoffset + offset_step;
    tmp = 3*dspsk_MOD(x(usr+1,:),ds_parm);  
    y = y +[tmp(timeoffset(usr)+1:end),zeros(1,timeoffset(usr))];% 異步多用戶信號
end;
ds_parm.pnoffset = offset; % 恢復(fù)參數(shù)

% channel
my = conv(y,h);
% yf = fft([y zeros(1,length(h)-1)]);     % only for speed up convolution
% hf = fft([h zeros(1,length(y)-1)]);
% my = real(ifft(yf.*hf));

msy = conv(sy,h);
% syf = fft([sy zeros(1,length(h)-1)]);     % only for speed up convolution
% shf = fft([h zeros(1,length(sy)-1)]);
% msy = real(ifft(syf.*shf));

snr = snr-10*log10(0.5)-10*log10(Fs/Fd); % db
ny=awgn(my,snr,'measured');
nsy = awgn(msy,snr,'measured');

% % output from matchfilter, sampled at 'chip' rate, synchronized with desired
% % user's signal 
% common direct matchfilter detection
[px,pskout] = dspsk_DEMOD(ny(mph:end),ds_parm);  
end;  % 如果全局變量中已經(jīng)有解調(diào)結(jié)果，則直接進入自適應(yīng)算法！

perror = length(find(x(1,:)-px(1:bitnum))~=0);
s = sprintf('\ndirect MUD(Matchfilters), error = %d /%d',perror,bitnum);
disp(s);

[cx,BWc,esthc] = Demod_my( pskout, ds_parm,1 );
bits = min(length(x(1,:)),length(cx));
cerror = length(find(x(1,1:bits)-cx(1:bits))~=0);
s = sprintf('\nAdaptive Constrant CMA multiuser detection, error = %d /%d',cerror,bits);
disp(s);

[dx,BWd,esthd] = Demod_my( pskout, ds_parm,2 );
bits = min(length(x(1,:)),length(dx));
derror = length(find(x(1,1:bits)-dx(1:bits))~=0);
s = sprintf('\nAdaptive DCMA multiuser detection, error = %d /%d',derror,bits);
disp(s);

[wx,BWw,esthw] = Demod_my( pskout, ds_parm,3 );
bits = min(length(x(1,:)),length(wx));
werror = length(find(x(1,1:bits)-wx(1:bits))~=0);
s = sprintf('\nAdaptive E4 multiuser detection, error = %d /%d',werror,bits);
disp(s);

[ox,BWm,esthm] = Demod_my( pskout, ds_parm,4 );
bits = min(length(x(1,:)),length(ox));
oerror = length(find(x(1,1:bits)-ox(1:bits))~=0);
s = sprintf('\nAdaptive CMOE multiuser detection, error = %d /%d',oerror,bits);
disp(s);

[rx,BWr,esthr] = Demod_my( pskout, ds_parm,5 );
bits = min(length(x(1,:)),length(rx));
rerror = length(find(x(1,1:bits)-rx(1:bits))~=0);
s = sprintf('\nAdaptive E4-RLS multiuser detection, error = %d /%d',rerror,bits);
disp(s);


return;


% functions----------------------------------------------------------------
function  [bx,W,esth] = Demod_my( E,ds_parm,method)
pn =ds_parm.pn;
pnoffset = ds_parm.pnoffset;
Fd = ds_parm.fd;
Fc = ds_parm.fc;
Fs = ds_parm.fs;
N = ds_parm.N;
mod = ds_parm.mod;
rdcdma = ds_parm.rdcdma;

pnseries = pncode(pn,pnoffset,N,0);
S1 = pnseries*2-1;
W=[];

[le,syms] = size( E ); 
L = N;
%calculate Maltipath----------------------------
global h;
ts = 0.25; %ns
arsv = 2; %ns
[hout,rsv] = shrinkMP(h',ts,arsv);
thres=0.1*max(abs(hout));                   % 調(diào)節(jié)多徑延遲解析范圍
for ml = 0:length(hout)-1
    if (abs(hout(end-ml))>thres)
        break;
    end;
end;
ML = length(hout)-ml;
% construct 'yy'------------------
yy = E(:,1:syms-3); % 3*Ns*Nf is the maxium length of maltipath profile
mml = ML-1;
while mml > 0
    if mml > L
        yy = [yy;E(:,1:syms-3)];
    else
        yy = [yy;E(1:mml,1:syms-3)];
    end;
    mml = mml - L ;
end;

% construct multipath C------------
C  = [];
for i=1:ML
    cc = [zeros(1,i-1),S1,zeros(1,ML-i)];
    C=[C,cc'];
end;

esth = [];
switch method
    case 1      % CMA algorithm with constrain
        [bx,W,esth]=Cma1(yy,S1,C);
    case 2      % DCMA algorithm
        [bx,W,esth]=DCma(yy,S1,C);
    case 3      % minimum Shalvi-Weistein modulus algorithm(SW)
        [bx,W,esth]=E4(yy,S1,C); 
    case 4      % CMOE
        [bx,W,esth]=CMOE(yy,S1,C);
    case 5      % E4-RLS algorithm
        [bx,W,esth]=E4RLS(yy,S1,C);
    
end;

return;

function  [bx,ww,esth] = Cma1(yy,S1,C)
L = length(S1);
[a,ML] = size(C);
[a,syms] = size( yy ); 

invcc = pinv(C'*C);
P  = eye(L+ML-1) - C * invcc * C';
G = zeros(ML,1);G(1)=1;

% algorithm--------------------
mu = 0.00001;%0.00004;
Em=[];
err=[];
W = zeros(size(C(:,1)));
ww=zeros(length(W),syms);

for k=1:syms
    r   = yy(:,k);	% new input vector
    y = W'*r;
    grad = y*(y^2-1);
    % W  = W - mu*P*grad*r;
    W=P*(W-mu*grad*r)+C*invcc*G;
    ww(:,k)=W;
    
    Z   = W' * r; 
    if( Z > 1e+5 )                      
        error('CMA1 cannot convergence with LMS, LMS stopped!'); return;
    end;
    err = [err,y^2];
    Em = [Em,Z];
end;

% Em = W'*yy;     % 如果考察收斂正確性
bx = round( (sign(-Em)+1)/2 );
esth = [];

return;

function  [bx,W,esth] = DCma(yy,S1,C)
L = length(S1);
[a,ML] = size(C);
[a,syms] = size( yy );

invcc = pinv(C'*C);
P  = eye(L+ML-1) - C * invcc * C';
G = zeros(ML,1);G(1)=1;

% algorithm--------------------
D =10;
mu = 0.00001;%0.00004;
Em=[];
err=[];
W = zeros(size(C(:,1)));

for k=1:syms
    r   = yy(:,k);	% new input vector
    y = W'*r;
    if k <= D
        yd = 1;
    else
        yd = Em(k-D);%W'*yy(:,k-D);%
    end;
    grad = y*(y^2-yd^2);
    W=P*(W-mu*grad*r)+C*invcc*G;
    Z   = W' * r; 
    if( Z > 1e+5 ) 
        error('DCMA cannot convergence with LMS, LMS stopped!'); return;
    end;
    err = [err,y^2];
    Em = [Em,Z];
end;

% Em = W'*yy;     % 如果考察收斂正確性
bx = round( (sign(-Em)+1)/2 );
esth = [];
return;

function  [bx,W,esth] = E4(yy,S1,C)
L = length(S1);
[a,ML] = size(C);
[a,syms] = size( yy );

invcc = pinv(C'*C);
P  = eye(L+ML-1) - C * invcc * C';
G = zeros(ML,1);G(1)=1;

% algorithm--------------------
mu = 0.00001;%0.000003;
Em=[];
zz=[];
W = zeros(size(C(:,1)));

for k=1:syms
    r   = yy(:,k);	% new input vector
    y = W'*r;
    grad = y^3;
%     
%     if abs(grad) > 10
%         grad =sign(grad)/10;
%     end;
    
    W=P*(W-mu*r*grad)+C*invcc*G;

    Z   = W' * r; 
    if( Z > 1e+5 ) 
        error('E4 cannot convergence with LMS, LMS stopped!'); return;
    end;    
    Em = [Em,Z];

end;

Em = W'*yy;     % 如果考察收斂正確性
bx = round( (sign(-Em)+1)/2 );
esth = [];
return;



function [bx,W,esth]=CMOE(yy,S1,C)

L = length(S1);
[a,ML] = size(C);
[a,syms] = size( yy ); 

mu = 0.000011;
invcc = pinv(C'*C);
P  = eye(L+ML-1) - C * invcc * C';
invc = pinv(C);

Em=[];
I=eye(ML,1);    
W = invc'*I;
for k=1:syms
    y   = yy(:,k);	% new input vector
    grad = y*y'*W;
    W = W - mu*P*grad;

    Z = W' * y;
    Em  = [Em,Z];
end;

Em = W'*yy;     % 如果考察收斂正確性
bx = round( (sign(-Em)+1)/2 );
esth = [];

return;


function  [bx,W,esth] = E4RLS(yy,S1,C)
L = length(S1);
[a,ML] = size(C);
[a,syms] = size( yy );

invcc = pinv(C'*C);
P  = eye(L+ML-1) - C * invcc * C';
G = zeros(ML,1);G(1)=1;

% algorithm--------------------
lambda = 1-0.0015;
Em=[];
zz=[];
W = zeros(size(C(:,1)));
invRn = 0.01*eye(L);

for k=1:syms
    x   = yy(:,k);	% new input vector
    r = x*x'*W;
    kn = invRn * r / ( lambda + r' * invRn * r );
    invRn = ( invRn - kn * r' * invRn )/lambda;
    hn = invRn * S1';
    W = hn / ( S1 * hn );
    
    Z   = W' * x; 
    if( Z > 1e+5 ) 
        error('E4-RLS cannot convergence, stopped!'); return;
    end;    
    Em = [Em,Z];

end;

Em = W'*yy;     % 如果考察收斂正確性
bx = round( (sign(-Em)+1)/2 );
esth = [];
return;