ica_cdma.m

程序对盲信号分离空时多用户检测器的仿真 但在程序中A用了随机矩阵

clear all
clc

%----------------初始量的声明------------------------------
K=4;                                        %输入信号的个数
len_seq=50;                                 %各个信号序列的长度
m=5;                                        %PN序列的阶数             
len_chip=2^m-1;                             %PN序列的长度
sgma=1;                                    	%噪声平均功率
snr=20;                                     %输入信噪比
Eb=2*sgma^2*10^(snr/10);                    %发送信号的功率
E_chip=Eb/len_chip;                         %单位片长的功率
initial_code=1;                             %差分编码的初始值
R=30;                                       % Analog Signal Simulation rate

%----------------随机生成K个源信号---------------------------
orig_seq=sign(rand(K,len_seq)*2-0.5);

%--------------对orig_seq进行差分编码-----------------------
x=(orig_seq>0);                             %双极性码转化为单极性码
y=zeros(K,len_seq);
%x
for k=1:K                                   %循环为差分编码
    temp_code=initial_code;
    for i=1:length(x(1,:))
            y(k,i)=xor(x(k,i),temp_code);
            temp_code=y(k,i);
    end
end
d_code_seq=(y==1)-(y==0)  ;                 %单极性转换为双极性，d_code_seq为差分编码后的序列

%---------------扩频编码----------------------------------------------------

%-----------产生K个GOLD码序列-----------------

gold_seq=gold_generate(m);
pn=(gold_seq==1)-(gold_seq==0);

%----------对各个序列进行扩频------------------
for i=1:K
    spreaded_seq(i,:)=ds_sp(d_code_seq(i,:),pn(i,:),len_seq,len_chip);
end

%---------进行波形整形，通过的是Rectangular Filter-----------------------------------------
for i=1:K
chips_out=spreaded_seq(i,:)';
chips = [chips_out, zeros(len_chip*len_seq, R-1)];	% Pulse generation
chips = reshape(chips.' , len_chip*len_seq*R, 1);

L = R;
L_2 = floor(L/2);
MF = ones(L, 1);
MF = MF/sqrt(sum(MF.^2));                            %MF为 mathed filter taps
TxOut = sqrt(R)*conv(MF, chips)/sqrt(2);
TxOut = TxOut(L_2+1: end - L_2);
tx(i,:)=TxOut';
end
spreaded_seq=tx;                              %这里的spreaded_seq为整形后的信号矩阵  

%*******************信号的混合加噪声**********************************

A=rand(K);                                     %A为混合矩阵
X=A*sqrt(E_chip)*spreaded_seq+sgma*randn(K,length(spreaded_seq));   %X为混迭信号矩阵

%**********************************************
for i=1:K
    X(i,:)=X(i,:)-1/(len_chip*len_seq*R)*sum(X(i,:));
end

%----------对混迭信号进行预白化---------------------------------------------

[y x]=eig(cov(X',1));
%---特征值按大小排，相应的特征矢量也进行排序----
for i=1:K
   for j=(i+1):K
       temp=x(i,i);
       temp_2=y(:,i);
       if temp<x(j,j)
           x(i,i)=x(j,j);
           x(j,j)=temp;
           y(:,i)=y(:,j);
           y(:,j)=temp_2;
       end
   end
end
d=x;
v=y;
Y=inv(d.^(1/2))*v'*X;           %Y为白化后的输出混迭信号矩阵

%----------用快速ICA算法实现信号分离-----------------------------------------

epsilon=0.001;
W=rand(K);
for p=1:K
    W(:,p)=W(:,p)/norm(W(:,p));
    exit=0;
    count=0;
    iter=1;
    while ((exit==0)&(iter<1000))
          count=count+1;
          temp=W(:,p); %记录上次迭代的值
         W(:,p)=1/(len_chip*len_seq*R)*Y*((temp'*Y).^3)'-3*temp;
          ssum=zeros(K,1);
          for counter=1:p-1
              ssum=ssum+(W(:,p)'*W(:,counter))*W(:,counter);   
          end
          W(:,p)=W(:,p)-ssum;  %正交化
          W(:,p)=W(:,p)/norm(W(:,p));
          if(abs((dot(W(:,p),temp)))<1+epsilon)&(abs((dot(W(:,p),temp)))>1-epsilon)    %判断是否收敛
                        exit=1;
          end
          iter=iter+1;
    end
  end
out=W'*Y;                               %out为输出矩阵

%**************进行匹配滤波后抽样******************
for i=1:K    
    RxIn=out(i,:)';
    L = length(MF);                      %MF为 mathed filter taps
    L_2 = floor(L/2);
    rr = conv(flipud(conj(MF)), RxIn);
    rr = rr(L_2+1: end - L_2);

    Rx = sign(real(rr(1:R:end))) ;
    rx(i,:)=Rx';
end

huang=sign(rx);                         %huang为匹配滤波抽样后的矩阵
%*******************通过判断与PN的相关性确定各个信号的顺序********************
y=zeros(1,4);
for i=1:K
    [x_1(i,1:K) y(i)]=xiang_guan(huang(i,:),pn,len_seq,len_chip,K);
end
 x_1                                %得到的相关性矩阵 x_1=[huang(1,:)*pn1 huang(1,:)*pn2 huang(1,:)*pn3 huang(1,:)*pn4
                                    %                     huang(2,:)*pn1 huang(2,:)*pn2 huang(2,:)*pn3 huang(2,:)*pn4
                                    %                     huang(3,:)*pn1 huang(3,:)*pn2 huang(3,:)*pn3 huang(3,:)*pn4
                                    %                     huang(4,:)*pn1 huang(4,:)*pn2 huang(4,:)*pn3 huang(4,:)*pn4]
                                    %上面的huang(i,:)*pnj,代表他们之间的相关值
zhuan=zeros(K);                     %zhuan为信号顺序转换矩阵
for i=1:K
    zhuan(y(i),i)=1;
end

seq_separated=zhuan*huang;          %seq_separated为已经变为与原来信号矩阵顺序相同

%**************************************************************************

%--------分离得到的信号进行解扩---------------------------------------------
rec_seq_temp=despread_spectrum(seq_separated,pn,len_seq,len_chip);

%--------为测误码率而进行修正相反的幅度----------------------------------------
rec_seq=correct(d_code_seq,rec_seq_temp);

%%-----------差分译码过程------------------------------------------------------
y=rec_seq>0;

for k=1:K
    temp_dcode=initial_code;
for i=1:length(y(1,:))
    z(k,i)=xor(temp_dcode,y(k,i));
    temp_dcode=y(k,i);
end
end

p=2*z-1;
rec_seq=p;

%-----------统计误差-------------------------------------------------------
num_error=sum((rec_seq~=orig_seq)');
p_error=sum(num_error)/(len_seq*K)                  %误码率