starchannelestimation.m

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

    H_pre_i_1  = (kron(eye(N), (J')^Nc))*H_i_1;  
    H_next_i_1 = (kron(eye(N),J^Nc))*H_i_1;    
    
    % Product of Hi and Beta vector of i'th user
    HB_pre = H_pre_i_1*Beta_i_1;  
    HB = H_i_1*Beta_i_1;  
    HB_nex = H_next_i_1*Beta_i_1;
    
    %% Arrange the data stream
    data_pre = M_1(1:No_symb-2,i).'; %% a_i[n-1]
    data_act = M_1(2:No_symb-1,i).'; %% a_i[n]
    data_nex = M_1(3:No_symb,i).';   %% a_i[n+1]
    
    X_i =   HB * data_act;
    ISI_i = [HB_pre  HB_nex]*[data_pre; data_nex]; 
        
    X = X + X_i + ISI_i;
    
    % Get the clean signal and RAKE weight for the Desired (1st) User
    if i == 1,
        Beta    = Beta_i_1;  % Store the Betas
        H_1     = H_i_1;     % Store the Hi
        X_des = X_i;         % First Term of Eq 10 of [1]
    end
    
    % Compute the Composite Channel Matrix for MU DECORR Receiver 
    Comp_Manifold_Matrix = [Comp_Manifold_Matrix HB_pre  HB  HB_nex]; 
    
    Ki_indx = Ki_indx + Ki(i); % Update Index Flag for No. of multipaths
end

% Average Desired User's signal power per symbol
Avg_Des_Sig_power=(1/(No_symb-2))*trace(X_des'*X_des)  
% Desired User's Signal power per sample
Des_Sig_power = Avg_Des_Sig_power/(2*N*Nc); 

SNR_abs = 10^(SNR/10);                   % SNR in absolute units
Noise_power =  Des_Sig_power / SNR_abs;  % Noise Power per sample 

noise =(sqrt(Noise_power)/sqrt(2))*...
                    (randn(2*N*Nc,No_symb-2)+j*randn(2*N*Nc,No_symb-2));

Avg_Noise_pow = (1/(No_symb-2))*trace(noise'*noise)  % Average Noise Power

Avg_SNR = 10*log10(Avg_Des_Sig_power/Avg_Noise_pow)  % Average SNR

X     = X + noise;		% Add Noise



function [DOAs_est, TOAs_est, MUSIC_Spec] = ...
                  STAR_Channel_Estimation(x_final, c1, Rx_array, Q, K_des);
%****************************************************
% Function to Perform the STAR Channel Estimation 
% Written by Farrukh Rashid
% OUTPUTS
% DOAs_est = Vector containing Estimated DOAs (K_des x 1)
% TOAs_est = Vector containing Estimated TOAs (K_des x 1)
% INPUTS 
% x_final  = Received Signal             (2*N*Nc x Lw)
% c1       = Ref Temp Manif Vector       (2*Nc x 1)
% Rx_array = Array Position Matrix
% Q        = Number of sumatrices to be averaged in temporal smoothing 
% Note that Q > K_des so that d < 2*Nc (we choose Q = 10 for simulation)
% K_des    = Number of multipaths of the desired user
%****************************************************
Lw = size(x_final,2);     % Number of Symbols
Nc = length(c1)/2;        % Code Length
N = size(Rx_array, 1);    % Number of array elements

%% STEP 1: PASS THE RECEIVED VECTOR THROUGH PREPROCESSOR (Equation 15)
[z1, D1, phi_vec] = preprocessor(x_final, c1, Nc, N); 

%% STEP 2: PERFORM SPATIAL SMOOTHING ON TOP OF TEMPORAL SMOOTHING
[R_z_STsmooth] = Spatial_Temporal_Smoothing(z1, N, Q); 

%% Spatial & Temporal Smoothing the Diagonal matrix for Generalized EVs 
[D1_STsmooth] = smoothing_Diagonal(D1, N, Q); 

d = size(R_z_STsmooth,1)/(N-1); %Length of subvector with Spatial Smothing 
%d=size(R_z_STsmooth1,1)/N;%Length of subvector without Spatial Smoothing 
phi_vec_d (:,1) = phi_vec(1:d,1); % Sub vector of phi_vec with length d       

%% Detect the size of signal subspace
[Mmdl,Maic,MDL,AIC] = detect(R_z_STsmooth, Lw);
SizeofSignalSubspace=Maic;

%% STEP 3: APPLY THE MUSIC-type COST FUNCTION (Equation 19)
MUSIC_Spec=music(R_z_STsmooth,Rx_array(1:4,:),...
                           SizeofSignalSubspace,0:180,0,0:Nc-1,phi_vec_d);

%% STEP 4: PLANAR SEARCH TO FIND DOA/TOA
best = minmatr(-1*MUSIC_Spec,0:180,0:Nc-1,K_des);
DOAs_est = best(:,1);
TOAs_est = best(:,2);

%% PLOTS FOR MUSIC SPECTRUM 
figure;
%colormap(cool);
surf([0:1:180],[0:1:(Nc-1)], MUSIC_Spec);
shading('interp');
xlabel('DOA (Deg)');
ylabel('TOA (Tc)');
zlabel('dB');

figure;
%colormap(cool);
contour([0:1:180],[0:1:(Nc-1)], MUSIC_Spec,40);
grid on;
xlabel('DOA (Deg)');
ylabel('TOA (Tc)');


function [z1, D, phi_vec] = preprocessor(X, c1, Nc, N);
%****************************************************
% Function to Preprocess the Received Signal (Equation 15)
% Written by Farrukh Rashid
% OUTPUTS
% z1= Preprocessed Signal (2*N*Nc x Lw)
% D = Spatio Temporal Smoothed Diagonal Matrix Z1*Z1' (2*N*Nc x 2*N*Nc)
% phi_vec = phi vector to be used in Music3D function (2*N*Nc x 1)

% INPUTS
% X = Received Signal to be preprocessed (2*N*Nc x Lw)
% c1 = Ref Temp Manif Vector (2*N*Nc x 1)
% Nc = Code Size 
% N = No of Array elements
%****************************************************
phi=exp(-j*(2*pi)/(2*Nc));          % Define the Twiddle Factor
phi_vec = [phi.^([0:(2*Nc-1)]')];   % Define the Basis Vector for FT

F=[]; % Fourier Transformation Matrix (Equation 11)
for i = 0 : 2*Nc-1,
    F=[F phi_vec.^i];
end

% Spatio Temporal Preprocessor (2*N*Nc x 2*N*Nc) (Equation 13)
Z1=kron(eye(N), ( inv( diag(F*c1) ) )*F ); 

z1=Z1*X; % Applying the preprocessor to the received signal (Equation 15)

D=Z1*Z1'; % Spatio Temporal Smoothed Diagonal Matrix


function [R_z_STsmooth] = Spatial_Temporal_Smoothing(z1, N, Q);
%****************************************************
% FUNCTION TO OVERLAY SPATIAL SMOOTING ON TOP OF TEMPORAL SMOOTHING
% OUTPUT
% Written by Farrukh Rashid
% OUTPUT
% R_z_STsmooth = (N-1)*d x (N-1)*d Spatial Smoothed Covariance Matrix 
% INPUTS
% z1 = (2*N*Nc x Lw) matrix
% Q = Number of sumatrices to be averaged, Q > Ki so that d < 2*Nc 
%****************************************************
Nc = size(z1,1) /(2*N);  % Code Size

Rz1 = (1/size(z1,2))*z1*z1';

%% FIRST PERFORM SPATIAL SMOOTHING
% R_z_Ssmooth is (N-1)*2*Nc x (N-1)*2*Nc Spatial Smoothed Covariance Matrix
[R_z_Ssmooth ] = spat_smooth_over (Rz1, 2, N); 

%% THEN PERFORM TEMPORAL SMOOTHING (Equation 18)
% R_z_STsmooth is (N-1)*d x (N-1)*d
[ R_z_STsmooth ] = temp_smooth_matrix(R_z_Ssmooth, Nc, Q); 



function [R_z_Ssmooth ] = spat_smooth_over (R_z, Q, N );
%****************************************************
% Function to do Spatial smoothing overlaid on Temporal Smoothing
% Written by Farrukh Rashid
% INPUTS
% R_z  = Covariance matrix which is to be spatially smoothed (N*d x N*d) 
% Q    = Number of submatrices for Spatial Smoothing (=2 in our case)
% N    = Number of antenna elements (5 in our case)

% OUTPUTS
% R_z_Ssmooth = Covariance of the spatially smoothed preprocessed signal 
%****************************************************
d = size(R_z, 1)/(N);       % Size of R_z is N*d  
m = d *(N-1);               

R_z_Ssmooth = zeros( m, m); 
for i = 0 : (Q-1), % Loop for 2 subarrays
    %% First Form the submatrix
	R_sub(:,:) = R_z( ((i*d) + 1) : m, ((i*d) + 1) : m  );
    R_z_Ssmooth = R_z_Ssmooth + R_sub;    
    m = m + d;
end

R_z_Ssmooth = R_z_Ssmooth/Q;  


function [R_z_Tsmooth ] = temp_smooth_matrix(R_mat, Nc, Q);
%****************************************************
% Function to do temporal smoothing on a Matrix  (Equation 18)
% Written by Farrukh Rashid
% INPUTS
% R_mat = Covariance Matrix which is to be smoothed 
% Nc = Code Size (31)
% Q = Number of submatrices ( >= Ki)

% OUTPUTS
% R_z_Tsmooth = Covariance Matrix
%****************************************************
d = 2*Nc - 1 - Q;    % Length of the overlapping submatrices 
N = size(R_mat,1)/(2*Nc); % (N = 5)

Num_segment = N*N; % Total number of segmented matrices 

for i = 0 : N - 1 , % Parse the Rows of R_mat (5 rows)
    for j = 0: N - 1, % Parse the columns of R_mat (5 columns)
        
        %% Extract the segmented matrix of size 2*Nc
     R_seg_mat(:,:)=R_mat((1+i*2*Nc):((i+1)*2*Nc),(1+j*2*Nc):((j+1)*2*Nc));
        
        Rz_acc_mat = zeros(d,d); 
        for k = 1:Q, % For each segment parse Average a set of Q submarices
            
            %% Extract the Submatrices (N_1 x d)
            Rz_sub_mat(:,:) = R_seg_mat( k:(k+d-1), k:(k+d-1)); 

            %% Accumulate the Submatrices
            Rz_acc_mat = Rz_acc_mat + Rz_sub_mat;
        end

        % Divide by the number of submatrices
        Rz_acc_mat = Rz_acc_mat/Q; 
        
        R_z_Tsmooth((1+i*d):((i+1)*d),(1+j*d):((j+1)*d))=Rz_acc_mat(:,:); 
    end  
end


function [D1_STsmooth, D2_STsmooth] = smoothing_Diagonal(D1, N, Q);
%****************************************************
% FUNCTION TO do the SPATIAL SMOOTING ON TOP OF TEMPORAL SMOOTHING for the
% ST Diagional Matrix
% Written by Farrukh Rashid
% OUTPUTS
% D1_STsmooth = (N-1)*d x (N-1)*d Spatial Smoothed Covariance Matrix 
% Q = Number of sumatrices to be averaged, Q > Ki so that d < 2*Nc 

% INPUT
% D1 = (2*N*Nc x 2*N*Nc) matrix
%****************************************************
Nc = size(D1,1) /(2*N);  % Code Size (Nc = 31)

%% FIRST PERFORM TEMPORAL SMOOTHING
[ D1_Tsmooth ] = temp_smooth_matrix(D1, Nc, Q); 

%% THEN PERFORM SPATIAL SMOOTHING
% (N-1)*d x (N-1)*d Spatial Smoothed Covariance Matrix
[D1_STsmooth ] = spat_smooth_over (D1_Tsmooth, 2, N); 


function [Z]=music ...
(Rxx,array,M,AZarea,ELarea,TOA_Area, phi_vec_d,mainlobe,ingain,inphase);
%****************************************************
% Z=music(Rxx,array,M,AZarea,ELarea,phi_vec_d, mainlobe,ingain,inphase);
% writen by Dr A.Manikas (IC)
% modified by Farrukh Rashid (September 2005)
% estimates the MuSIC spectrum
% ref.:
% N.B.: if f and c are given then  use array*2*f/c
%****************************************************
N=length(array(:,1));

if nargin<9; g=ones(N,length(AZarea));
      else g=repc(ingain.*exp(j*inphase),length(AZarea));
;end;
if nargin<8
   mainlobe=[];
end;

[MEIG,D] = eig(Rxx);
[lamda,k]=sort(diag(D));
MEIG=MEIG(:,k);

EE = MEIG(:,1:length(Rxx(:,1))-M);
for I=1:N-M;
    EE(:,I)=EE(:,I)/sqrt(abs(EE(:,I)'*EE(:,I)));
end;

for el=ELarea;
      SOURCES=[AZarea',el*ones(size(AZarea'))];
      for k = TOA_Area,
          S=kron(g.*spv(array,SOURCES,mainlobe), phi_vec_d.^k);
          Z(k+1,:)=-10*log10(real(diag(S'*EE*EE'*S)))';
     end; 
end;
  

function B_est = ...
            EstimateFadingCoeff(DOAs_est,TOAs_est,DOA_orig,TOA_orig,B_1);
%****************************************************
% Function to estimate the fading coefficients from DOA/TOA channel 
% estimates (Here we assume that Fading Coefficients have already been 
% computed by using another algorithm & we just re-order them according to 
% DOA/TOA estimates)
% Written by Farrukh Rashid
%****************************************************
for k = 1:length(DOAs_est),
    B_est(k,1) = 0;
    for cnt = 1:length(DOAs_est),
        if DOA_orig(cnt)==round(DOAs_est(k)) & ...
                                    TOA_orig(cnt)== round(TOAs_est(k)),
            B_est(k,1) = B_1(cnt);
        end
    end
    
end


function Dat_out = Decoder(Sym_in, Mod)
%****************************************************
% Function to decode the BPSK or QPSK symbols (Eq 20 in [1])
% Written by Farrukh Rashid
% OUTPUT  
% Dat_out = Vector of binary decoded data (Lw x 1 or 2*Lw x 1)
% INPUT
% Sym_in = Vector of symbols to be decoded (Lw x 1)
%****************************************************
len = length(Sym_in);
if Mod == 'BPSK',
    Dat_out = 0.5*(1+sign(real(Sym_in))); 
else