gmc_16qam_turbodd_pic.m

B3g_phase2_C语言_Matlab程序及说明

clear

% Define the slot structure
Spread_factor=8;                                                          % Spread factor
Subslot_number=8;                                                         % Number of subslots
Symbol_number_subslot=13;                                                 % Symbol number in a subslot
Symbol_number_slot=Subslot_number*Symbol_number_subslot;                  % Symbol number in a slot
Pilot_length=16;                                                          % Length of pilot signal in chips
Gaurd_length=8;                                                           % Length of gaurd time in chips
GaurdPilot_length=Gaurd_length+Pilot_length;                              % Total length of gaurd time and pilot signal in chips
SubslotData_length=Symbol_number_subslot*Spread_factor;                   % Length of data in a subslot in chips
Subslot_length=SubslotData_length+GaurdPilot_length;                      % Length of a subslot in chips
SlotData_length=Subslot_number*SubslotData_length;                        % Length of data in a slot in chips
Slot_length=Subslot_number*Subslot_length+GaurdPilot_length;              % Length of a slot in chips 

%Define the channel
Antenna_number=4;                                                         % Number of receive antennas
Path_number=8;                                                            % Number of paths  
% Path_Gain=[0.8084 0.462 0.253 0.259 0.0447 0.01];                        % Profile of channel model
% delay=[0 18 36  54 72  90]*2;                                            % Delays of paths

Delay_spread=10;                                                           % in us
Path_Gain=exp(-(0:1/(Path_number-1):1)*3*log(10));
Path_Gain=Path_Gain/sqrt(Path_Gain*Path_Gain');
delay=round((0:1/(Path_number-1):1)*Delay_spread*1.28);

Fc = 3.2e9;                                                               % Carrier frequency 
V =250;                                                                   % moving speed in km/h
Tc = 1/1.28e6/Sampling_factor;                                            % Chip width 
Time_Begin = 0;                                                           % Initializing the time
% Phase = 2*pi*rand(1,2);                                                  % Initializing the phase

% Define pilot signal
Pilot_sequence=[1 j 1 -1 1 1 -1 -j 1 -j 1 1 1 -1 -1 j];

%Pilot_sequence=[ 1  1  1 1  (1-j)/sqrt(2) -1  (1+j)/sqrt(2) -j...
%                 j -1 -j 1  (1+j)/sqrt(2)  1  (1-j)/sqrt(2) -j...
%                -1  1 -1 1 (-1+j)/sqrt(2) -1 (-1-j)/sqrt(2) -j...
%                -j -1  j 1 (-1-j)/sqrt(2)  1 (-1+j)/sqrt(2) -j];

for m=1:Path_number
    Pilot_matrix(m,:)=[Pilot_sequence(Pilot_length-m+2:Pilot_length) Pilot_sequence(1:Pilot_length-m+1)];
end
Gaurd_Pilot=[Pilot_sequence(Pilot_length-Gaurd_length+1:Pilot_length) Pilot_sequence];
Gaurd_Pilot1=[zeros(1,Gaurd_length) Pilot_sequence];

% turbo code
poly_g1 = 11;			                                                % 1 0 1 1
poly_g2 = 13;			                                                % 1 1 0 1
MemLen = 3;
puncture = 0;		                                                	% puncture = 0,rate=1/2;puncture = 1,no puncture rate =1/3
nIter = 7;
CodeRate = 1/(2+puncture);

%trel=poly2trellis(9,[561 753]); 
%tblen=40;

InputN=1; OutputN=2;
Code_Rate=InputN/OutputN;
Symbol_bitN=4;                                                           % 2 for QPSK and 4 for 16QAM

nIterDD=4;

%msg_L=Symbol_bitN*SlotData_length/OutputN*InputN*Slot_number-InputN*tblen;    % Number of message bits
Slot_number=64;
Packet_number=16;
Packet_Slot_number=Slot_number/Packet_number
Packet_msg_L=Symbol_bitN*Code_Rate*SlotData_length*Packet_Slot_number-MemLen;               
Packet_code_L=Symbol_bitN*SlotData_length*Packet_Slot_number;
Slot_code_L=Symbol_bitN*SlotData_length;

[Inner_intl_table] = interleave(Packet_msg_L+MemLen, floor(sqrt((Packet_msg_L+MemLen)/2)));
[Outer_intl_table] = interleave(Packet_code_L*Packet_number,floor(sqrt((Packet_code_L*Packet_number)/2)))+1;


% Bit error rate
ber=zeros(nIterDD,5);
format long

% Main loop

for Antenna_number=2.^(0:0);                                                    % Number of receive antennas
    
    SNR1=11; SNR2=12;
    
    Phase = 2*pi*rand(Antenna_number,2);                                        % Initializing the phase
    
    
    for SNR=SNR1:SNR2
        
        errors=zeros(1,nIterDD);
        
        for k=1:2000      
            SNR
            k
            % Transmiter
            % Transmiter: Generate transmitted bit stream
            msg=randint(Packet_number,Packet_msg_L,2);                          % Random data
            
            % Transmiter:Turbo encode and interleave
            
            for np=0:packet_number-1 
                code_p(np,:) = TuEnc(msg(np,:), puncture, Inner_intl_table, poly_g1, poly_g2);
            end
            code=reshape(code_p,1,Packet_code_L*Packet_number);
            code=code(Outer_intl_table);
            
            % Transmiter:Generate single-carrier baseband signals 
            
            s=[]; 
            for sn=0:Slot_number-1
                
                [x] = Modulation(code(sn*Slot_code_L+(1:Slot_code_L)), 4);
                
                Index_x=0;                                                                               
                for kk=1:Subslot_number
                    tmp=x(Index_x+(1:SubslotData_length));
                    tmp=[Gaurd_Pilot tmp]; 
                    s=[s tmp];
                    Index_x=Index_x+SubslotData_length;
                end
                s=[s Gaurd_Pilot];
            end
            
            % Channel
            Signal_length=Slot_number*Slot_length;                                 % Length of baseband signal 
            
            % Channel: The transmitted signal passes through multipath channel 
            for p=0:Path_number-1
                ss(p+1,:)=[zeros(1,p) s(1:Signal_length-p)];                            
            end
            
            for na=1:Antenna_number
                CH_Data = MultiCHannel(Path_Gain,Fc,V,Tc,Signal_length,Time_Begin,Phase(na,:));
                Sm(na,:)=sum(ss.*CH_Data);                                         
            end 
            
            Time_Begin = Time_Begin+Signal_length; 
            
            % Channel: Add Gauss noise
            Rm=awgn(Sm,SNR);%,'measured');
            
            Noise_variance=sum(sum(abs(Rm-Sm).^2))/length(Sm)/Antenna_number;
            
            %    Receiver
            
            for sn=0:Slot_number-1
                R=Rm(:,sn*Slot_length+(1:Slot_length));
                
                %    Receiver: Channel estimation 
                Index_R=Gaurd_length;Fades=zeros(Antenna_number,Path_number*(Subslot_number+1));Index_Fades=0;
                for kk=1:Subslot_number+1
                    Tmp_RP=R(:,Index_R+(1:Pilot_length));
                    Tmp_Fades=Tmp_RP*Pilot_matrix'/Pilot_length;
                    Fades(:,Index_Fades+(1:Path_number))=Tmp_Fades;
                    Index_R=Index_R+Subslot_length;
                    Index_Fades=Index_Fades+Path_number;
                end
                
                for p=1:Path_number
                    tmp_Fades=Fades(:,p:Path_number:end);
                    for na=1:Antenna_number
                        [Coefficients,Structure]=polyfit(0:Subslot_number,tmp_Fades(na,:),3);
                        [tmp_FadesP(na,:),delta]=polyval(Coefficients,0:0.5:Subslot_number,Structure);
                    end
                    Fades(:,p:Path_number:length(Fades))=tmp_FadesP(:,1:2:end);
                    FadesI(:,p:Path_number:length(Fades)-Path_number)=tmp_FadesP(:,2:2:end);
                end
                
                Index_R=Gaurd_length;Index_Fades=0;
                for kk=1:Subslot_number+1
                    Tmp_Fades=Fades(:,Index_Fades+(1:Path_number));
                    Tmp_Noise = R(:,Index_R+(1:Pilot_length))-Tmp_Fades(:,1:Path_number)*Pilot_matrix;
                    Nv(kk)=sum(sum(abs(Tmp_Noise.*Tmp_Noise)))/Pilot_length/Antenna_number;                                    % Estimate of noise variance 
                    Index_R=Index_R+Subslot_length;
                    Index_Fades=Index_Fades+Path_number;
                end
                Tmp_Nv=sum(Nv)/(Subslot_number+1);
                %    Receiver: Equalization in DFT domain
                Index_R=Pilot_length; Index_Fades=0;R_EQ=[];
                for kk=0:Subslot_number-1
                    Tmp_Fades=(FadesI(:,Index_Fades+(1:Path_number)));
                    RM=R(:,Index_R+(1:Subslot_length)); 
                    
                    RMC((sn*Subslot_number+kk)*Subslot_length+(1:Subslot_length))=sum(cyc_cov(RM,Tmp_Fades),1);
                    
                    h_half((sn*Subslot_number+kk)*Path_number+(1:Path_number))=sum(self_cov(Tmp_Fades),1);
                    
                    Index_Fades=Index_Fades+Path_number;
                    Index_R=Index_R+Subslot_length;
                end
            end   
            
            for n=1:nIterDD            
                for sn=0:Slot_number-1
                    
                    for kk=0:Subslot_number-1
                        
                        Tmp_RMC=RMC((sn*Subslot_number+kk)*Subslot_length+(1:Subslot_length));
                        
                        Tmp_h_half=h_half((sn*Subslot_number+kk)*Path_number+(1:Path_number));                        
                        Tmp_h=[conj(Tmp_h_half(end:-1:2)) 0 Tmp_h_half(2:end)];
                        
                        if n==1
                            
                            x_rake=Tmp_RMC(Gaurd_length+(1:SubslotData_length));
                            Tmp_R_EQ=x_rake;
                            rou=Tmp_h_half(1);
                            x_sigma=ones(1,SubslotData_length)*(Noise_variance*rou+abs(Tmp_h*Tmp_h'));
                            for nn=1:3
                                [x_symbol,x_variance] =Symbol_Decision(Tmp_R_EQ,rou,x_sigma);
                                
                                x_mean=[Gaurd_Pilot(Pilot_length+(Gaurd_length-Path_number+2:Gaurd_length)) x_symbol Gaurd_Pilot(1:Pilot_length)];
                                x_variance=[zeros(1,Path_number-1) x_variance zeros(1,Pilot_length)];
                                
                                [x_intf, x_sigma]=Intf_Sigma(x_mean,x_variance,Tmp_h,SubslotData_length);
                                
                                Tmp_R_EQ=x_rake-x_intf;
                                
                                x_sigma=Noise_variance*rou+x_sigma;
                            end
                            LLR_D = Soft_Demod(Tmp_R_EQ, rou,x_sigma ,zeros(4,SubslotData_length), SubslotData_length);
                            
                        else
                            LLR_P=reshape(LLR_DI(nc,(sn*Subslot_number+kk)*SubslotData_length*4+(1:SubslotData_length*4)),4,SubslotData_length); 
                            
                            [x_mean0, x_variance]=Mean_Var(LLR_P);
                            
                            x_mean=[Gaurd_Pilot(Pilot_length+(Gaurd_length-Path_number+2:Gaurd_length)) x_mean0 Gaurd_Pilot(1:Pilot_length)];
                            x_variance=[zeros(1,Path_number-1) x_variance zeros(1,Pilot_length)];
                            rou =Tmp_h_half(1);
                            
                            [x_intf, x_sigma]=Intf_Sigma(x_mean,x_variance,Tmp_h,SubslotData_length);
                            
                            x_rake=Tmp_RMC(Gaurd_length+(1:SubslotData_length));
                            x_rec=x_intf+rou*x_mean0;
                            alpha=abs((x_rake*x_rec')/(x_rake*x_rake'));
                            
                            Tmp_R_EQ=x_rake-alpha*x_intf;
                            
                            x_sigma=Noise_variance*rou+alpha^2*x_sigma;
                            
                            LLR_D = Soft_Demod(Tmp_R_EQ, rou,x_sigma ,LLR_P, SubslotData_length);
                            
                        end
                        
                        [Dem_signal((sn*Subslot_number+kk)*SubslotData_length*4+(1:SubslotData_length*4))]=LLR_D;
                    end
                    
                    Dem_signal(Outer_intl_table)=Dem_signal;
                    Dem_signal_p=vec2mat(Dem_signal,Packet_number)';
                    
                    for np=1:Packet_number
                        
                        [decoded, LLR_all(np,:)] = TuDecLogMapNew(Dem_signal_p(np,:), puncture, nIter, Inner_intl_table, 1, 1, poly_g1, poly_g2);
                        errors(n)=errors(n)+sum(abs(decoded(1:msg_L)-msg(np,1:msg_L)));
                        
                    end
                    
                    LLR_DI=reshape(LLR_all,1,Packet_code_L*Packet_number);
                    LLR_DI=LLR_DI(Outer_intl_table);
                    
                end 
            end
            
            errors
            ber(:,SNR-SNR1+1)=errors'/k/msg_L/packet_N
            
            
            if (errors(nIterDD)>300 & k>3) 
                break;
            end
        end
        
        if ber(1,SNR-SNR1+1)<1.0*10^(-6)
            break;
        end
        
    end 
    save ber_16QAM200v_1 ber
end
semilogy(SNR1:SNR2,ber(:,1:SNR2-SNR1+1)')
grid
xlabel('SNR of Received Signal(in dB)')
ylabel('Bit Error Rate')
pause(0.2)
save ber_16QAM200v_1 ber