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