% BEGIN, main function
function channel_estimation_main(SNRs, pilot_time_spacing, pilot_frequency_spacing, speed_v_kmh, B_decay, window_t, window_f, method_sel, mode ,webdemo_en, filename)
% BEGIN, set up graphical output

% Create an invisible figure.
if (webdemo_en)
    fig = figure(1); set(fig, 'visible', 'off');
    set(fig,'PaperUnits','inches');
    set(fig,'PaperPosition',[0 0 8.889 6.667]); % actual size is determined by 8.889, 6.667 x 72 dpi = 640 x 480
    set(fig,'PaperPositionMode','manual');
else
    fig = figure(1); set(fig, 'visible', 'on');
end
% set default Fontsize and name

set(gca, 'defaultlinelinewidth', 4);
set(gca, 'defaulttextfontname', 'Times New Roman');
set(gca, 'fontsize', 23);

rng('shuffle'); % true random signal and channel coefficients

    fftSize = 256;
    nPRB = 12; % Physical Resource Block
    NsymbsperTTI = 201; % number of total symbols on time domain
    %SNRs = 20; % signal to Noise ratio in dB
    Ncp = 17; % cyclic prefix length
    nSubband = 6; % number of subbands
    Totactivesub = nSubband*nPRB;  %total numberof active sub carriers
    CPOFDM_length = fftSize+Ncp;
    bitsPerSymbol = 2; % M = 4 Size of constellation M, number of bits for one symbol
    GrayLabeling = 1;
    mapPamQam = 1; % QAM
    [symlabel,QAM_Symbols] = get_mapping_and_labeling(bitsPerSymbol,GrayLabeling,mapPamQam);
    deltaF = 15e3;                        % Frequency spacing between adjacent sub carriers
    bandwidth = fftSize*deltaF;           % bandwidth to plot, in Hz
    
    %% multipath channel models & %EDPDP parameters
    %B_decay = 0.6; % decay factor of EDPDP channel
    ch_Length = Ncp+1; % the delay spread length
    
    %% channel estimation pilot symbols
    %pilot_time_spacing = 10; % the interval between two pilots on time domain
    %pilot_frequency_spacing = 4; % the interval between two pilots on frequency domain
    
    %% doppler effect parameters
    %speed_v_kmh = 100;                    % Doppler shift defined through vehicle speed: in km/h
    fc = 2.0e9;                           % carrier frequency in Hz
    
    enableBF = 1;                         % with Block Fading or not
    
    fD = fc*(speed_v_kmh*1000/3600.0)/3e8;  % calculate the maximum Doppler shift fD
    T = 1/bandwidth;     % symbol duration for a modulated carrier.
    fD_T = fD*T;
    fD_T_LMMSE = fD_T*CPOFDM_length;
    
    %% ofdm system simulation
    
    pilots_nb_time = ceil(NsymbsperTTI/pilot_time_spacing); % number of pilots on time domain
    pilots_nb_fre = ceil(Totactivesub/pilot_frequency_spacing); % number of pilots on frequency domain
    
    NsymbsperTTI = pilots_nb_time*pilot_time_spacing - pilot_time_spacing + 1; % ensure last symbol to be pilot symbol on time-domain
    Totactivesub_temp = pilots_nb_fre*pilot_frequency_spacing - pilot_frequency_spacing + 1; % ensure last subcarrier to be pilot subcarrier on fre-domain
    
    pilots_time_instances = 1:pilot_time_spacing:NsymbsperTTI;           % pilot positions in the time domain
    pilots_frequency_instances = 1:pilot_frequency_spacing:Totactivesub; % pilot positions in the frequency domain
    
    R_Pilot_t = 1/((pilot_time_spacing+1) * T * CPOFDM_length); % pilot rate in the time domain
     
    N_sym_OFDM = CPOFDM_length * NsymbsperTTI; % total number of OFDM signals
    
    %% correlated sequences generation
    fadingDrop = CPOFDM_length;  % no. of fading sample drops
    t_discrete_r = (2*pi*fD_T).*(0:fadingDrop:N_sym_OFDM-1);
    J_1 = besselj(0,t_discrete_r); % Rayleich fading process model ,bessel func with first kind and order 0
    R_design = toeplitz(J_1);
    R_chol = cholcov(R_design); % cholesky decomposition
    [NoDoF,~] = size(R_chol);
    
    s_received = zeros(1,N_sym_OFDM); % the signal should be cleared to 0
    
    %%  data and pilot symbol generation
    [ s_mapped,pilot_symbols] = pilot_data_generation( nPRB,nSubband,Totactivesub,NsymbsperTTI, bitsPerSymbol,symlabel,QAM_Symbols, pilots_frequency_instances, pilots_time_instances );
    
    %% ifft processing
    x_ofdm = ifft(s_mapped,fftSize)*sqrt(fftSize);
    
    %% Add cyclic prefix of length Ncp
    x_ofdm_vec = ofdm_add_cp( x_ofdm,Ncp,fftSize,NsymbsperTTI );
    
    %% channel gerneration
    [ h_channel,tap_delay_f,h_weights_pdp] = EDPDP_channel_model( ch_Length,B_decay,N_sym_OFDM,R_chol,fadingDrop);   
    h_channel_pdp = [h_weights_pdp;zeros(fftSize-length(h_weights_pdp),1)];
    
    %% generate the AWGN noise
    sqrt_nvar = 1/sqrt(10^(0.1*SNRs));
    noise_power = sqrt_nvar^2;
    n = 1/sqrt(2)*sqrt_nvar*(randn(1,N_sym_OFDM) + 1i*randn(1,N_sym_OFDM));
    
    %% signal go through the channel
    for i = 1:length(tap_delay_f)
        h_channel_r = h_channel(i,(1:end-tap_delay_f(i)+1));
        h_channel_rr = [zeros(1,tap_delay_f(i)-1),h_channel_r.*(x_ofdm_vec(1:end-tap_delay_f(i)+1).')];
        s_received = s_received + h_channel_rr;
    end
    s_received = s_received + n; % received signal with noise
    
    %% cp removal and fft processing
    y_ofdm = reshape(s_received,[CPOFDM_length,NsymbsperTTI]);
    y_ofdm_nocp = y_ofdm(Ncp+1:end,:);
    y_ofdm_fft = fft(y_ofdm_nocp,fftSize)/sqrt(fftSize);
    
    % h_channel estimation by zero-forcing and mmse method
    [ ~,h_channel_est_zf ] = channel_get_estimation( y_ofdm_fft, pilots_time_instances, pilots_frequency_instances, pilot_symbols, noise_power  );
    
    %% original channel calculation
    h_channel_BF = h_channel(:,1:fadingDrop:N_sym_OFDM);
    h_channel_current = zeros(fftSize,NsymbsperTTI);
    for i = 1:length(tap_delay_f)
        h_channel_current(tap_delay_f(i),:) = h_channel_BF(i,:);
    end
    h_channel_current_fft = fft(h_channel_current,fftSize);
    h_channel_fft_diag = h_channel_current_fft(1:Totactivesub_temp,:);
    
    %% interpolation and MSE calculation
    % sliding window parameters
    enable_SW = 1; % with sliding window or not
    
    %window_t = 6; % sliding window length on time domain
    %window_f = 6; % sliding window length on frequency domain
    
    if(enable_SW == 1) % degree selection for LS legendry polynomials
        LE_degree_t = floor(window_t/2)+1;
        LE_degree_f = floor(window_f/2);
    else
        LE_degree_t = 9;
        LE_degree_f = 9;
    end
    
    DFT_Degree_t = min(NoDoF,pilots_nb_time); % degree selection for DFT polynomials
    DFT_Degree_f = length(tap_delay_f);
    
    t_axis = 1:1:NsymbsperTTI;
    t_axis_pilots = pilots_time_instances;
    
    f_axis = 1:1:Totactivesub_temp;
    f_axis_pilots = pilots_frequency_instances;
    
    % linear interpolation
    if (method_sel(1))
        est_channel_li_time = zeros(pilots_nb_fre,NsymbsperTTI);
        est_channel_li = zeros(Totactivesub_temp,NsymbsperTTI);
        for i = 1:pilots_nb_fre
            est_channel_li_time(i,:) = interp1(pilots_time_instances, h_channel_est_zf(i,:),1:1:NsymbsperTTI,'linear');
        end
        
        for j = 1:NsymbsperTTI
            est_channel_li(:,j) = interp1(pilots_frequency_instances,est_channel_li_time(:,j),1:1:Totactivesub_temp,'linear');
        end
        MSE_li_t = mean((abs(h_channel_fft_diag(1,:) - est_channel_li(1,:))).^2);
        MSE_li_f = mean((abs(h_channel_fft_diag(:,1) - est_channel_li(:,1))).^2);
        li_dist_mat = 10*log10(abs(h_channel_fft_diag(:,:) - est_channel_li(:,:)).^2);
        MSE_li = mean(mean(abs(h_channel_fft_diag(:,:) - est_channel_li(:,:)).^2));
        
        MSE_li_db = 10*log10(MSE_li);
        MSE_li_t_db = 10*log10(MSE_li_t);
        MSE_li_f_db = 10*log10(MSE_li_f);
        
        if (mode)
            clims = [-30 -5];
            imagesc(1:1:NsymbsperTTI,1:1:Totactivesub_temp,li_dist_mat,clims)
            colorbar('FontSize',14)
            title(colorbar,'in dB','FontSize',14)
            xlabel('Time Lag')
            ylabel('Subcarrier Index')
            title(['MSE: ',num2str(MSE_li_db,4),'dB',', MAX Doppler f_D = ',num2str(fD,4),'Hz',', Pilot Rate R_P = ',num2str(R_Pilot_t,5),'Hz'])
            
            %set(gca,'units','pixels','outerposition',[-15,-5,640,430])
            print(filename,'-dpng','-r72');
        else
            subplot(2,1,1)
            plot(t_axis,real(h_channel_fft_diag(1,:)),'r-','linewidth',1.5)
            hold on
            plot(t_axis,real(est_channel_li(1,:)),'b-.','linewidth',1.5);
            plot(t_axis_pilots,real(h_channel_est_zf(1,:)),'m*');
             xlabel('Time Lag','FontSize',14)
            ylabel('real part of H','FontSize',14)
            xlim([0 205]);
            legend({'Ori','LI','Pilot'},'FontSize',14,'Location','southeast')
            title(['Time Domain, MSE: ',num2str(MSE_li_t_db,4),'dB'],'FontSize',15)
            
            subplot(2,1,2)
            plot(f_axis,real(h_channel_fft_diag(:,1)),'r-','linewidth',1.5)
            hold on
            plot(f_axis,real(est_channel_li(:,1)),'b-.','linewidth',1.5);
            plot(f_axis_pilots,real(h_channel_est_zf(:,1)),'m*');
             xlabel('Subcarrier Index','FontSize',14)
            ylabel('real part of H','FontSize',14)
            legend({'Ori','LI','Pilot'},'FontSize',14,'Location','southeast')
            title(['Frequency Domain, MSE: ',num2str(MSE_li_f_db,3),'dB'],'FontSize',15)
            
            %set(fig, 'Position', [0 0 640 480])
            %set(gca,'units','pixels','outerposition',[-15,-5,640,430])
            print(filename,'-dpng','r72');
        end
        
        % Least square with Legendre polynomials
    elseif (method_sel(2))
        est_channel_ls_le_time = zeros(pilots_nb_fre,NsymbsperTTI);
        est_channel_ls_le = zeros(Totactivesub_temp,NsymbsperTTI);
        
        for i = 1:pilots_nb_fre
            est_channel_ls_le_time(i,:) = ls_time_legendry( pilots_time_instances, NsymbsperTTI, h_channel_est_zf(i,:),pilot_time_spacing, LE_degree_t,enable_SW, window_t);
        end
        
        for j = 1:NsymbsperTTI
            est_channel_ls_le(:,j) = ls_fre_legendry( pilots_frequency_instances, Totactivesub_temp, est_channel_ls_le_time(:,j),pilot_frequency_spacing,LE_degree_f,enable_SW, window_f);
        end
        
        MSE_ls_le_t = mean((abs(h_channel_fft_diag(1,:) - est_channel_ls_le(1,:))).^2);
        MSE_ls_le_f = mean((abs(h_channel_fft_diag(:,1) - est_channel_ls_le(:,1))).^2);
        ls_le_dist_mat = 10*log10(abs(h_channel_fft_diag(:,:) - est_channel_ls_le(:,:)).^2);
        MSE_ls_le = mean(mean(abs(h_channel_fft_diag(:,:) - est_channel_ls_le(:,:)).^2));
        
        MSE_ls_le_t_db = 10*log10(MSE_ls_le_t);
        MSE_ls_le_f_db = 10*log10(MSE_ls_le_f);
        MSE_ls_le_db = 10*log10(MSE_ls_le);
                
        if (mode)
            clims = [-30 -5];
            imagesc(1:1:NsymbsperTTI,1:1:Totactivesub_temp,ls_le_dist_mat,clims)
            colorbar('FontSize',14)
            title(colorbar,'in dB','FontSize',14)           
            xlabel('Time Lag')
            ylabel('Subcarrier Index')
            title(['MSE: ',num2str(MSE_ls_le_db,4),'dB',', MAX Doppler f_D = ',num2str(fD,4),'Hz',', Pilot Rate R_P = ',num2str(R_Pilot_t,5),'Hz'])
            
            %set(gca,'units','pixels','outerposition',[-15,-5,640,430])
            print(filename,'-dpng','-r72');         
        else
            subplot(2,1,1)
            plot(t_axis,real(h_channel_fft_diag(1,:)),'r-','linewidth',1.5)
            hold on
            plot(t_axis,real(est_channel_ls_le(1,:)),'b-.','linewidth',1.5);
            plot(t_axis_pilots,real(h_channel_est_zf(1,:)),'m*');
            xlabel('Time Lag','FontSize',14)
            ylabel('real part of H','FontSize',14)
            xlim([0 205]);
            legend({'Ori','LS LE','Pilot'},'FontSize',14,'Location','southeast')
            title(['Time Domain, MSE: ',num2str(MSE_ls_le_t_db,4),'dB'],'FontSize',15)
            
            subplot(2,1,2)
            plot(f_axis,real(h_channel_fft_diag(:,1)),'r-','linewidth',1.5)
            hold on
            plot(f_axis,real(est_channel_ls_le(:,1)),'b-.','linewidth',1.5);
            plot(f_axis_pilots,real(h_channel_est_zf(:,1)),'m*');
            xlabel('Subcarrier Index','FontSize',14)
            ylabel('real part of H','FontSize',14)
            legend({'Ori','LS LE','Pilot'},'FontSize',14,'Location','southeast')
            title(['Frequency Domain, MSE: ',num2str(MSE_ls_le_f_db,4),'dB'],'FontSize',15)
            
            %set(fig, 'Position', [-10 -10 640 480])
            print(filename,'-dpng','r72');
        end
        
        % Least squares with DFT polynomials
    elseif (method_sel(3))
        est_channel_ls_dft_time = zeros(pilots_nb_fre,NsymbsperTTI);
        est_channel_ls_dft = zeros(Totactivesub_temp,NsymbsperTTI);
        
        for i = 1:pilots_nb_fre
            est_channel_ls_dft_time(i,:) = ls_time_dft(pilots_time_instances,h_channel_est_zf(i,:),R_chol,DFT_Degree_t,enable_SW,NsymbsperTTI,window_t,pilot_time_spacing,CPOFDM_length,fadingDrop,enableBF);
        end
        
        for j = 1:NsymbsperTTI
            est_channel_ls_dft(:,j) = ls_fre_dft( pilots_frequency_instances,est_channel_ls_dft_time(:,j),Totactivesub_temp,fftSize,tap_delay_f,DFT_Degree_f,window_f,pilot_frequency_spacing,enable_SW );
        end
        
        MSE_ls_dft_t = mean((abs(h_channel_fft_diag(1,:) - est_channel_ls_dft(1,:))).^2);
        MSE_ls_dft_f = mean((abs(h_channel_fft_diag(:,1) - est_channel_ls_dft(:,1))).^2);

        ls_dft_dist_mat = 10*log10(abs(h_channel_fft_diag(:,:) - est_channel_ls_dft(:,:)).^2);
        MSE_ls_dft = mean(mean(abs(h_channel_fft_diag(:,:) - est_channel_ls_dft(:,:)).^2));
        
        MSE_ls_dft_t_db = 10*log10(MSE_ls_dft_t);
        MSE_ls_dft_f_db = 10*log10(MSE_ls_dft_f);
        MSE_ls_dft_db = 10*log10(MSE_ls_dft);
        
        if (mode)
            clims = [-30 -5];
            imagesc(1:1:NsymbsperTTI,1:1:Totactivesub_temp,ls_dft_dist_mat,clims)
            colorbar('FontSize',12)
            title(colorbar,'in dB','FontSize',14) 
            xlabel('Time Lag')
            ylabel('Subcarrier Index')
            title(['MSE: ',num2str(MSE_ls_dft_db,4),'dB',', MAX Doppler f_D = ',num2str(fD,4),'Hz',', Pilot Rate R_P = ',num2str(R_Pilot_t,5),'Hz'])
            
            set(gca,'units','pixels','outerposition',[-15,-5,640,430])
            print(filename,'-dpng','-r72');
        else
            subplot(2,1,1)
            plot(t_axis,real(h_channel_fft_diag(1,:)),'r-','linewidth',1.5)
            hold on
            plot(t_axis,real(est_channel_ls_dft(1,:)),'b-.','linewidth',1.5);
            plot(t_axis_pilots,real(h_channel_est_zf(1,:)),'m*');
            xlabel('Time Lag','FontSize',14)
            ylabel('real part of H','FontSize',14)
            xlim([0 205]);
            legend({'Ori','Eig','Pilot'},'FontSize',14,'Location','southeast')
            title(['Time Domain, MSE: ',num2str(MSE_ls_dft_t_db,4),'dB'],'FontSize',15)
            
            subplot(2,1,2)
            plot(f_axis,real(h_channel_fft_diag(:,1)),'r-','linewidth',1.5)
            hold on
            plot(f_axis,real(est_channel_ls_dft(:,1)),'b-.','linewidth',1.5);
            plot(f_axis_pilots,real(h_channel_est_zf(:,1)),'m*');
            xlabel('Subcarrier Index','FontSize',14)
            ylabel('real part of H','FontSize',14)
            legend({'Ori','Eig','Pilot'},'FontSize',14,'Location','southeast')
            title(['Frequency Domain, MSE: ',num2str(MSE_ls_dft_f_db,4),'dB'],'FontSize',15)
            
            %set(fig, 'Position', [0 0 640 480])
            print(filename,'-dpng','r72');
        end
        
        % LMMSE method
    elseif (method_sel(4))
        est_channel_lmmse_time = zeros(pilots_nb_fre,NsymbsperTTI);
        est_channel_lmmse = zeros(Totactivesub_temp,NsymbsperTTI);
        for i = 1:pilots_nb_fre
            est_channel_lmmse_time(i,:) = lmmse_time( h_channel_est_zf(i,:), pilots_time_instances,NsymbsperTTI,fD_T_LMMSE,SNRs,window_t,pilot_time_spacing,enable_SW);
        end
        
        for j = 1:NsymbsperTTI
            est_channel_lmmse(:,j) = lmmse_fre(est_channel_lmmse_time(:,j),fftSize,Totactivesub_temp,h_channel_pdp,pilot_frequency_spacing,SNRs,window_f,pilots_frequency_instances,enable_SW);
        end
        
        MSE_lmmse_t = mean((abs(h_channel_fft_diag(1,:) - est_channel_lmmse(1,:))).^2);
        MSE_lmmse_f = mean((abs(h_channel_fft_diag(:,1) - est_channel_lmmse(:,1))).^2);
        lmmse_dist_mat = 10*log10(abs(h_channel_fft_diag(:,:) - est_channel_lmmse(:,:)).^2);
        MSE_lmmse = mean(mean(abs(h_channel_fft_diag(:,:) - est_channel_lmmse(:,:)).^2));
        
        MSE_lmmse_t_db = 10*log10(MSE_lmmse_t);
        MSE_lmmse_f_db = 10*log10(MSE_lmmse_f);
        MSE_lmmse_db = 10*log10(MSE_lmmse);
        
        if (mode)
            clims = [-30 -5];
            imagesc(1:1:NsymbsperTTI,1:1:Totactivesub_temp,lmmse_dist_mat,clims)
            colorbar('FontSize',12)
            title(colorbar,'in dB','FontSize',14)
            xlabel('Time Lag')
            ylabel('Subcarrier Index')
            title(['MSE: ',num2str(MSE_lmmse_db,4),'dB',', MAX Doppler f_D = ',num2str(fD,4),'Hz',', Pilot Rate R_P = ',num2str(R_Pilot_t,5),'Hz'])
            
            set(gca,'units','pixels','outerposition',[-15,-5,640,430])
            print(filename, '-dpng','-r72');           
        else
            subplot(2,1,1)
            plot(t_axis,real(h_channel_fft_diag(1,:)),'r-','linewidth',1.5)
            hold on
            plot(t_axis,real(est_channel_lmmse(1,:)),'b-.','linewidth',1.5);
            plot(t_axis_pilots,real(h_channel_est_zf(1,:)),'m*');
            xlabel('Time Lag','FontSize',14)
            ylabel('real part of H','FontSize',14)
            xlim([0 205]);
            legend({'Ori','LMMSE','Pilot'},'FontSize',14,'Location','southeast')
            title(['Time Domain, MSE: ',num2str(MSE_lmmse_t_db,4),'dB'],'FontSize',15)
            
            subplot(2,1,2)
            plot(f_axis,real(h_channel_fft_diag(:,1)),'r-','linewidth',1.5)
            hold on
            plot(f_axis,real(est_channel_lmmse(:,1)),'b-.','linewidth',1.5);
            plot(f_axis_pilots,real(h_channel_est_zf(:,1)),'m*');
            xlabel('Subcarrier Index','FontSize',14)
            ylabel('real part of H','FontSize',14)
            legend({'Ori','LMMSE','Pilot'},'FontSize',14,'Location','southeast')
            title(['Frequency Domain, MSE: ',num2str(MSE_lmmse_f_db,4),'dB'],'FontSize',15)
            
            %set(fig, 'Position', [100 100 640 480])
            print(filename,'-dpng','-r72'); % generate the picture results
        end               
    end
end


