%%%新方案，接收端先進(jìn)行合并，將受到j(luò)amming影響的子載波去除

clear all
clc
format long
%本次仿真載頻為2GHz，帶寬1MHz，子載波數(shù)64個(gè)，cp為8
%子載波間隔為15.625kHz
%一個(gè)ofdm符號長度為64us，cp長度為8us
%系統(tǒng)調(diào)制級別定為QPSK
%最大doppler頻率為50Hz,信道是慢變的
%多徑信道為4徑，功率延遲譜服從負(fù)指數(shù)分布~exp(-t/trms),trms=(1/5)*tmax＝1.2us時(shí)長，各徑延遲取為delay=[0 3e-6 6e-6]

cp_length=16;%cp長度為16
N_carrier=64;%OFDM子載波個(gè)數(shù)
bandwidth=1000000;%系統(tǒng)帶寬為1MHz
t_interval=(1/bandwidth)*N_carrier/(cp_length+N_carrier);%采樣間隔64/72us，加上循環(huán)前綴后，采樣率增加
delta_f=bandwidth/N_carrier;%Hz
% SNR_dB=[0 4 8 12 16 20 24 28 32 35];%Eb/N0
% SJR_dB=[-10 -5 -3 0 3 5 10 15 20 25];%符號信干比
%SNR_dB=[8 12 16 20 24 28 32 35];
SNR_dB=[30];
Pe=1e-3;
groupnumber1=[4 8 16];
ber_snr_persjr_ofdm=zeros(length(groupnumber1),length(SNR_dB)); 
ber_snr_reallocatpow=zeros(length(groupnumber1),length(SNR_dB)); 
for kk=1:length(groupnumber1)
    groupsize=N_carrier/groupnumber1(kk);
    groupnumber=groupnumber1(kk);
    for i=1:length(SNR_dB)%每個(gè)SNR點(diǎn)上仿真若干次
    snr=10^(SNR_dB(i)/10);
    error_bit_ofdm=0;
    total_bit_num=0;%發(fā)送總比特?cái)?shù)統(tǒng)計(jì)
    error_bit_ofdm1=0;%for 二次功率分配
    total_bit_num1=0;%for 二次功率分配
    loop_num=10; %共仿真1000次

     for l=1:loop_num
    ofdm_symbol_num=12;%每次仿真產(chǎn)生10個(gè)ofdm符號,則每次仿真共有200×64個(gè)星座映射符號；QPSK調(diào)制下，1個(gè)星座映射符號包含2個(gè)bit
    %%%%信道參數(shù)%%%%%%
    num=6;

    delay=[0 1e-6 2e-6 3e-6 4e-6 6e-6];
    trms=1.5e-6;
    var_pow=10*log10(exp(-delay/trms));%各徑功率衰減,以dB形式給出
    fd=50;%最大doppler頻率為50Hz
    %t_interval=0.8889e-6;%采樣間隔64/72us，加上循環(huán)前綴后，采樣率增加
    counter=10000;%200000000;%各徑信道的采樣點(diǎn)間隔，應(yīng)該大于信道采樣點(diǎn)數(shù)。由以上條件現(xiàn)在信道采樣點(diǎn)數(shù)
    count_begin=(l-1+100000000)*5*counter;%每次仿真信道采樣的開始位置
    %%%%%%以上為信道參數(shù)%%%%%%
             %%%%以下過程為了得到準(zhǔn)確的信道信息%%%%%%
            map_out1=rand(2*N_carrier,ofdm_symbol_num)>0.5;
            map_out1=map_module(map_out1,2);
            ofdm_modulation_out1=sqrt(N_carrier)*ifft(map_out1,N_carrier);
            ofdm_cp_out1=insert_cp(ofdm_modulation_out1,cp_length);
            [passchan_ofdm_symbol1,Hk]=multipath_chann(ofdm_cp_out1,num,var_pow,delay,fd,t_interval,counter,count_begin,cp_length);
            sig_temp=cut_cp(passchan_ofdm_symbol1,cp_length);
            ofdm_demodulation_out1=fft(sig_temp,N_carrier)/sqrt(N_carrier);
            HHk=ofdm_demodulation_out1./map_out1;
            [nn,mm]=size(ofdm_cp_out1);
            spow2=0;
            for k=1:nn
                for b=1:mm
                    spow2=spow2+real(ofdm_cp_out1(k,b))^2+imag(ofdm_cp_out1(k,b))^2;
                end
            end
            %spow1=spow2/(nn*mm);%信號平均能量
            spow1=1;
            %%%%以上過程為了得到準(zhǔn)確的信道信息%%%%%%
            
            
            
            
    %%%%%%%%%以下為不做二次功率分配，原參考文獻(xiàn)方法%%%%%%%%%%%%
    [poweralloctpower,bitalloctvector,gama]=bitandpoweralloct2(N_carrier,groupnumber,HHk,snr,Pe,ofdm_symbol_num);%自適應(yīng)分配比特和功率
    
    bit_maxnum_ofdmsig=groupsize*max(sum(bitalloctvector));
    bit_source=zeros(bit_maxnum_ofdmsig,ofdm_symbol_num);
    bit_num_ofdmsig=zeros(1,ofdm_symbol_num);
    map_out=zeros(N_carrier,ofdm_symbol_num);
           
            for nn=1:ofdm_symbol_num
                for v=1:1:groupnumber
                    map_flag=bitalloctvector(v,nn);
                    sourcebit=zeros(1,bitalloctvector(v,nn)*groupsize);
                    for w=1:1:groupsize
                        if bitalloctvector(v,nn)>6
                            input=zeros(1,bitalloctvector(v,nn));
                        elseif bitalloctvector(v,nn)<1
                            input=[];
                        else
                            input=(rand(1,bitalloctvector(v,nn)))>0.5;
                        end
                        sourcebit(1,1+(w-1)*bitalloctvector(v,nn):(w-1)*bitalloctvector(v,nn)+bitalloctvector(v,nn))=input;
                        %按照map_flag指示完成各種星座映射,input為輸入比特塊
                    end
                    
                    bit_source((v-1)*groupsize*map_flag+1:(v-1)*groupsize*map_flag+groupsize*map_flag,nn)=sourcebit';
                    bit_num_ofdmsig(nn)=bit_num_ofdmsig(nn)+bitalloctvector(v,nn)*groupsize;
                    if length(sourcebit)==0
                        map_out((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_out((v-1)*groupsize+1:(v-1)*groupsize+groupsize,nn)=map_module_adp(sourcebit',bitalloctvector(v,nn));
                    end
                end
            end
            map_out1=zeros(size(map_out));
            %%%%%%%%以下根據(jù)功率分配值做分配功率%%%%%%%%%%
            for nn=1:ofdm_symbol_num
                for v=1:groupnumber
                    if poweralloctpower(v,nn)~=0
                       for w=1:groupsize
                         map_out1((v-1)*2+w,nn)=sqrt(poweralloctpower(v,nn))*map_out((v-1)*2+w,nn);    
                       end
                    end
                end
            end
    total_bit_num=total_bit_num+sum(bit_num_ofdmsig);


   ofdm_modulation_out=sqrt(N_carrier)*ifft(map_out1,N_carrier);%作64點(diǎn)逆FFT運(yùn)算，完成ofdm調(diào)制,前面乘系數(shù)sqtr(64)是為了保持ifft前后的符號能量不變
   
      ofdm_cp_out=insert_cp(ofdm_modulation_out,cp_length);%插入循環(huán)前綴 
      %%%%%%%%以下計(jì)算噪聲方差%%%%%%%%%%
            
            map_flag=sum(bit_num_ofdmsig)/(N_carrier*ofdm_symbol_num);
           % sgma=sqrt(spow1/(2*snr)/map_flag);
            if map_flag==0
                sgma=sqrt(spow1/(2*snr));
            else
            %sgma=sqrt(spow1/(2*snr)/map_flag);%sgma如何計(jì)算，與當(dāng)前SNR和信號平均能量有關(guān)系
            sgma=sqrt(spow1/(2*snr));
            end


    [passchan_ofdm_symbol,Hk]=multipath_chann(ofdm_cp_out,num,var_pow,delay,fd,t_interval,counter,count_begin,cp_length);
    
    passnoise_ofdm_symbol=add_noise(sgma,passchan_ofdm_symbol);%加入隨機(jī)高斯白噪聲，receive_ofdm_symbol為最終接收機(jī)收到的ofdm符號塊
  
     
    cutcp_ofdm_symbol=cut_cp(passnoise_ofdm_symbol,cp_length);%去除循環(huán)前綴
    
    ofdm_demodulation_out=fft(cutcp_ofdm_symbol,N_carrier)/sqrt(N_carrier);%作128點(diǎn)FFT運(yùn)算，完成ofdm解調(diào)
    %receive_ofdm_symbol=ofdm_demodulation_out;
    receive_ofdm_symbol=ofdm_demodulation_out./HHk;
    receive_ofdm_symbol1=zeros(size(receive_ofdm_symbol));
        for nn=1:ofdm_symbol_num
            for v=1:groupnumber
                if poweralloctpower(v,nn)~=0
                    for w=1:groupsize
                        receive_ofdm_symbol1((v-1)*2+w,nn)=receive_ofdm_symbol((v-1)*2+w,nn)/sqrt(poweralloctpower(v,nn));    
                    end
                end
            end
        end
  receive_bit_sig=de_map_module_adp(receive_ofdm_symbol1,bitalloctvector,groupsize,groupnumber,bit_num_ofdmsig);
 
    
    %以下過程統(tǒng)計(jì)接收信號中的錯(cuò)誤比特?cái)?shù)
           [mn,nn]=size(bit_source);
           err_num=sum(sum(rem(bit_source+receive_bit_sig,2)));
           error_bit_ofdm=error_bit_ofdm+err_num;
 %%%%%%%%%以上為不做二次功率分配，原參考文獻(xiàn)方法%%%%%%%%%%%%         
           
           
           
           
           
           
%%%%%%%%%%%以下為做功率二次分配的方法%%%%%%%%%%%%%%%%

    [poweralloctpower,bitalloctvector,gama]=reallocat_power(N_carrier,groupnumber,HHk,snr,Pe,ofdm_symbol_num);%做功率二次分配
    bit_maxnum_ofdmsig=groupsize*max(sum(bitalloctvector));
    bit_source=zeros(bit_maxnum_ofdmsig,ofdm_symbol_num);
    bit_num_ofdmsig=zeros(1,ofdm_symbol_num);
    map_out=zeros(N_carrier,ofdm_symbol_num);
           
            for nn=1:ofdm_symbol_num
                for v=1:1:groupnumber
                    map_flag=bitalloctvector(v,nn);
                    sourcebit=zeros(1,bitalloctvector(v,nn)*groupsize);
                    for w=1:1:groupsize
                        if bitalloctvector(v,nn)>6
                            input=zeros(1,bitalloctvector(v,nn));
                        elseif bitalloctvector(v,nn)<1
                            input=[];
                        else
                            input=(rand(1,bitalloctvector(v,nn)))>0.5;
                        end
                        sourcebit(1,1+(w-1)*bitalloctvector(v,nn):(w-1)*bitalloctvector(v,nn)+bitalloctvector(v,nn))=input;
                        %按照map_flag指示完成各種星座映射,input為輸入比特塊
                    end
                    
                    bit_source((v-1)*groupsize*map_flag+1:(v-1)*groupsize*map_flag+groupsize*map_flag,nn)=sourcebit';
                    bit_num_ofdmsig(nn)=bit_num_ofdmsig(nn)+bitalloctvector(v,nn)*groupsize;
                    if length(sourcebit)==0
                        map_out((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_out((v-1)*groupsize+1:(v-1)*groupsize+groupsize,nn)=map_module_adp(sourcebit',bitalloctvector(v,nn));
                    end
                end
            end
            map_out1=zeros(size(map_out));
            %%%%%%%%以下根據(jù)功率分配值做分配功率%%%%%%%%%%
            for nn=1:ofdm_symbol_num
                for v=1:groupnumber
                    if poweralloctpower(v,nn)~=0
                       for w=1:groupsize
                         map_out1((v-1)*2+w,nn)=sqrt(poweralloctpower(v,nn))*map_out((v-1)*2+w,nn);    
                       end
                    end
                end
            end
    total_bit_num1=total_bit_num1+sum(bit_num_ofdmsig);


   ofdm_modulation_out=sqrt(N_carrier)*ifft(map_out1,N_carrier);%作64點(diǎn)逆FFT運(yùn)算，完成ofdm調(diào)制,前面乘系數(shù)sqtr(64)是為了保持ifft前后的符號能量不變
   ofdm_cp_out=insert_cp(ofdm_modulation_out,cp_length);%插入循環(huán)前綴 
      %%%%%%%%以下計(jì)算噪聲方差%%%%%%%%%%
            
            map_flag=sum(bit_num_ofdmsig)/(N_carrier*ofdm_symbol_num);
           % sgma=sqrt(spow1/(2*snr)/map_flag);
            if map_flag==0
                sgma=sqrt(spow1/(2*snr));
            else
            %sgma=sqrt(spow1/(2*snr)/map_flag);%sgma如何計(jì)算，與當(dāng)前SNR和信號平均能量有關(guān)系
            sgma=sqrt(spow1/(2*snr));
            end
         


    [passchan_ofdm_symbol,Hk]=multipath_chann(ofdm_cp_out,num,var_pow,delay,fd,t_interval,counter,count_begin,cp_length);
    
    passnoise_ofdm_symbol=add_noise(sgma,passchan_ofdm_symbol);%加入隨機(jī)高斯白噪聲，receive_ofdm_symbol為最終接收機(jī)收到的ofdm符號塊
  
     
    cutcp_ofdm_symbol=cut_cp(passnoise_ofdm_symbol,cp_length);%去除循環(huán)前綴
    
    ofdm_demodulation_out=fft(cutcp_ofdm_symbol,N_carrier)/sqrt(N_carrier);%作128點(diǎn)FFT運(yùn)算，完成ofdm解調(diào)
    %receive_ofdm_symbol=ofdm_demodulation_out;
    receive_ofdm_symbol=ofdm_demodulation_out./HHk;
    receive_ofdm_symbol1=zeros(size(receive_ofdm_symbol));
        for nn=1:ofdm_symbol_num
            for v=1:groupnumber
                if poweralloctpower(v,nn)~=0
                    for w=1:groupsize
                        receive_ofdm_symbol1((v-1)*2+w,nn)=receive_ofdm_symbol((v-1)*2+w,nn)/sqrt(poweralloctpower(v,nn));    
                    end
                end
            end
        end
  receive_bit_sig=de_map_module_adp(receive_ofdm_symbol1,bitalloctvector,groupsize,groupnumber,bit_num_ofdmsig);

    
    %以下過程統(tǒng)計(jì)接收信號中的錯(cuò)誤比特?cái)?shù)
           [mn,nn]=size(bit_source);
           err_num=sum(sum(rem(bit_source+receive_bit_sig,2)));
           error_bit_ofdm1=error_bit_ofdm1+err_num;
%%%%%%%%%%%以上為做功率二次分配的方法%%%%%%%%%%%%%%%%                                  
end%for l=1:loop_num

ber_snr_persjr_ofdm(kk,i)=error_bit_ofdm/total_bit_num;
ber_snr_reallocatpow(kk,i)=error_bit_ofdm1/total_bit_num1;
fre_eff_persnr(kk,i)=total_bit_num/(loop_num*ofdm_symbol_num*(1/bandwidth)*N_carrier)/bandwidth;
fre_eff_persnr_reallocatpow(kk,i)=total_bit_num1/(loop_num*ofdm_symbol_num*(1/bandwidth)*N_carrier)/bandwidth;

save bitalloc_data ber_snr_persjr_ofdm ber_snr_reallocatpow fre_eff_persnr fre_eff_persnr_reallocatpow
end%for i=1:length(SNR_dB)

end %for kk
 
% SNR_dB_theo=0:0.1:30;%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
% % % for i=1:length(SNR_dB_theo)
% % % SNR_theo=10.^(SNR_dB_theo(i)/10);
% % % ber_theo(1,i)=Qfunct(sqrt(2*SNR_theo));
% % % end
% semilogy(SNR_dB_theo,ber_theo,'b')


%semilogy(SJR_dB,ber_snr_persjr_new(1,:),'b-o',SJR_dB,ber_snr_persjr_new(2,:),'b-+')



%  save precode_ts_new_data ber_snr_persjr_new ber_snr_persjr_ofdm