ofdm_idma_ray_perfect_comp.m

这是一个关于未编码IDMA系统的仿真程序

% This is a simulation program for uncoded IDMA system with BPSK modulation
% and real channel.

clear all;
clc;

rand('state',54321);
randn('state',12345);

%********************** idma parameters ***************************
NUser = 1;
DataLen	= 256;
SpreadLen = 16;
ChipLen	= ( DataLen .* SpreadLen );
ItNum = 10;
Block = 40;
EbNoNum = 4;
EbNoStart = 0;
EbNoStep = 3;
% STEP(x)		( (x) > 0 ? 0 : 1 )

%********************** ofdm parameters ***************************

para=128;   % Number of parallel channel to transmit (points)
fftlen=128; % FFT length
noc=128;    % Number of carrier
ml=2;       % Modulation level : QPSK
nd=ChipLen./ml./para;       % Number of information OFDM symbol for one loop
br=256000;  % Bit rate per carrier
sr=br./ml;  % QPSK symbol rate
sr_ofdm=sr./para;  % OFDM symbol rate
gilen=32;   % Length of guard interval (points)

%********************** rayleigh parameters ***********************
tstp=1./sr_ofdm./SpreadLen./(fftlen+gilen);
itn_add=nd.*(fftlen+gilen).*100;  % Update time
itn=1000;
fd=200;
nowave=6;
% itn=1000;
flat=1;

Data = zeros(NUser,DataLen);
Data_RS = zeros(NUser,DataLen);
% Data_RP = zeros(NUser,DataLen);
ScrambleRule = zeros(NUser,ChipLen);
% Ber = zeros(1,EbNoNum);
H = zeros(1,NUser);
H2 = zeros(1,NUser);	% H2[.] = H[.] * H[.]%
Chip = zeros(NUser,ChipLen);
Tx = zeros(NUser,ChipLen);
%    Tx_i = zeros(NUser,(ChipLen./2));
%    Tx_q = zeros(NUser,(ChipLen./2));
Rx = zeros(NUser,ChipLen);
   Rx_i = zeros(1,(ChipLen./2));
   Rx_q = zeros(1,(ChipLen./2));
SpreadSeq = zeros(1,SpreadLen);
TotalMean = zeros(1,ChipLen);
TotalVar = zeros(1,ChipLen);
Mean = zeros(NUser,ChipLen);
Var = zeros(NUser,ChipLen);
AppLlr = zeros(NUser,DataLen);
Ext = zeros(NUser,ChipLen);

%************************* Energy Profile *********************************
% All users use the same energy level, i.e. no energy allocation.
for ( nuser_EP = 1:NUser )
    H2( nuser_EP ) = 1;
end
% 	 The following energy levels are for the uncoded IDMA system wih
% 	 BPSK modulation and real channel. The parameters of the system
% 	 are: No. of user = 64 and spreading length =16.
% 	
% 	 For details of the energy allocation technique, please refer to
% 	 this paper: Li Ping and Lihai Liu, "Analysis and design of IDMA
% 	 systems based on SNR evolution and power allocation," in Proc.
% 	 VTC' 2004-Fall, Los Angeles, CA. Sept. 2004 or 
% 	 http://www.ee.cityu.edu.hk/~liping/Research/Conference/IDMA-PA.pdf
% 	
if( ( NUser == 64 ) && ( SpreadLen == 16 ) )
    H2(64)	=   1.0000; %  0.000dB
    H2(63)	=   1.0000; %  0.000dB
	H2(62)	=   1.0000; %  0.000dB
	H2(61)	=   1.0000; %  0.000dB
	H2(60)	=   1.0000; %  0.000dB
	H2(59)	=   1.0000; %  0.000dB
	H2(58)	=   1.0000; %  0.000dB
	H2(57)	=   1.0000; %  0.000dB
	H2(56)	=   1.0000; %  0.000dB
	H2(55)	=   1.0000; %  0.000dB
	H2(54)	=   1.0000; %  0.000dB
	H2(53)	=   1.0000; %  0.000dB
	H2(52)	=   1.0000; %  0.000dB
	H2(51)	=   1.0000; %  0.000dB
	H2(50)	=   1.0000; %  0.000dB
	H2(49)	=   1.0000; %  0.000dB
	H2(48)	=   1.0000; %  0.000dB
	H2(47)	=   1.0000; %  0.000dB
	H2(46)	=   1.0000; %  0.000dB
	H2(45)	=   1.0000; %  0.000dB
	H2(44)	=   1.0000; %  0.000dB
	H2(43)	=   1.0000; %  0.000dB
	H2(42)	=   1.0000; %  0.000dB
	H2(41)	=   1.0000; %  0.000dB
	H2(40)	=   6.2748; %  7.976dB
	H2(39)	=   6.2748; %  7.976dB
	H2(38)	=   6.2748; %  7.976dB
	H2(37)	=   6.2748; %  7.976dB
	H2(36)	=   6.2748; %  7.976dB
	H2(35)	=   6.2748; %  7.976dB
	H2(34)	=   6.2748; %  7.976dB
	H2(33)	=   6.2748; %  7.976dB
	H2(32)	=   6.2748; %  7.976dB
	H2(31)	=   6.2748; %  7.976dB
	H2(30)	=   6.2748; %  7.976dB
	H2(29)	=   6.2748; %  7.976dB
	H2(28)	=   8.1583; %  9.116dB
	H2(27)	=   8.1583; %  9.116dB
	H2(26)	=  23.2970; % 13.673dB
	H2(25)	=  23.2970; % 13.673dB
	H2(24)	=  23.2970; % 13.673dB
	H2(23)	=  23.2970; % 13.673dB
	H2(22)	=  23.2970; % 13.673dB
	H2(21)	=  23.2970; % 13.673dB
	H2(20)	=  23.2970; % 13.673dB
	H2(19)	=  23.2970; % 13.673dB
	H2(18)	=  30.2901; % 14.813dB
	H2(17)	=  30.2901; % 14.813dB
	H2(16)	=  30.2901; % 14.813dB
	H2(15)	=  30.2901; % 14.813dB
	H2(14)	=  86.4968; % 19.370dB
	H2(13)	=  86.4968; % 19.370dB
	H2(12)	=  86.4968; % 19.370dB
	H2(11)	=  86.4968; % 19.370dB
	H2(10)	=  86.4968; % 19.370dB
	H2(9)	=  86.4968; % 19.370dB
	H2(8)	= 112.4605; % 20.510dB
	H2(7)	= 112.4605; % 20.510dB
	H2(6)	= 112.4605; % 20.510dB
	H2(5)	= 112.4605; % 20.510dB
	H2(4)	= 112.4605; % 20.510dB
	H2(3)	= 112.4605; % 20.510dB
	H2(2)	= 112.4605; % 20.510dB
	H2(1)	= 112.4605; % 20.510dB
end

% Normalization of the energy level of all users.
tmp_EP = 0;
for( nuser_EP = 1:NUser )
    tmp_EP = tmp_EP + H2(nuser_EP);
end
for( nuser_EP=1:NUser )
    H2(nuser_EP) = NUser .* H2(nuser_EP) ./ tmp_EP;
end

% Produce the amplitude scaling level
for( nuser_EP = 1:NUser )
    H(nuser_EP) = sqrt( H2(nuser_EP) );
end

% Generate the spreading sequence: {+1, -1, +1, -1, ... }.
for( i=1:SpreadLen )
    SpreadSeq(i) = 2 .* mod(i,2) - 1;
end

% Generate the interleaving index.
for( nuser=1:NUser ) 
	for( i=1:ChipLen )
        ScrambleRule(nuser,i) = i;
    end
	for( i=1:(ChipLen - 1) ) 
		j = i + mod(( randint(1,1,32768) .* randint(1,1,32768) ),( ChipLen - i ));
		tmp = ScrambleRule(nuser,i);
		ScrambleRule(nuser,i) = ScrambleRule(nuser,j);
		ScrambleRule(nuser,j) = tmp;
    end
end

%************************* The Simulation Process *************************
for( i=1:EbNoNum )
	ebno = EbNoStart + (i-1) .* EbNoStep;
%     EsNo = ebno + 10 .* log10(2);
% 	snr = 10^( EsNo ./ 10 ) ./ SpreadLen;
% 	sigma = sqrt( 0.5 ./ snr );
    snr = 10^( ebno ./ 10 ) ./ SpreadLen;
% 	sigma = sqrt( 0.5 .*0.5./ snr );
    sigma = sqrt( 0.5 ./ snr );
	error = 0;
    itn=1000;
	for( block=1:Block )
	
%******************************* Transmitter( ) ***************************
% int i, j, m, nuser;
% double tmp;

% Data Generation.
for( nuser_Tr=1:NUser )
    for( i_Tr=1:DataLen )
        Data(nuser_Tr,i_Tr) = mod( randint(1,1,32768) , 2 );
    end
end
% Spreading process.
for( nuser_Tr=1:NUser )    
    m_Tr = 1;
    for( i_Tr=1:DataLen )
        tmp_Tr = 1 - ( 2 .* Data(nuser_Tr,i_Tr) ) ;
        for( j_Tr=1:SpreadLen )
            Chip(nuser_Tr,m_Tr) = tmp_Tr .* SpreadSeq(j_Tr);
            m_Tr = m_Tr + 1;
        end
    end
end

% Interleaving process.
for( nuser_Tr=1:NUser )
    for( i_Tr=1:ChipLen )
        Tx(nuser_Tr,i_Tr) = Chip(nuser_Tr,ScrambleRule(nuser_Tr,i_Tr));
    end
end

% ich4=zeros(1,fftlen.*nd./ml);    %  if without GI
% qch4=zeros(1,fftlen.*nd./ml);

ich4=zeros(1,(fftlen+gilen).*nd);    %  if within GI
qch4=zeros(1,(fftlen+gilen).*nd);
for( nuser_Tr=1:NUser )
    
%****************** Serial to parallel conversion ***********************

Tx_para=reshape(Tx(nuser_Tr,:),para,nd.*ml); %  reshape : built in function

%************************** QPSK modulation ***************************** 

% kmod = 1/sqrt(2);
Tx_i = H(nuser_Tr) .* Tx_para(:,1:2:end);% .* kmod;
Tx_q = H(nuser_Tr) .* Tx_para(:,2:2:end);% .* kmod;

% [ich,qch]=qpskmod(paradata,para,nd,ml);
% kmod=1/sqrt(2); %  sqrt : built in function
% ich1=ich.*kmod;
% qch1=qch.*kmod;

%******************* IFFT ************************

tx_ifft=complex(Tx_i,Tx_q);
ty=ifft(tx_ifft);      %  ifft : built in function
ich2=real(ty);   %  real : built in function
qch2=imag(ty);   %  imag : built in function

%********* Gurad interval insertion **********

[ich3(nuser_Tr,:),qch3(nuser_Tr,:)]= giins(ich2,qch2,fftlen,gilen,nd);
fftlen2=fftlen+gilen;

%********* OFDM symbols in total **********

ich4=ich4+ich3(nuser_Tr,:);
qch4=qch4+qch3(nuser_Tr,:);
end

% %******************************* Channel_AWGN( sigma ) ********************
% % 	int i, nuser;
% %   separate into i and q 
% 
% for( i_ChAWGN = 1:(ChipLen./2) )
% 	% Additive white Gaussian noise.
%     Rx_i(i_ChAWGN) = sigma .* randn();
%     Rx_q(i_ChAWGN) = sigma .* randn();
% 
% 	% Multiple access channel and energy scaling.
% 	for( nuser_ChAWGN=1:NUser )
%         Rx_i(i_ChAWGN) = Rx_i(i_ChAWGN) + H(nuser_ChAWGN) .* Tx_i(nuser_ChAWGN,i_ChAWGN);
%         Rx_q(i_ChAWGN) = Rx_q(i_ChAWGN) + H(nuser_ChAWGN) .* Tx_q(nuser_ChAWGN,i_ChAWGN);
%     end
% end

%********* Attenuation Calculation *********

spow=sum(sum(ich3.^2+qch3.^2))./NUser./nd./fftlen2./ml;;  %  sum : built in function
attn=0.5*spow./snr;
attn=sqrt(attn);

%******************************* Channel_RaySingle(  ) ********************
for( nuser_Tr=1:NUser )
    [ich33(nuser_Tr,:),qch33(nuser_Tr,:),ramp(nuser_Tr,:)]=fade(ich3(nuser_Tr,:),qch3(nuser_Tr,:),fftlen2.*nd,tstp,fd,nowave,itn,flat);
    itn=itn+itn_add;
end
ich44=sum(ich33,1);
qch44=sum(qch33,1);

%***************** AWGN addition ********* 

[ich5,qch5]=comb(ich44,qch44,attn);

%***************************  Receiver  *****************************
for (nuser_RO=1:NUser)
%****************** Perfect compensation *********

ifade2=1./ramp(nuser_RO,:).*ich5;
qfade2=1./ramp(nuser_RO,:).*qch5;

%****************** Guard interval removal *********

[ich6,qch6]= girem(ifade2,qfade2,fftlen2,gilen,nd);

%******************  FFT  ******************

rx_fft=complex(ich6,qch6);
ry_fft=fft(rx_fft);   	% fft : built in function
ich7=real(ry_fft);	% real : built in function
qch7=imag(ry_fft);	% imag : built in function

%***************** demoduration *******************

% ich8=ich7;%./ kmod;
% qch8=qch7;%./ kmod;
% combine i and q
Rx_para(:,1:2:nd.*ml)=ich7;
Rx_para(:,2:2:nd.*ml)=qch7;
% [demodata]=qpskdemod(ich8,qch8,para,nd,ml);   

%**************  Parallel to serial conversion  *****************

Rx(nuser_RO,:)=reshape(Rx_para,1,ChipLen);
end

%******************************* Receiver Serial( sigma ) *****************
% int i, j, m, nuser, it;

% Initialization
% int i, nuser;
% double s2 
s2_Ini = sigma .* sigma;

for( i_Ini=1:ChipLen )
	TotalMean(i_Ini) = 0;
	TotalVar(i_Ini) = s2_Ini;
end
for( nuser_Ini=1:NUser ) 
    for( i_Ini=1:ChipLen )
        Mean(nuser_Ini,i_Ini) = 0;
        Var(nuser_Ini,i_Ini) = H2(nuser_Ini);
        TotalVar(i_Ini) = TotalVar(i_Ini) + H2(nuser_Ini);			
    end
end

% % combine i and q
% Rx(1:2:end)=Rx_i;
% Rx(2:2:end)=Rx_q;

% end of Initialization

for( it_RS=1:ItNum ) 
    for( nuser_RS=1:NUser )
        % Produce the LLR values for the de-interleaved chip sequence.
        for( i_RS=1:ChipLen )	
            TotalMean(i_RS) = TotalMean(i_RS) - Mean(nuser_RS,i_RS);
            TotalVar(i_RS) = TotalVar(i_RS) - Var(nuser_RS,i_RS);
            Chip(nuser_RS,ScrambleRule(nuser_RS,i_RS)) = 2 .* H(nuser_RS) .* ( Rx(nuser_RS,i_RS) - TotalMean(i_RS) ) ./ TotalVar(i_RS);
        end

        % De-spreading operation.
        % int i, j, m;
        
        m_DeS = 1;
        for( i_DeS=1:DataLen )
            AppLlr(nuser_RS,i_DeS) = SpreadSeq(1) .* Chip(nuser_RS,m_DeS);
            m_DeS = m_DeS + 1;
            for( j_DeS=2:SpreadLen )
                AppLlr(nuser_RS,i_DeS) = AppLlr(nuser_RS,i_DeS) + SpreadSeq(j_DeS) .* Chip(nuser_RS,m_DeS);
                m_DeS = m_DeS + 1;
            end
        end
        
        % end of De-spreading

        % Feed the AppLlr to decoder, if there is a FEC codec in the system.

        % Spreading operation: Produce the extrinsic LLR for each chip
        m_RS = 1;
        for( i_RS=1:DataLen ) 
            for( j_RS=1:SpreadLen )
                Ext(nuser_RS,m_RS) = SpreadSeq(j_RS) .* AppLlr(nuser_RS,i_RS) - Chip(nuser_RS,m_RS);
                m_RS = m_RS + 1;
            end
        end

        % Update the statistic variable together with the interleaving process.
        for( i_RS=1:ChipLen )
	
            Mean(nuser_RS,i_RS) = H(nuser_RS) .* tanh( Ext(nuser_RS,ScrambleRule(nuser_RS,i_RS)) ./ 2 );
            Var(nuser_RS,i_RS) = H2(nuser_RS) - Mean(nuser_RS,i_RS) .* Mean(nuser_RS,i_RS);
            TotalMean(i_RS) = TotalMean(i_RS) + Mean(nuser_RS,i_RS);
            TotalVar(i_RS) = TotalVar(i_RS) + Var(nuser_RS,i_RS);
        end

    end
end

Data_RS = (AppLlr <= 0);
%******************************* Receiver Parallel( sigma ) ***************

%******************************* Error Counter ****************************
    error = error + sum(sum(abs(Data - Data_RS)));

    end
    
    Ber(i) = error ./ (NUser .* DataLen .* Block);
    Ber

end

figure(3)
semilogy([EbNoStart:EbNoStep:(EbNoNum-1).*EbNoStep], Ber, 's-b');
axis([EbNoStart (EbNoNum.*EbNoStep) min(Ber) 1]);
title(['OFDM\_IDMA QPSK RAY ',int2str(NUser),'USER  ',int2str(DataLen),'DATA ',int2str(SpreadLen),'SP  ',int2str(ItNum),'IT']);
xlabel('Eb/No(dB)');
ylabel('BER');
legend('BER');
grid on;
hold off;