邮电通讯系统

开发平台：
Matlab

vblast_ofdm.m：源码内容
							clear all;
close all;
fprintf('**** Simulation starts at: ');
disp(datestr(now,0));
% pdp : Multipath channel power delay profile.
% 相邻径之间的相对时延为一个采样周期
% L 多径数目: 0至L-1径。
L=3;
var0=1;
%%%% var0=(1-exp(-1))/(1-exp(-L)); % 多径总功率归一化
ll=0:L-1;
pdp=var0.*exp(-ll);
%%%%%%%%%%%%%%%%%%%%%%%%%%
% 在每个子载波的范围内分配能量与比特
%%%%%%%%%%%%%%%%%%%%%%%%%%%
num=4;  % 平均每个天线的每个子载波分配的比特数
load parameters.txt;   % 装入仿真参数
cases=parameters(1);
caseparas=parameters(2);
for kk=1:cases
    % The number of transmit antennas.
    Nt=parameters((kk-1)*caseparas+3);
    % Nr: The number of receive antennas.
    Nr=parameters((kk-1)*caseparas+4);
    % N: The number of subcarriers.
    N=parameters((kk-1)*caseparas+5);
    % GI: CP length.
    GI=parameters((kk-1)*caseparas+6);
    % framenumber: The number of frames to simulate.
    framenumber=parameters((kk-1)*caseparas+7);
    % minSNR: Minimum SNR in dB to simulate begin with.
    minSNR=parameters((kk-1)*caseparas+8);
    % maxSNR: Maximum No in dB to simulate end with.
    maxSNR=parameters((kk-1)*caseparas+9);
    % stepSNR: SNR step in dB to increment.
    stepSNR=parameters((kk-1)*caseparas+10);
    % state: Modulation format: Adaptive=0; BPSK=1; QPSK=2; 16QAM=4; 64QAM=6.
    state=parameters((kk-1)*caseparas+11);
    % detection: Detection method: ZF=1; MMSE=2.
    detection=parameters((kk-1)*caseparas+12);
    
    
    
    % 线形星座预编码矩阵
    % LCPindication=0; % 为零的时候不使用LCP预编码，用串行干扰删除；
                     % 若为5表示不用串行干扰删除，且不用LCP预编码。
    % if LCPindication~=0
        % switch LCPindication
            % case 2
                % lcp=[exp(-j*pi/4) exp(-j*5*pi/4)].';
                % lcpint=[1 1].';
                % LCP=1/sqrt(2)*[lcpint lcpint.*lcp];
            % case 4
                % lcp=[exp(-j*pi/8) exp(-j*5*pi/8) exp(-j*9*pi/8) exp(-j*13*pi/8)].';
                % lcpint=[1 1 1 1].';
                % LCP=1/2*[lcpint lcpint.*lcp lcpint.*(lcp.^2) lcpint.*(lcp.^3)];
            % case 5
                % LCP=eye(Nt);
        % end
        % invLCP=inv(LCP);
    % end
 
    
    fprintf('==== Case %d ====n',kk);
    disp('---- Simulation parameter list:');
    fprintf('Transmit antennas: %d.n', Nt);
    fprintf('Receive antennas: %d.n', Nr);
    fprintf('The number of subcarriers : %d.n', N);
    fprintf('CP length : %d.n', GI);
    fprintf('Total frames: %d.n', framenumber);
    fprintf('SNR range to simulate: [%d:%1.1f:%d] dB.n', minSNR,stepSNR,maxSNR);
    
    switch(state)
        case 0
            disp('Adaptive modulation.');
            fmod='ADPM';
        case 1
            disp('BPSK modulation.');
            fmod='BPSK';
        case 2
            disp('QPSK modulation.');
            fmod='QPSK';       
        case 4
            disp('16QAM modulation.');
            fmod='QAM16';
        case 6
            disp('64QAM modulation.');
            fmod='QAM64';
    end        
    switch detection
        case 1
            disp('ZF detection.');
        case 2
            disp('MMSE detection.');
    end
    
    % fprintf('LCPindication: %d.n', LCPindication);
    
    disp('----------------------------');
    
    
    SNR=[minSNR:stepSNR:maxSNR];
    randn('state',0);
    rand('state',1);
    
    frame_length=N+GI; % 每帧的长度;
    err=zeros(1,length(SNR));
    averBER=zeros(1,length(SNR));
    Ptotmax=Nt*N;
    if state==0
        bitstotal=num*Nt*N*framenumber; % 每个发天线的每个子载波上平均分num个比特。对应比较2.^num阶调制。
    else
        bitstotal=Nt*N*framenumber*state; % 发送的总比特数
    end
     
    for SNRIndex=1:length(SNR) 
        
        sig2=Nt/(10.^(SNR(SNRIndex)/10));    % 此时的噪声功率为按整个ofdm符号平均能量归一化时的
        % 每个发送天线每次发送一帧
        last_symbol=zeros(Nt,length(pdp)-1); %信道记忆产生的影响
    
        for kk=1:framenumber 
            % 完成一个帧的发送与接收
            % 假设信道为准静态信道，在一个ofdm符号内保持不变
            % 在符号间隔处发生跃变
            [H H_f]=create_channel(Nt,Nr,pdp,N,GI); 
            % H_f为每个子载波对应的MIMO信道矩阵
            
          
           
            % if state==0
                % for ii=1:N
                %    Hnorm2(ii)=1/Nt/Nr*trace(H_f{ii}*H_f{ii}');
                % end
                % [Psubc,Rsubc]=BPallocRNew(Hnorm2/N/sig2,N,Ptotmax,N*num); % 每个发天线的每个子载波上平均分配的功率与比特数
            % end 
           
            w=cell(1,N);
            korder=cell(1,N);
            
            % G2=cell(1,N);    % 有LCP编码的时候用
            
       
            if detection==2
                eyenew=eye(Nr);  % Detection method: ZF=1; MMSE=2.
            end
            
            for f=1:N
                Hnew=H_f{f};
                % if LCPindication
                  
                  %  if detection==1  
                        % ZF
                   %     G2{f}=pinv(Hnew);    %pinv is the pseudoinverse
                       
                    % else   
                        % MMSE
                      %  G2{f}=Hnew'*pinv(Hnew*Hnew'+N*sig2*eyenew);
 
                    % end
                % else
                    
                korder{f}=zeros(1,Nt);
                krest=[1:Nt];
                for nt=1:Nt
                    if detection==1  
                            % ZF
                        G=pinv(Hnew);    %pinv is the pseudoinverse
                        if state==0
                            [k,postSNR]=argmaxzf(G,N*sig2,krest);  % find the maximum norm form krest columes of G
                        else
                            [k,postSNR]=argminzf(G,N*sig2,krest);  % find the minimum norm form krest columes of G
                        end
                    else   
                        % MMSE
                        G=Hnew'*pinv(Hnew*Hnew'+N*sig2*eyenew);
                        if state==0
                            [k,postSNR]=argmaxmmse(G,Hnew,N*sig2,krest);
                        else
                            [k,postSNR]=argminmmse(G,Hnew,N*sig2,krest);
                        end
                            eyenew(:,k)=zeros(Nr,1);
                    end
                        % 查找出最小的SNR所在的列,用零替代,产生新的信道伪逆
                    korder{f}(nt)=k;
                    w{f}(k,:)=G(k,:);
                    krest=setdiff(krest,k); % delete the selected k
                    Hnew(:,k)=zeros(Nr,1);
                    if state==0
                        % averpostSNR(SNRIndex,f)=averpostSNR(SNRIndex,f)+postSNR;
                        SNRsubch(k,f)=postSNR;
                    end
                end
            end
            % 产生发送数据
            
            
            tr_data=cell(Nt,N);
            symbol=zeros(Nt,N);
            
            if state==0
                
                SNRsubch2=reshape(SNRsubch.',1,Nt*N); % 变形后每个发天线的子载波信道对应信噪比占连续的一段。
                [P1,R1]=BPallocRNew(SNRsubch2,Nt*N,Ptotmax,num*Nt*N);
                P=reshape(P1,N,Nt).';    % 每个发天线的各子载波的能量占一行
                R=reshape(R1,N,Nt).';    % 每个发天线的各子载波的比特数占一行
                for nt=1:Nt
                    for f=1:N
                        tr_data{nt,f}=randn(1,R(nt,f))>0;
                          
                        switch R(nt,f)
                            case 0
                                symbol(nt,f)=zeros(1,1);
                            case 1           % 根据自适应算法比特增量为2，此调制方式没有被使用
                                symbol(nt,f)=BPSKMod(tr_data{nt,f});
                            case 2
                                symbol(nt,f)=QPSKMod(tr_data{nt,f});
                            case 4
                                symbol(nt,f)=QAM16Mod(tr_data{nt,f});
                            case 6
                                symbol(nt,f)=QAM64Mod(tr_data{nt,f});
                        end
                    end
                end
                
                Proot=sqrt(P);
                symbol=Proot.*symbol;
             
            else
                tr_data2=randn(Nt,N*state)>0;
                for nt=1:Nt
                    symbol(nt,:)=modulate(tr_data2(nt,:),state);
                end
            end
                
            % if LCPindication
                 % symbol=LCP*symbol;
            % end  
            
            % ofdm调制，进行ifft变换
          
            ofdm_symbol=ifft(symbol,[],2); 
           
            % 添加CP
            tr_symbol=[ofdm_symbol(:,end-GI+1:end) ofdm_symbol];
            
            % 发送信号平均功率归一化
            tr_symbol = tr_symbol*N/sqrt(N);
            
           
        
            % 通过信道
            re_cp_symbol= channel(H,tr_symbol,sig2,last_symbol,Nr,Nt,pdp);
        
            last_symbol=tr_symbol(:,end-length(pdp)+2:end); 
        
            % 恢复信号的功率水平
            re_cp_symbol=re_cp_symbol/sqrt(N);
        
            % 去除cp
            re_symbol=re_cp_symbol(:,GI+1:end);
    
            % ofdm解调，进行fft
            
            de_ofdm=fft(re_symbol,[],2);
         
            
   
            % 对每个载波进行vblast检测
            bitRx=cell(Nt,N);       % 存放自适应的接收比特
            y=zeros(Nt,N);
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            % if LCPindication
              %  for f=1:N
               %     y(:,f)=invLCP*G2{f}*de_ofdm(:,f);
                 %   for nt=1:Nt
                  %      if state==0
                   %         switch Rsubc(f)
                    %            case 0
                        
                     %           case 2
                      %              bitRx{nt,f}=QPSKDemod(y(nt,f));
                       %         case 4
                        %            bitRx{nt,f}=QAM16Demod(y(nt,f));
                         %       case 6
                          %          bitRx{nt,f}=QAM64Demod(y(nt,f));
                           % end
                        %else           
                         %   bitRx{nt,f}=demodulate(y(nt,f),state);
                        %end
                    %end
                %end
            %end
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
            for f=1:N  
                symbolRx=de_ofdm(:,f);
                for nt=1:Nt  
                    y=w{f}(korder{f}(nt),:)*symbolRx;
                    % decision and demodulation
                    if state==0
                        switch R(korder{f}(nt),f)
                            case 0
                                a=0;   
                            case 2
                                bitRx{korder{f}(nt),f}=QPSKDemod(y/Proot(korder{f}(nt),f));
                                a=QPSKMod(bitRx{korder{f}(nt),f})*Proot(korder{f}(nt),f); 
                            case 4
                                bitRx{korder{f}(nt),f}=QAM16Demod(y/Proot(korder{f}(nt),f));
                                a=QAM16Mod(bitRx{korder{f}(nt),f})*Proot(korder{f}(nt),f);
                            case 6
                                bitRx{korder{f}(nt),f}=QAM64Demod(y/Proot(korder{f}(nt),f));
                                a=QAM64Mod(bitRx{korder{f}(nt),f})*Proot(korder{f}(nt),f);
                        end
                    else           
                        bitRx{korder{f}(nt),f}=demodulate(y,state);
                        a=modulate(bitRx{korder{f}(nt),f},state);
                    end
                        symbolRx=symbolRx-H_f{f}(:,korder{f}(nt))*a; 
                end
            end
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            tr=[];
            if state==0
                for kk=1:Nt
                    for f=1:N
                        tr=[tr tr_data{kk,f}];
                    end
                end
            else
                tr=reshape(tr_data2.',1,Nt*N*state);
            end
            Rx=[];
            for kk=1:Nt
                for f=1:N
                    Rx=[Rx bitRx{kk,f}];
                end
            end
            err(SNRIndex)=err(SNRIndex)+sum(tr~=Rx); % 加上一帧中误比特数
        end   % 对应一帧检测完
        
    end  % 对应SNRIndex所有帧都检测完
    
    averBER=err/bitstotal;
    if state==0
        adap4BER11=averBER;
        save 'adap4BER11.mat' SNR adap4BER11; 
    else 
        qam16BER11=averBER;
        save 'qam16BER11.mat' SNR qam16BER11; 
    end
end  % End one case
% switch LCPindication
  %  case 2
   %     lcp22BER=averBER;
    %    save 'lcp22BER.mat' SNR averBER
    % case 4
      %  lcp44BER=averBER;
       % save 'lcp44BER.mat' SNR averBER
    % case 5
      %  nolcp_canBER44=averBER;
       % save 'nolcp_canBER44' SNR averBER
        
        
% end
figure;
semilogy(SNR,qam16BER11,'bo--',SNR,adap4BER11,'r^-');
axis([0 30 10.^-4 1]);
xlabel('average SNR (dB)');
ylabel('Raw BER');
legend('conventional','adptive');
grid on;
fprintf('n**** Ending time: ');
disp(datestr(now,0));
  
  
    
  
    
    
    
    
    
    
    

						