%product received vector 
function [p]=BP(snr_in_dB)
Lc=31;                  %Lc is frame size 碼長
snr=10^(snr_in_dB/10);   %信噪比由dB形式轉化               也有的寫為snr=exp(snr_in_dB*log(10)/10)
sgma=1;                  % noise standard deviation is fixed 定義方差
Eb=2*sgma^2*snr;         % signal level required to achieve the given 
                         % signal to noise ratio 
E_chip=Eb/Lc;            % energy per chip 每個碼片能量
A=[sqrt(Eb),sqrt(Eb),sqrt(Eb)];     %sqrt(20*Eb),sqrt(20*Eb),sqrt(20*Eb),sqrt(20*Eb)  
% A為干擾用戶的幅度，0-30，最多不要超過31，干擾強度為20dB,BPSK中
% 干擾強度=10*log(Ak^2/E1),E1為期望用戶的能量,若想實現干擾強度隨snr_in_dB變化,則干擾用戶為sqrt(Eb*exp(snr_in_dB/10))
N=1000;                  % number of bits transmitted  發送的信息位，如果圖效果不好可以改大                 
K=length(A);             % number of users 用戶數，即A的長度
 
% gold sequence for the duration of the bit is generated next 
connections1=[1 0 1 0 0];
connections2=[1 1 1 0 1];
gold_seq=get_gold(connections1,connections2);  %產生31x31的擴頻碼
user1code=(1/sqrt(Lc))*gold_seq(:,31);   %定義了用戶1的擴頻碼，1/sqrt(Lc)表示歸一化，把用戶的擴頻碼定義在第31列，可以減少干擾用戶的影響

% The generation of the data, noise, interference, decoding process and error 
% counting is performed all together in order to decrease the run time of
% the program. This is accomplished by avoiding very large sized vectors. 


%算法賦初值
 num_of_err=0;  %錯誤碼數先賦0
 S=zeros(31,K+1);  %定義擴頻碼矩陣
 S(:,K+1)=user1code; 
 w=[A,sqrt(Eb)];
 W=diag(w);
 
%從此位開始一直到最后一個end，是循環，作用是發送并檢測，用接收與發送比較計算錯誤碼個數

for i=1:10, 
  
  b=zeros(K+1,1);   %定義發送端發送的信息列
  y=zeros(K+1,1);   %定義經過匹配濾波器組相關處理后的信息列  
 
% generate the next data bit  產生新的數據位
  temp=rand;  
  if (temp<0.5),
      data=-1;
  else
      data=1;
  end
  % the transmitted signal is 
  trans_sig=sqrt(Eb)*data*user1code; %與公式相比sqrt(Eb)相當于Ak data相當于bk user1code為sk   sqrt(Eb)認為是期望用戶，
%   與A中干擾用戶幅度定義不同，為了突出最后圖的效果
  
% AWGN with variance sgma^2 
  noise=sgma*randn(Lc,1);   %產生一列31行的隨機噪聲，加性高斯白噪聲  不是sgma^2是因為sgma本身為1（或許可以這么認為）
 
  % interference user data is dat
  % 產生K個干擾用戶，dat是干擾用戶數據，真正的用戶為K+1個，期望用戶直接定義為sqrt(Eb)
  inter_sig=zeros(31,1);

  for k=1:K,
      temp=rand;  
      if (temp<0.5),
          dat=-1;
      else
          dat=1;
      end
      inter_sig=inter_sig+A(k)*dat*(1/sqrt(Lc))*gold_seq(:,k); %干擾用戶數據的擴頻碼放在1到第K列，期望用戶的在第31列
      S(:,k)=(1/sqrt(Lc))*gold_seq(:,k); 
      b(k)=dat;
  end;
  b(K+1)=data;
  
  R=corrcoef(S);  %擴頻碼的互相關矩陣
 
 % received signal  
  rec_sig=trans_sig+noise+inter_sig;   %接收到的數據是31x1位的
 
 

 % making detect                       開始檢測
  for k=1:K+1,
      y(k)=rec_sig'*S(:,k);               %求出經過匹配濾波器組相關處理后的信息列 
  end
 % 建立Bp神經網絡，并加以訓練
  net=newff(minmax(y'),[3,1],{'tansig','purelin'},'traingd','learngd');
  net.trainParam.goal=0.001;
  net.trainParam.lr=0.05;
  net.trainParam.eporchs=300;
  net.trainParam.show=NaN;
  b=b';
  net=train(net,y',b);
end
% 網絡訓練好以后，開始正式檢測
for i=1:N, 
  
  b=zeros(K+1,1);   %定義發送端發送的信息列
  y=zeros(K+1,1);   %定義經過匹配濾波器組相關處理后的信息列  
 
% generate the next data bit  產生新的數據位
  temp=rand;  
  if (temp<0.5),
      data=-1;
  else
      data=1;
  end
  % the transmitted signal is 
  trans_sig=sqrt(Eb)*data*user1code; %與公式相比sqrt(Eb)相當于Ak data相當于bk user1code為sk   sqrt(Eb)認為是期望用戶，
%   與A中干擾用戶幅度定義不同，為了突出最后圖的效果
  
% AWGN with variance sgma^2 
  noise=sgma*randn(Lc,1);   %產生一列31行的隨機噪聲，加性高斯白噪聲  不是sgma^2是因為sgma本身為1（或許可以這么認為）
 
  % interference user data is dat
  % 產生K個干擾用戶，dat是干擾用戶數據，真正的用戶為K+1個，期望用戶直接定義為sqrt(Eb)
  inter_sig=zeros(31,1);

  for k=1:K,
      temp=rand;  
      if (temp<0.5),
          dat=-1;
      else
          dat=1;
      end
      inter_sig=inter_sig+A(k)*dat*(1/sqrt(Lc))*gold_seq(:,k); %干擾用戶數據的擴頻碼放在1到第K列，期望用戶的在第31列
      S(:,k)=(1/sqrt(Lc))*gold_seq(:,k); 
      b(k)=dat;
  end;
  b(K+1)=data;
  
  R=corrcoef(S);  %擴頻碼的互相關矩陣
 
 % received signal  
  rec_sig=trans_sig+noise+inter_sig;   %接收到的數據是31x1位的
 
 

 % making detect                       開始檢測
  for k=1:K+1,
      y(k)=rec_sig'*S(:,k);               %求出經過匹配濾波器組相關處理后的信息列 
  end
 % 調用神經網絡，得到輸出估計值
  output=sim(net,y');
  output(K+1)=sign(output(K+1));  
  if  (output(K+1)<=0),
           bit=-1;
      else
           bit=1;
      end;
     
      if (bit~= data), 
        num_of_err=num_of_err+1; 
      end;
end;       % 對應for i=1:N,
p=num_of_err/N;