afd_ofdm_system_j64qam2oldnew.m

好东西

clear all
close all
clc
%format long
%本次仿真载频为2GHz，带宽1MHz，子载波数64个，cp为8
%子载波间隔为15.625kHz
%一个ofdm符号长度为64us，cp长度为8us
%系统调制级别定为QPSK
%最大doppler频率为50Hz,信道是慢变的
%多径信道为4径，功率延迟谱服从负指数分布~exp(-t/trms),trms=(1/5)*tmax＝1.2us时长，各径延迟取为delay=[0 3e-6 6e-6]

bandwidth=2000000;%系统带宽为1MHz
num=4;
Pe=1e-6;
cp_length=16;%cp长度为16
groupnumber1=[ 16];
N_carrier=16;%OFDM子载波个数

 delay=[0 2e-6 3e-6 6e-6];
% delay=[0 1.5e-6 3e-6 4e-6 5e-6 6e-6];
%delay=[0 0.2e-6 0.5e-6 1.6e-6 2.3e-6 5e-6];
% delay=[0 5e-5 11e-5 17e-5 29e-5 31e-5];
% delay=[0 1e-6 2e-6 3e-6 4e-6 6e-6];
trms=2e-6;%
 var_pow=10*log10(exp(-delay/trms));%各径功率衰减,以dB形式给出
% averagepow=sum(var_pow)/num;
% A=[1 exp(-1) exp(-2) exp(-3) exp(-4) exp(-5)];
% var_pow=10*log10(A);
%var_pow=[-3 0 -2 -6 -8 -10 ];
% var_pow=[ 0 -3 -10 -18 -26 -32];
fd=50;


t_interval=(1/bandwidth)*N_carrier/(cp_length+N_carrier);%采样间隔64/72us，加上循环前缀后，采样率增加
delta_f=bandwidth/N_carrier;%Hz
SNR_dB=[5 10 15 20 25 30 35 ];%Eb/N0
% SNR_dB = 10.^(1:(1.2/10):1.7);
% SNR_dB=[20];
ber_snr_persjr=zeros(length(groupnumber1),length(SNR_dB));%最终给出各种SJR下误码率随SNR变化的曲线
% counter_begin=[1000];
counter_d=100000000;%各loop信道抽样跳动间隔
%counter=[1000,2000,3000];
% no=[6,7,6,7,6,7];
no=[10,10,10,10,10,10];
%   no=[20,20,20,20,20,20];
for ii=1:length(SNR_dB)%每个SNR点上仿真若干次
    
    snr=10^(SNR_dB(ii)/10);
    loop_num=100; %共仿真30次
    for kk=1:length(groupnumber1)
        error_bitold=0;
        error_bitnew=0;
        %错误比特数统计
        total_bit_numold=0;%发送总比特数统计
        total_bit_numnew=0;%发送总比特数统计
        groupsize=N_carrier/groupnumber1(kk);
        groupnumber=groupnumber1(kk);
        for l=1:loop_num
          
            N_ofdm=1;%每次仿真产生200个ofdm符号,则每次仿真共有200×64个星座映射符号；QPSK调制下，1个星座映射符号包含2个bit
%             counter_d=800*N_ofdm;
%             counter_begin=counter_begin+counter_d;
            counter_begin=(l+100000)*counter_d;


%                 counter_begin=l*counter_d;
            [Hk]=channel(N_carrier,cp_length,num,t_interval,var_pow,delay,fd,counter_begin,counter_d,no,N_ofdm);
%             Hk=abs(Hk).^2;
% %             plot
%             mesh(Hk);
%             figure;
%             contour(Hk);
            %调用比特功率分配函数得到各子载波组的调制方式和发送功率
%             [poweralloctpower,bitalloctvector]=bitandpoweralloct2(N_carrier,groupnumber,Hk,snr,Pe,N_ofdm);
               [poweralloctpowerold,bitalloctvectorold,poweralloctpowernew,bitalloctvectornew]=bitandpoweralloct2fenpei(N_carrier,groupnumber,Hk,snr,Pe,N_ofdm);
            powerold=sum(poweralloctpowerold);
            powernew=sum(poweralloctpowernew);
            %       if bitalloctvector==0
            %       continue
            %       end
            %       bitalloctvector;
            bit_maxnum_ofdmsigold=groupsize*max(sum(bitalloctvectorold));
            bit_maxnum_ofdmsignew=groupsize*max(sum(bitalloctvectornew));
            bit_sourceold=zeros(bit_maxnum_ofdmsigold,N_ofdm);
            bit_sourcenew=zeros(bit_maxnum_ofdmsignew,N_ofdm);
            bit_num_ofdmsigold=zeros(1,N_ofdm);
            bit_num_ofdmsignew=zeros(1,N_ofdm);
            map_outold=zeros(N_carrier,N_ofdm);
            map_outnew=zeros(N_carrier,N_ofdm);
            for nn=1:N_ofdm
            %nn=nn
                for v=1:1:groupnumber
                    %v=v
                    map_flagold=bitalloctvectorold(v,nn);
                    map_flagnew=bitalloctvectornew(v,nn);
                    %bitalloctvector(v)=bitalloctvector(v)
                    sourcebitold=zeros(1,bitalloctvectorold(v,nn)*groupsize);
                    sourcebitnew=zeros(1,bitalloctvectornew(v,nn)*groupsize);
                    for w=1:1:groupsize
                    % input=round(rand(1,bitalloctvector(v)));
                        if bitalloctvectorold(v,nn)>6
                            inputold=zeros(1,bitalloctvectorold(v,nn));
                        elseif bitalloctvectorold(v,nn)<1
                            inputold=[];
                        else
                            inputold=(rand(1,bitalloctvectorold(v,nn)))>0.5;
                        end
                        
                        
                        if bitalloctvectornew(v,nn)>6
                            inputnew=zeros(1,bitalloctvectorold(v,nn));
                        elseif bitalloctvectornew(v,nn)<1
                            inputnew=[];
                        else
                            inputnew=(rand(1,bitalloctvectornew(v,nn)))>0.5;
                        end
                        
                        
                        
                        
                        
                        
                        
                        
                        sourcebitold(1,1+(w-1)*bitalloctvectorold(v,nn):(w-1)*bitalloctvectorold(v,nn)+bitalloctvectorold(v,nn))=inputold;
                        sourcebitnew(1,1+(w-1)*bitalloctvectornew(v,nn):(w-1)*bitalloctvectornew(v,nn)+bitalloctvectornew(v,nn))=inputnew;
                        %按照map_flag指示完成各种星座映射,input为输入比特块
                    end
                    
                    bit_sourceold((v-1)*groupsize*map_flagold+1:(v-1)*groupsize*map_flagold+groupsize*map_flagold,nn)=sourcebitold';
                    bit_sourcenew((v-1)*groupsize*map_flagnew+1:(v-1)*groupsize*map_flagnew+groupsize*map_flagnew,nn)=sourcebitnew';

                    
                    
                    bit_num_ofdmsigold(nn)=bit_num_ofdmsigold(nn)+bitalloctvectorold(v,nn)*groupsize;
                    bit_num_ofdmsignew(nn)=bit_num_ofdmsignew(nn)+bitalloctvectornew(v,nn)*groupsize;
                    if length(sourcebitold)==0
                        map_outold((v-1)*groupsize+1:(v-1)*groupsize+groupsize,nn)=zeros(groupsize,1);
                    else
                    %map_out((v-1)*groupsize+1:(v-1)*groupsize+groupsize,nn)=sqrt(poweralloctpower(v,nn))*map_module(sourcebit',bitalloctvector(v,nn));
                    map_outold((v-1)*groupsize+1:(v-1)*groupsize+groupsize,nn)=map_module(sourcebitold',bitalloctvectorold(v,nn));
                    end
                    
                     if length(sourcebitnew)==0
                        map_outnew((v-1)*groupsize+1:(v-1)*groupsize+groupsize,nn)=zeros(groupsize,1);
                    else
                    %map_out((v-1)*groupsize+1:(v-1)*groupsize+groupsize,nn)=sqrt(poweralloctpower(v,nn))*map_module(sourcebit',bitalloctvector(v,nn));
                    map_outnew((v-1)*groupsize+1:(v-1)*groupsize+groupsize,nn)=map_module(sourcebitnew',bitalloctvectornew(v,nn));

                    end
                    
            end
     
        end
    
        %对一次仿真符号块进行映射,映射方式或者说调制方式由map_flag表征
        ofdm_modulation_outold=ifft(map_outold,N_carrier);%作64点逆FFT运算，完成ofdm调制,前面乘系数sqtr(64)是为了保持ifft前后的符号能量不变
        ofdm_modulation_outnew=ifft(map_outnew,N_carrier);
        ofdm_cp_outold=insert_cp(ofdm_modulation_outold,cp_length);%插入循环前缀 
        ofdm_cp_outnew=insert_cp(ofdm_modulation_outnew,cp_length);%插入循环前缀 

        % ofdm_cp_out=[ofdm_modulation_out((N_carrier-cp_length+1):N_carrier,:); ofdm_modulation_out];%插入循环前缀 
      
        [nnold,mmold]=size(ofdm_cp_outold);
        spow2=0;
        for k=1:nnold
            for b=1:mmold
                spow2=spow2+real(ofdm_cp_outold(k,b))^2+imag(ofdm_cp_outold(k,b))^2;
            end
        end
        spow1=spow2/(nnold*mmold);%信号平均能量
        map_flagold=sum(bit_num_ofdmsigold)/(N_carrier*N_ofdm);
        if map_flagold==0
            sgmaold=sqrt(spow1/(2*snr));
        else
            %sgma=sqrt(spow1/(2*snr)/map_flag);%sgma如何计算，与当前SNR和信号平均能量有关系
            sgmaold=sqrt(spow1/(2*snr));
        end
        
        [nnnew,mmnew]=size(ofdm_cp_outnew);
        spow2=0;
        for k=1:nnnew
            for b=1:mmnew
                spow2=spow2+real(ofdm_cp_outnew(k,b))^2+imag(ofdm_cp_outnew(k,b))^2;
            end
        end
        spow1=spow2/(nnnew*mmnew);%信号平均能量
        map_flagnew=sum(bit_num_ofdmsignew)/(N_carrier*N_ofdm);
        if map_flagnew==0
            sgmanew=sqrt(spow1/(2*snr));
        else
            %sgma=sqrt(spow1/(2*snr)/map_flag);%sgma如何计算，与当前SNR和信号平均能量有关系
            sgmanew=sqrt(spow1/(2*snr));
        end
        
        
        
    %********************** 以下过程为ofdm符号通过频率选择性多径信道 *************************
    %num=3;
    %假设功率延迟谱服从负指数分布~exp(-t/trms),trms=(1/4)*cp时长；
    %t在0~cp时长上均匀分布
    %若cp时长为16e-6s，可以取5径延迟如下
    %delay=[0 3e-6 6e-6];
    %trms=1.2e-6;
    %var_pow=10*log10(exp(-delay/trms));%各径功率衰减,以dB形式给出
    %fd=50;%最大doppler频率为50Hz
    %t_interval=0.8889e-6;%采样间隔64/72us，加上循环前缀后，采样率增加

    %counter_begin=(l-1)*counter_d;
%     trms_1=trms/t_interval;
%     t_max=16e-6/t_interval;
    %信道采样点数，每个调制符号采一个点


        [passchan_ofdm_symbolold,Hk]=multipath_chann(ofdm_cp_outold,num,var_pow,delay,fd,t_interval,counter_begin,counter_d,cp_length,no);
        [passchan_ofdm_symbolnew,Hk]=multipath_chann(ofdm_cp_outnew,num,var_pow,delay,fd,t_interval,counter_begin,counter_d,cp_length,no);

%         Hk=abs(Hk);
%         mesh(Hk);
    %[passchan_ofdm_symbol,Hk]=multipath_chann1(ofdm_cp_out,num,var_pow,delay,t_interval,cp_length);
%    counter=counter+counter_d;%更新各径信道取样开始位置
    %********************** 以上过程为ofdm符号通过频率选择性多径信道 *************************
    
     %********************** 以下过程为ofdm符号加高斯白噪声 *************************
%      passchan_ofdm_symbol=cut_cp(passchan_ofdm_symbol,cp_length);
%      Hk=passchan_ofdm_symbol./map_out;
        passnoise_ofdm_symbolold=add_noise(sgmaold,passchan_ofdm_symbolold);%加入随机高斯白噪声，receive_ofdm_symbol为最终接收机收到的ofdm符号块
        passnoise_ofdm_symbolnew=add_noise(sgmanew,passchan_ofdm_symbolnew);%加入随机高斯白噪声，receive_ofdm_symbol为最终接收机收到的ofdm符号块

    
    %********************** 以下过程为ofdm符号加BPJ干扰 *************************
%     jpow=spow1/sjr;%计算十进制信干比
        partial_fraction=0.25;
        jfre_begin=37;
     %receive_ofdm_symbol=add_pbj_jamming_new(passnoise_ofdm_symbol,jpow,delta_f,t_interval,partial_fraction,jfre_begin,cp_length);
    
    %********************** 以上过程为ofdm符号加高斯白噪声 *************************
        cutcp_ofdm_symbolold=cut_cp(passnoise_ofdm_symbolold,cp_length);%去除循环前缀
        cutcp_ofdm_symbolnew=cut_cp(passnoise_ofdm_symbolnew,cp_length);%去除循环前缀
        ofdm_demodulation_outold=fft(cutcp_ofdm_symbolold,N_carrier);%作128点FFT运算，完成ofdm解调
        ofdm_demodulation_outnew=fft(cutcp_ofdm_symbolnew,N_carrier);%作128点FFT运算，完成ofdm解调
        ofdm_demodulation_outold=ofdm_demodulation_outold./Hk;
        ofdm_demodulation_outnew=ofdm_demodulation_outnew./Hk;
         %ofdm_demodulation_out=ofdm_demodulation_out;
        receive_bit_sigold=de_map_module(ofdm_demodulation_outold,bitalloctvectorold,groupsize,groupnumber,bit_num_ofdmsigold);%解映射
        receive_bit_signew=de_map_module(ofdm_demodulation_outnew,bitalloctvectornew,groupsize,groupnumber,bit_num_ofdmsignew);%解映射

    %receive_bit_sig=de_map_module(ofdm_demodulation_out,map_flag);%解映射
 
    
    %以下过程统计接收信号中的错误比特数
        [mold,nold]=size(bit_sourceold);
        [mnew,nnew]=size(bit_sourcenew);
        %err_num=error_count(bit_source,receive_bit_sig);
        err_numold=sum(sum(rem(bit_sourceold+receive_bit_sigold,2)));
        err_numnew=sum(sum(rem(bit_sourcenew+receive_bit_signew,2)));
        error_bitold=error_bitold+err_numold;
        error_bitnew=error_bitnew+err_numnew;
        total_bit_numold=total_bit_numold+sum(bit_num_ofdmsigold);
        total_bit_numnew=total_bit_numnew+sum(bit_num_ofdmsignew);
    end
% %计算各种信噪比下的误比特率
%     error_bit;
%     total_bit_num;
    ber_snr_persjrold(kk,ii)=error_bitold/total_bit_numold;
    ber_snr_persjrnew(kk,ii)=error_bitnew/total_bit_numnew;
    fre_eff_persnrold(kk,ii)=total_bit_numold/(loop_num*N_ofdm*(1/bandwidth)*N_carrier)/bandwidth;
    fre_eff_persnrnew(kk,ii)=total_bit_numnew/(loop_num*N_ofdm*(1/bandwidth)*N_carrier)/bandwidth;
  end

end
 
 
SNR_dB_theo=0:0.1:35;%flat rayleigth fade
for i=1:length(SNR_dB_theo)
    SNR_theo=10.^(SNR_dB_theo(i)/10);
    ber_theo(1,i)=(1/2)*(1-sqrt(SNR_theo/(1+SNR_theo)));
end 

figure(1)
% semilogy(SNR_dB_theo,ber_theo(1,:),'b')
% hold on
semilogy(SNR_dB,ber_snr_persjrold(1,:),'r-*'),hold on
semilogy(SNR_dB,ber_snr_persjrnew(1,:),'g-*')
% semilogy(SNR_dB,ber_snr_persjr(2,:),'g-o')
% semilogy(SNR_dB,ber_snr_persjr(3,:),'m-+')
hold off
xlabel('SNR(dB)');
ylabel('仿真结果，实际Pe=1e-3');

figure(2)

plot(SNR_dB,fre_eff_persnrold(1,:),'r-*'),hold on
plot(SNR_dB,fre_eff_persnrnew(1,:),'g-*')
% plot(SNR_dB,fre_eff_persnr(1,:),'r-*',SNR_dB,fre_eff_persnr(2,:),'g-o',SNR_dB,fre_eff_persnr(3,:),'m-+')
xlabel('SNR(dB)');
ylabel('Spectral Efficiency(bps/Hz)');
% plot(SNR_dB,fre_eff_persnr(2,:),'g-o')
% plot(SNR_dB,fre_eff_persnr(3,:),'m-')
hold off