%參數：
% M-ary UWB System
M=8
k=log2(M)
SNR_dB=[20];                      %這是 DIRMA receiver 的output SNR
data_size=200;                     %信源數據數
re_imp_width=0.7*1e-9;              %脈沖寬度 以秒為單位
re_imp_width_para=0.2877*1e-9;      %以秒為單位
frame_length=100*1e-9;              
data_ex=0.7*1e-9;         %數據的時間偏移量等于脈沖寬度
time_slot=fix(frame_length/re_imp_width);
if rem(time_slot,2)~=0;
    time_slot=time_slot-1;
end;
Ns=M;                  % 這里讓 Ns=M;
User_num=[1,10];              %最好小于time_slot/2;這樣沖突概率大大減小.
a=length(SNR_dB);
b=length(User_num)
SNR=10.^(SNR_dB/10);
BER=zeros(b,a);
for iu=1:b
    for i=1:a
        err_num=0;
        for ii=1:data_size
            ii
            inf_source=randint(User_num(iu),k);
            %將信源信息轉換成10進制
            base=2.^(0:k-1);
            inf_dec=inf_source*base';
            % M*M 的 Walsh 矩陣
            Walsh=Hadamard(M);
            inf_Walsh=-1*ones(User_num(iu),M);
            for i1=1:User_num(iu)
                inf_Walsh(i1,:)=Walsh(inf_dec(i1)+1,:);   
            end;
            %擴頻
            PN_code=fix(rand(User_num(iu),Ns)*time_slot/2)*2+1;   %要求time_slot是偶數,這是產生的偽隨機碼序列
            spread_data=-1*ones(User_num(iu),Ns*time_slot);       % "-1"表示此處沒有脈沖傳送
            for i1=1:User_num(iu)
                for i2=1:Ns
                    if inf_Walsh(i1,i2)==1;
                        spread_data(i1,PN_code(i1,i2)+1+(i2-1)*time_slot)=inf_Walsh(i1,i2);
                    end;
                    if inf_Walsh(i1,i2)==0;
                        spread_data(i1,PN_code(i1,i2)+(i2-1)*time_slot)=inf_Walsh(i1,i2);
                    end;
                end;
            end;

            %理想接收的波形,抽樣數為8;
            step=0.1;
            x=0:step:0.7;
            rec_w=(1-4.*pi.*((x-0.35)/0.2877).^2).*exp(-2.*pi.*((x-0.35)/0.2877).^2);

            %接收的信號(含波形抽樣);
            mod_data=zeros(User_num(iu),Ns*time_slot*length(x));
            rec_sig=zeros(1,Ns*time_slot*length(x));
            for i1=1:User_num(iu)
                for i2=1:Ns*time_slot
                    if spread_data(i1,i2)~=-1;
                        mod_data(i1,(i2-1)*length(x)+1:i2*length(x))=rec_w;
                    end;
                end;
            end;
        
            local_sig_n_PN=zeros(1,2*length(x));
            mono_fir=zeros(1,2*length(x));
            local_sig_n_PN(1,1:length(x))=rec_w;
            local_sig_n_PN(1,(length(x)+1):end)=-rec_w;
            mono_fir(1,(length(x)+1):end)=rec_w;
            sig_pow=(Ns*local_sig_n_PN*mono_fir'*step)^2;               % 算DIRMA receiver 的output signal power
            noise_pow_add_mul=sig_pow/SNR(i);
            %%%%%%%%%%%%%% 用符號計算計算多址噪聲%%%%%%%%%%%%%
            syms t r
            if ii==1
                f=(1-4*pi*((t-r-0.35)/0.2877)^2)*exp(-2*pi*((t-r-0.35)/0.2877)^2);
                g=(1-4*pi*((t-0.35)/0.2877)^2)*exp(-2*pi*((t-0.35)/0.2877)^2);
                h=(1-4*pi*((t-0.35-0.7)/0.2877)^2)*exp(-2*pi*((t-0.35-0.7)/0.2877)^2);
                z=(g-h)*f;
                d=int(int(z,t,-inf,inf)^2,r,-inf,inf);
                b=double(vpa(d)/100);
            end;
            
            noise_pow=noise_pow_add_mul-(User_num(iu)-1)*b*Ns;
            rec_sig=sum(mod_data,1)+sqrt(noise_pow)*randn(1,Ns*time_slot*length(x));
    
            %%%%%%%%%%%%%%%%  檢測信號 %%%%%%%%%%%%%%%%%%
            %本地產生的模版信號
            %注意:這里只計算第一個用戶的誤碼率
            local_data=-1*ones(1,Ns*time_slot); % "-1"表示此處沒有脈沖傳送
            for i1=1:Ns
                local_data(1,PN_code(1,i1)+(i1-1)*time_slot)=0;
                local_data(1,PN_code(1,i1)+1+(i1-1)*time_slot)=1;
            end;
            local_sig=zeros(1,Ns*time_slot*length(x));
            for i1=1:Ns*time_slot
                if local_data(1,i1)==0;
                    local_sig(1,(i1-1)*length(x)+1:i1*length(x))=rec_w;
                end;
                if local_data(1,i1)==1;
                    local_sig(1,(i1-1)*length(x)+1:i1*length(x))=-rec_w;
                end;
            end;
    
            %產生判決量
            corr=zeros(1,Ns);
            temp=time_slot*length(x);
            for i1=1:Ns
                corr(i1)=rec_sig(1,((i1-1)*temp+1):i1*temp)*(local_sig(1,((i1-1)*temp+1):i1*temp))'*step;
            end;
            dec=zeros(1,M);
            for i1=1:M
                dec(i1)=corr*Walsh(i1,:)';
            end;
            [temp1,dec_data_decimal]=max(dec);
            dec_data_decimal=dec_data_decimal-1;
            %將十進制轉換成二進制
            dec_data=zeros(1,k);
            for i1=1:k
                dec_data(i1)=fix(dec_data_decimal/(2^(k-i1)));  %將state_dec轉換成二進制
                dec_data_decimal=dec_data_decimal-dec_data(i1)*(2^(k-i1));
            end
            dec_data=dec_data(k:-1:1);
            %計算誤碼個數
            for i1=1:k
                if dec_data(i1)~=inf_source(1,i1);
                    err_num=err_num+1;
                end;
            end;
        end;
        BER(iu,i)=err_num/(k*data_size);
    end;
end;
semilogy(SNR_dB,BER(1,:),'k-p');
hold on;    
semilogy(SNR_dB,BER(2,:),'b-*');
legend('10 User','30 User');