starchannelestimation.m

DS-CDMA通信系统中的空时移动信道估计

function STARChannelEstimation();
%************************************************************************
% STARChannelEstimation()
%
% Comments: The purpose of this function is to perform STAR Channel  
% Estimation and the comparison of different types of receivers.
%
% Version 1.0, date : September 2005.
% Written by Farrukh Rashid, Array Communications Research Group,
% Imperial College London, 2005
%
% Ref.: [1] Manikas, A.; Huang, L.K. 
%           STAR channel estimation in DS-CDMA communication systems
%           Communications, IEE Proceedings-, Vol.151, Iss.4, 21 Aug. 2004
%           Pages: 387- 393
%     
%       [2] Manikas, A.; Sethi, M. 
%           A space-time channel estimator and single-user receiver for 
%           code-reuse DS-CDMA systems
%           Signal Processing,IEEE Transactions on, Vol.51, Iss.1,Jan. 2003
%           Pages: 39- 51
%
% LIST OF FUNCTIONS USED:
%
% function gold = Gold_Seq(pol_1,pol_2);            
% function [seq] = genMSeq(conns, initVal); 
% function [gSeq] = genGSeq(seq1, seq2);
% function sseq = shift(seq, k);
% function [Tx_Data_All_Users, Des_User_Data] = Tx_Data( M, Lw, Mod);
% function dataout = Encoder(Bin_data, Mod)
% function [X, X_des, Comp_Manifold_Matrix, H_1, Beta] = ...
%    Received_Signal(M_1, drcts_1, l_1, Beta_1, gold, array, J, Ki, SNR);
% function [DOAs_est, TOAs_est, MUSIC_Spec] = ...
%              STAR_Channel_Estimation(x_final, c1, Rx_array, Q, K_des);
% function [z1, D, phi_vec] = preprocessor(X, c1, Nc, N);
% function [R_z_STsmooth] = Spatial_Temporal_Smoothing(z1, N, Q);
% function [R_z_Ssmooth ] = spat_smooth_over (R_z, Q, N );
% function [R_z_Tsmooth ] = temp_smooth_matrix(R_mat, Nc, Q);
% function [D1_STsmooth, D2_STsmooth] = smoothing_Diagonal(D1, N, Q);
% function [Z]=music ...
% (Rxx,array,M,AZarea,ELarea,TOA_Area, phi_vec_d,mainlobe,ingain,inphase);
% function B_est = ...
%           EstimateFadingCoeff(DOAs_est,TOAs_est,DOA_orig,TOA_orig,B_1);
% function Dat_out = Decoder(Sym_in, Mod)
% function [w_STAR] = ...
%             STAR_Eq_21_Ref_1_Weight(X,DOAs_est,TOAs_est,array,PN_1,B_1); 
% function [w_RAKE] = RAKE_Eq_22_Ref_2_Weight(H_1, B_1);
% function [w_DecD] = DecD_Eq_23_Ref_2_Weight(Comp_Manifold_Matrix);
% function [Sym_Rcvd, dec_data] = Receiver(Rx_wt, X, Mod);
% function [SNIR, BER] = ...
%       Compute_SNIR_BER(Rx_wt, X, X_des, dec_data, Des_User_Data, Mod);
% function S=spv(array,directions,mainlobe);
% function x=frad(degrees);
% function KI=fki(az,el);
% function B=repc(A,m);
% function [Mmdl,Maic,MDL,AIC] = detect(Rxx,L)
% function best = minmatr(Z, AZarea, ELarea, n)
%************************************************************************
clear all; close all;

% SIMULATION PARAMETERS (ACCORDING TO SECTION 5 OF [1])
%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('SIMULATION PARAMETERS.');

% Set the array configuration
array = [-2 0 0;-1 0 0; 0 0 0 ; 1 0 0 ; 2 0 0 ];   % Uniform Linear Array
N = size(array,1);                                 % No. of Array elements

% Set the number of multipaths for all users
Ki = [10 8 6]'; % Number of multipaths for each user(e.g. 1st user 10, 
                % 2nd user 8, 3rd user 6)
M = length(Ki); % Number of users (i.e. M = 3 & 1st user = Desired User)

% Gold Codes used as PN sequence
All_PN_Codes = Gold_Seq([1 0 0 1 0 1],[1 1 1 1 0 1]); % Codes of length 31 
PN_1 = All_PN_Codes(:,10);      % PN sequence for user 1 
PN_2 = All_PN_Codes(:,11);      % PN sequence for user 2
PN_3 = All_PN_Codes(:,12);      % PN sequence for user 3
PNseqs = [PN_1 PN_2 PN_3];      % Matrix of PN sequences of all users
Nc=length(PN_1);                % Processing Gain or Length of each PN Code
J =[zeros(1,(2*Nc-1)) 0;eye(2*Nc-1) zeros((2*Nc-1),1)]; % Eq 8 in [1]

% Multipath DOAs, TOAs & Fading Coefficients (Table 1: Channel Parameters)
DOA_1=[50 0;119 0;100 0;20 0;65 0;20 0;80 0;90 0;70 0;40 0];
DOA_2=[100 0;60 0;27 0;85 0;45 0;30 0;118 0;105 0];
DOA_3=[50 0;10 0;129 0;15 0;25 0;70 0];

TOA_1 = [2 4 5 10 10 15 19 23 24 28]';
TOA_2 = [5 6 12 15 18 24 27 29]';
TOA_3 = [3 4 7 8 11 12]';

FadCoeff_1 = [-0.081 -0.122 -0.009 -0.043 -0.234 0.315 0.107 0.312 ...
                -0.121 -0.273]' + j*[0.197 0.01 -0.184 -0.022 -0.489 ...
                                        0.264 -0.239 -0.168 0.327 -0.222]';
FadCoeff_2 = [-0.969 -1.187 3.799 0.189 -2.536 -2.638 4.072 4.545]' + ...
                j*[3.836 -2.945 -0.501 -0.467 0.721 -1.343 -2.293 -1.048]';
FadCoeff_3 = [-1.362 -1.067 -5.234 -0.845 0.428 1.136]' + ...
                                  j*[5.21 1.267 2.883 3.396 3.396 -3.580]';
DOAs = [DOA_1; DOA_2; DOA_3];
TOAs = [TOA_1; TOA_2; TOA_3];
FadCoeffs = [FadCoeff_1; FadCoeff_2; FadCoeff_3];

% Set the Modulation, Observation Interval and Input Signal to Noise Ratio
Mod = 'BPSK'; % Modulation can be 'BPSK' or 'QPSK' 
Lw = 200;     % Number of Symbols/Observation Interval for 1 realization 
SNR = 40;     % Input Signal to Noise Ratio in dBs
%%%%%%%%%%%%%%%%%%%%%%%%%%%


% STEP 1: GENERATE THE TRANSMITTED DATA OF ALL USERS 
disp('TRANSMITTED DATA OF ALL USERS.');
[Tx_Data_All_Users, Des_User_Data] = Tx_Data( M, Lw, Mod); 


% STEP 2: MODEL THE RECEIVED SIGNAL (SECTION 2 OF [1], Eqns 1 - 10)
disp('RECEIVED SIGNAL MODELLING.');
[X, X_des, Comp_Manifold_Matrix, H_1, B_1] = Received_Signal ...
    (Tx_Data_All_Users, DOAs, TOAs, FadCoeffs, PNseqs, array, J, Ki, SNR);


% STEP 3: STAR CHANNEL ESTIMATION (SECTION 3 OF [1], Eqns 11 - 19) 
disp('STAR CHANNEL ESTIMATION.');
Q = 10;                     % No. of submatrices for temporal smoothing
c1=[PN_1;zeros(Nc,1)];        % Ref Temp Manif Vector for the desired user
[DOAs_est,TOAs_est,MUSIC_Spec]=STAR_Channel_Estimation(X,c1,array,Q,Ki(1));
B_est = EstimateFadingCoeff(DOAs_est,TOAs_est,DOA_1(:,1),TOA_1,FadCoeff_1);


% STEP 4: STAR SUBSPACE RECEPTION (SECTION 4 OF [1], Eqns 20 - 23)
disp('STAR RECEIVER: REPRESENTATIVE EXAMPLE.');
[w_STAR]=STAR_Eq_21_Ref_1_Weight(X,DOAs_est,TOAs_est,array,PN_1,B_est);
[STAR_plot, STAR_dec_data] = Receiver(w_STAR, X, Mod);  % Eq 20,23 in [1]
[SNIR_STAR, BER_STAR] = ...
  Compute_SNIR_BER(w_STAR, X, X_des, STAR_dec_data, Des_User_Data, Mod);


% TWO OTHER RECEIVERS ARE GIVEN FOR COMPARISON (Refer to [2])
disp('ST-RAKE Reception.');
% ST-RAKE Reception according to Eq 22 in [2]
[w_RAKE] = RAKE_Eq_22_Ref_2_Weight(H_1, B_1);
[RAKE_plot, RAKE_dec_data] = Receiver(w_RAKE, X, Mod);
[SNIR_RAKE, BER_RAKE] = ...
  Compute_SNIR_BER(w_RAKE, X, X_des, RAKE_dec_data, Des_User_Data, Mod);

disp('ST-Decorrelating Detector Reception.');
% ST-Decorrelating Detector Recpetion according to Eq 23 in [2]
[w_DecD] = DecD_Eq_23_Ref_2_Weight(Comp_Manifold_Matrix);
[DecD_plot, DecD_dec_data] = Receiver(w_DecD, X, Mod);
[SNIR_DecD, BER_DecD] = ...
  Compute_SNIR_BER(w_DecD, X, X_des, DecD_dec_data, Des_User_Data, Mod);

SNIR_RAKE
SNIR_DecD
SNIR_STAR

% Plot the Constellation Diagrams of all Receivers
figure; 
subplot(311); plot(RAKE_plot,'o'); title('ST-RAKE Receiver');grid on; 
xlabel('Real Part'); ylabel('Imag Part'); axis([-2 2 -2 2]);
subplot(312);plot(DecD_plot,'o'); title('ST-Decorrelating Receiver');
grid on; xlabel('Real Part'); ylabel('Imag Part'); axis([-2 2 -2 2]);
subplot(313); plot(STAR_plot,'o'); title('STAR Subspace Receiver');grid on; 
xlabel('Real Part'); ylabel('Imag Part'); axis([-2 2 -2 2]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



function gold = Gold_Seq(pol_1,pol_2);  
%****************************************************
% Function to generate unbalanced Gold Seqs using Primitive Polynomials
% Written by Farrukh Rashid
% OUTPUT
% gold = Gold sequences for All users (Nc x NoOfGoldCodesofGivenLength) 
% INPUT
% pol_1 = Number of users
% pol_2 = Code Length
%****************************************************
[mSeq1] = genMSeq(pol_1, [1 1 1 1 1]);
[mSeq2] = genMSeq(pol_2, [1 1 1 1 1]);

N = length(mSeq1); % Period of the m sequence
gold=[mSeq1 mSeq2];

% Generate N Gold sequences using the shifted m sequences
for i=1:N,
    temp = genGSeq(mSeq1, shift(mSeq2, i-1)); 
    gold =[gold temp];
end


function [seq] = genMSeq(conns, initVal); 
%****************************************************
% Function to generate m-sequence by maximal length shift register
% Written by Farrukh Rashid
%****************************************************
conns = conns(1:length(conns)-1);         % discard C0
m     = length(conns);                    % length of shift reg   
L     = 2^m-1;                            % length of sequence  
reg   = initVal;                          % initial reg contents 
conns = fliplr(conns);                    % C1 along MSB reg bit is needed 
seq   = zeros(L,1);                       % for optimazation puposes
for var = 1:L
    temp     = mod( sum(conns & reg), 2 ); % calculate the ith value 
    seq(var) = temp;                      % store the ith value 
    reg      = [temp reg(1:m-1)];         % Shifting
end
seq   = ones(size(seq))-2*seq;            % Converting to bipolar
seq   = flipud(seq); 


function [gSeq] = genGSeq(seq1, seq2);
%****************************************************
% Function to generates a gold sequence from two input m-sequences
% Written by Farrukh Rashid
%****************************************************
seq1  = (1 - seq1)/2;            % Convert to bipolar
seq2  = (1 - seq2)/2;            % Convert to bipolar
gSeq  = mod((seq1 + seq2), 2);   % Modulo-2 addition to obtain  0,1.
gSeq  = ones(size(gSeq))-gSeq*2; % Convert to +1,-1 sequence from 0,1


function sseq = shift(seq, k);
%****************************************************
% Function to shift 'seq' down by 'k' bits
% Written by Farrukh Rashid
%****************************************************
L= length(seq);   % Find full length of the sequence
k = mod(k, L);    
sseq = [seq((L-k+1):L); seq(1:(L-k))];   % shift the sequence


function [Tx_Data_All_Users, Des_User_Data] = Tx_Data( M, Lw, Mod);
%****************************************************
% Function to Generate the Transmitted Data of all users
% Written by Farrukh Rashid
% OUTPUTS
% Tx_Data_All_Users,  = Matrix of Transmitted Data of all users (Lw+2 x M)
% Des_User_Data       = Binary Data Stream for the desired user (Lw+2 x 1)
% INPUTS
% M                   = Number of users
% Lw                  = Number of Symbols or Observation Interval
% Mod                 = Flag indicating the Modulation scheme(BPSK or QPSK)   
%****************************************************
Tx_Data_All_Users = []; 
for i = 1:M,
    if Mod == 'BPSK',
        Curr_bin_data = randint(Lw+2,1); % Generate BPSK binary data
    else
        Curr_bin_data = randint(2*(Lw+2),1); % Generate QPSK binary data
    end
    
    if i == 1, 
        Des_User_Data = Curr_bin_data;  % Store the desired user's data
    end
    
    Curr_sym_stream = Encoder(Curr_bin_data, Mod); % Encode into symbols        
    
    Tx_Data_All_Users  = [Tx_Data_All_Users Curr_sym_stream]; 
end


function dataout = Encoder(Bin_data, Mod)
%****************************************************
% Function to encode the data by BPSK or QPSK
% Written by Farrukh Rashid
% OUTPUT
% dataout = Vector of complex encoded data symbols (Lw+2 x 1)
% INPUT
% Bin_data = Binary data bits to be encoded (Lw+2 x 1 or 2*(Lw+2) x 1)
%****************************************************
if Mod == 'BPSK',
    dataout = 2*Bin_data - 1; % 0-> -1, 1 -> 1
else
    odd = Bin_data(1:2:end-1); even = Bin_data(2:2:end);
    odd_data = 2*odd - 1;    even_data = 2*even - 1;
    dataout =  odd_data + j*even_data;
end


function [X, X_des, Comp_Manifold_Matrix, H_1, Beta] = ...
    Received_Signal(M_1, drcts_1, l_1, Beta_1, gold, array, J, Ki, SNR);
%****************************************************
% Function to Model the first 3 terms of Rcvd Signal (Equation 10)
% Written by Farrukh Rashid
% OUTPUTS
% X          = Received Signal Modelled (2*N*Nc x Lw )
% X_des      = First Term of Eq 10 of [1]    (2*N*Nc x Lw )
% Comp_Manifold_Matrix   = Composite Channel Matrix for MU Decorr Receiver  
% H1         = Hi matrix for the Desired User
% Beta       = Beta vector for the Desired User 
% INPUTS
% M_1     = Data of All Users Tx 1            (Lw+2 x M)
% drcts_1 = Directions for all users     (sum(Ki) x 2)
% l_1     = Delays for all users         (sum(Ki) x 1)
% Beta_1  = Betas for all users          (sum(Ki) x 1)
% gold    = Gold codes for all users          (Nc x M)
% array   = Array Position Matrix             (N x 3)
% J       = Shift Matrix                      (2*Nc x 2*Nc)  
% Ki      = No. of multipaths for each user   (M x 1)
% SNR     = Input SNR in dB
%****************************************************

N = size(array,1);        % Array size
Nc = size(J,1)/2;         % Code size  
M = size(M_1,2);          % Number of users
No_symb = size(M_1,1);    % Number of data symbols

S_1 = spv(array, drcts_1) ;   % S_1 will be N x sum(Ki) matrix

Comp_Manifold_Matrix  = [];         % Stores the Composite Channel Matrix
X       = zeros(2*N*Nc, No_symb-2); % Stores the composite signal
X_des   = zeros(2*N*Nc, No_symb-2); % First Term of Eq 10 of [1]
Ki_indx = 1;                        % Index Flag for No. of multipaths
for i = 1:M, %% Loop for all users
    
    % Extract the parameter, data & gold code of the i'th user
    S_i_1 = S_1(:, Ki_indx:(Ki_indx + Ki(i)-1));
    l_i_1 = l_1(Ki_indx:(Ki_indx + Ki(i)-1),1);
    Beta_i_1 = Beta_1(Ki_indx:(Ki_indx + Ki(i)-1));
    data_1 = M_1(:,i);           
    c1 = [ gold(:,i); zeros(Nc,1)]; 
    
    % Collect the STAR Manifold Vectors of i'th user
    H_i_1 = []; H_pre_i_1 = []; H_next_i_1 = [];
    for k = 1:Ki(i), %% Loop for all multipaths of the current user
        H_i_1      = [ H_i_1  kron(S_i_1(:,k),(J^l_i_1(k))*c1)];	   	
    end

    % Matrix Hi for the previous and next symbol of i'th user