starchannelestimation.m

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

    odd = 0.5*(1+sign(real(Sym_in)));  even = 0.5*(1+sign(imag(Sym_in)));
    Dat_out(1:2:2*len) = odd;          Dat_out(2:2:2*len) = even;
end


function [w_STAR] = ...
              STAR_Eq_21_Ref_1_Weight(X,DOAs_est,TOAs_est,array,PN_1,B_1);   
%****************************************************
% Function to compute the Weight Vector for STAR Receiver (Eq 21 in [1])
% Written by Farrukh Rashid
% OUTPUTS
% w_STAR        = Weight Vector for the STAR Receiver   (2*Nc*N x 1)
% INPUTS
% X             = Received Signal                       (2*Nc*N x Lw)
% DOAs_est      = Vector containing Estimated DOAs      (K_des x 1)
% TOAs_est      = Vector containing Estimated TOAs      (K_des x 1)
% B_1           = Desired user's Fading Coefficients    (K_des x 1)
%****************************************************
Lw = size(X,2);   % Total number of symbols/observation interval
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]
c1 = [PN_1; zeros(Nc,1)];

% Compute the Desired user's STAR manifold matrix   (2*N*Nc x K_1)
H_1 = [];
for k = 1:length(DOAs_est), %% Loop for all multipaths of the current user
    H_1      = [ H_1  kron(spv(array,[DOAs_est(k) 0]),(J^TOAs_est(k))*c1)];
end
    
Rxx = X * X'/Lw;  % Received Signal Covariance Matrix
R_Sig = H_1 * B_1 * B_1' * H_1'; % Desired Signal Covariance Matrix
R_N_I = Rxx - R_Sig; % Noise/Interference Covariance Matrix (Eq 22 of [1])

[E,D] = eig(R_N_I); Dia = diag(D);
%% Sort the Eigenvalues and their corresponding eigenvectors
[Dia_sorted, I] = sort(Dia);
Dia_sorted = flipud(Dia_sorted);
I = flipud(I);
for i = 1 : length(Dia),
    E_sorted(:,i) = E(:,I(i));
end

% Identify the Size of Signal Subspace by taking derivative
der = diff( (abs(Dia_sorted)) );
I = min(find(abs(der) < 100));
SizeofSignalSubspace = I-1;

% Noise Subspace of R(n + MAI/ISI) Matrix
En = E_sorted(:, SizeofSignalSubspace+1:size(E));
Pn = En*pinv(En'*En)*En';

w_STAR = Pn*H_1*pinv(H_1'*Pn*H_1)*B_1; %Equation 21 in [1]
w_STAR = w_STAR/norm(w_STAR) ; 


function [w_RAKE] = RAKE_Eq_22_Ref_2_Weight(H_1, B_1);
%****************************************************
% Function to compute the ST-RAKE Receiver Weight (Eq 22 in [2])
% Written by Farrukh Rashid
% OUTPUTS
% w_RAKE        = Weight Vector for the RAKE Receiver (2*Nc*N x 1)
% INPUTS
% H_1           = Desired user's STAR manifold matrix (2*N*Nc x K_1)
% B_1           = Desired user's Fading Coefficients  (K_1 x 1)
%****************************************************
w_RAKE = H_1* B_1; %  RAKE Weight vector for user 1, Tx Ant 1 (2*N*Nc x 1)
w_RAKE = w_RAKE / norm(w_RAKE);


function [w_DecD] = DecD_Eq_23_Ref_2_Weight(Comp_Manifold_Matrix);
%****************************************************
% Function to compute the ST-Decorrelating Detector Weight( Eq 23 in [2])
% Written by Farrukh Rashid
% OUTPUTS
% w_DecD = Weight Vector for the ST-Decorrelating Receiver (2*Nc*N x 1)
% INPUTS
% Comp_Manifold_Matrix = Composite Manifold Matrix (2*Nc*N x 3*M)
%****************************************************
Comp_1 = pinv(Comp_Manifold_Matrix);
w_DecD = Comp_1(2,:)';
w_DecD = w_DecD/norm(w_DecD);


function [Sym_Rcvd, dec_data] = Receiver(Rx_wt, X, Mod);
%****************************************************
% Function to Impelement the Receiver (Eqns 20 and 23 in [1])
% Written by Farrukh Rashid
% OUTPUTS
% Sym_Rcvd      = Normalized Received Symbols at the output of the receiver 
% dec_data      = Decoded binary data stream at the output of the receiver
% INPUTS
% Rx_wt         = Receiver Weight (2*N*Nc x 1)
% X             = Received Signal (2*N*Nc x Lw)
% Mod           = Modulation scheme used (BPSK or QPSK)
%****************************************************
Sym_Rcvd = Rx_wt'* X;               % Eq 23 in [1]
dec_data = Decoder(Sym_Rcvd, Mod);  % Eq 20 in [1]
scale  = mean(abs(Sym_Rcvd)); Sym_Rcvd = Sym_Rcvd / scale; % Normalize


function [SNIR, BER] = ...
 Compute_SNIR_BER(Rx_wt, X, X_des, dec_data, Des_User_Data, Mod);
%****************************************************
% Function to Compute the SNIR and BER at the Rx Output 
% Written by Farrukh Rashid
% OUTPUTS
% SNIR          = SNIR at the output of the receiver 
% BER           = BER at the output of the receiver
% INPUTS
% Rx_wt         = Receiver Weight (2*N*Nc x 1)
% X             = Received Signal (2*N*Nc x Lw)
% X_des         = Desired User's signal (First Term of Eq 10 of [1])
% dec_data      = Decoded binary stream (Lw x 1 or 2*Lw x 1)
% Des_User_Data = Desired User's Binary Stream (Lw+2 x 1 or 2*(Lw+2) x 1)
% Mod           = Modulation scheme used (BPSK or QPSK)
%****************************************************
Lw = size(X,2);   % Total number of symbols/observation interval

Rxx = X * X'/Lw;  % Received Signal Covariance Matrix

R_Sig = X_des*X_des'/Lw; % Desired Signal Covariance Matrix

R_N_I = Rxx - R_Sig; % Noise/Interference Covariance Matrix (Eq 22 of [1])

SNIR = 10*log10(abs((Rx_wt'*R_Sig*Rx_wt)/(Rx_wt'*R_N_I*Rx_wt)));

% Strip off 1st and last symbol
if Mod == 'BPSK',
    Des_Data = (Des_User_Data(2:(end-1))).'; 
else
    Des_Data = (Des_User_Data(3:(end-2))).'; 
end

BER = length(find((Des_Data == dec_data)~=1)); 


function S=spv(array,directions,mainlobe);
%****************************************************
% S=spv(array,directions,mainlobe);
% Source Position Vectors:
% NB: if f and c are given then use array*2*f/c;
% written by Prof A.Manikas
% Modified	 5.6.01 - Jason W.P. Ng
%****************************************************
if nargin<=2, mainlobe=[]; end;

SOURCES=frad(directions);
N=length(array(:,1));
M=length(SOURCES(:,1));

if nargin<=2 | isempty(mainlobe)
   U0=zeros(N,M);
else
   saz=frad(mainlobe(1));sel=frad(mainlobe(2));
   k0 = fki(saz,sel);
   U0=  repc(array*k0,M);
end;

KI = fki(SOURCES(:,1),SOURCES(:,2));

S = exp(-j*(array*KI -U0));


function x=frad(degrees);
%****************************************************
% written by Prof A.Manikas
%degrees to rad
%****************************************************
x=degrees*pi/180.0;


function KI=fki(az,el);
%****************************************************
% written by Prof A.Manikas
%wavenumber vector in half wavelengths;
%****************************************************
KI = pi*[cos(az).*cos(el)   sin(az).*cos(el)    sin(el)]';


function B=repc(A,m);
%****************************************************
% written by Prof A.Manikas
%Duplicates m columns of column vector A
%****************************************************
B=A*ones(1,m);


function [Mmdl,Maic,MDL,AIC] = detect(Rxx,L)
%*************************************************************************
% by Naushad Dowlut, 8 Dec 1993
% Reference : 	M. Wax and T. Kailath
%		"Detection of Signals by Information Theoretic Criteria"
%		IEEE Trans. ASSP, vol. ASSP-33, pp. 387-392, Apr. 1985              
%----------------------------------------------------------------------
% Estimates the number of sources using the AIC and MDL criteria
%	Mmdl, Maic - no. of sources as estimated by AIC and MDL
%	MDL, AIC   - vector of criteria values
%	Rxx        - practical data covariance matrix
%	L 	   - no. of snapshots
%*************************************************************************

N = size(Rxx,1); 		% no. of sensors

eigvals = sort(abs(eig(Rxx)));	% sorted in ascending order

cp = flipud(log(cumprod(eigvals)));
cs = flipud(log(cumsum(eigvals)));

M = [0:(N-1)]';		% M - parameter used for the estimation 
r = N - M;


AIC = zeros(N,1);
MDL = zeros(N,1);

AIC = -2*L*(cp + r.*(log(r) - cs)) + 2*M.*(2*N-M);
MDL = -L*(cp + r.*(log(r) - cs)) + (1/2)*log(L)*M.*(2*N-M);
AIC=1./AIC;MDL=1./MDL;
[Y,Maic] = max(AIC);
[Y,Mmdl] = max(MDL);

Maic = Maic-1;		% since M starts @ 0
Mmdl = Mmdl-1;


function best = minmatr(Z, AZarea, ELarea, n)
%****************************************************
%  best = minnmatr(Z, AZarea, ELarea, n)
%  The Lowest n-Mins of a Matrix Z
%  Z      ... Result of MUSIC
%  AZarea ... x axis (Azimuth Area)
%  ELarea ... y axis (Elevation area)
%  n      ... The number of minimum points required
%  Example best = minnmatr(Z, 1:1:180, 0:5:90, 3)
%
%  Digital Communications Section, EEED, IC
%  Written by : Nobuyuki Takai
%  Supervisor : Dr A. Manikas
%  18 April 1994
%****************************************************

Z     = Z';best  = [];
ON    = 1;
OFF   = 0;
flagN = OFF;
flagS = OFF;
flagW = OFF;
flagE = OFF;
[lenaz,lenel] = size(Z);
pits=[];

for el=1:lenel
    for az=1:lenaz
%        home; fprintf('az=%f\nel=%f\n',AZarea(az),ELarea(el));
                    C = Z(az,el);
        if az>1     N = Z(az-1,el); end;
        if az<lenaz S = Z(az+1,el); end;
        if el>1     W = Z(az,el-1); end;
        if el<lenel E = Z(az,el+1); end;
% North
        if az>1
            if C<=N
                flagN = ON;
            end;
        else
            flagN = ON;
        end;
% South
        if az<lenaz
            if C<=S
                flagS = ON;
            end;
        else
            flagS = ON;
        end;
% WZ
        if el>1
            if C<=W
                flagW = ON;
            end;
        else
            flagW = ON;
        end;
% East
        if el<lenel
            if C<=E
                flagE = ON;
            end;
        else
            flagE = ON;
        end;
% Decision
        if (flagN)&(flagS)&(flagW)&(flagE) 
            pits = [pits; [az el]];
        end;
% Initialize Flags
        flagN = OFF;
        flagS = OFF;
        flagW = OFF;
        flagE = OFF;
    end;
end;
if isempty(pits)
 best=[];
 return;
end 
%vals=zeros(1,length(pits(:,1)));
for a=1:length(pits(:,1))
    vals(a) = Z(pits(a,1),pits(a,2));
end;

[sorted index] = sort(vals);
if n<length(index)
    index  = index (1:n);
    sorted = sorted(1:n);
end;
pointer = pits(index,:);

az   = AZarea(pointer(:,1));
el   = ELarea(pointer(:,2));

best = [az' el' sorted'];