function [channel_matrix,channel_matrix_equivalent,channel_matrix_upsampling,MIMO_OFDM_link_par]=...
    MIMO_OFDM_channel_generation(MIMO_OFDM_system_par,MIMO_OFDM_link_par,...
    channel_0_range_min,channel_0_range_max)
% Function MIMO_OFDM_channel_generation +++++++++++++++++++++++++++++++++++++++++++++++++
% Generate the channel of MIMO OFDM system

% Input data:  MIMO_OFDM_system_par - parameters of MIMO OFDM system
%              MIMO_OFDM_link_par - parameters of MIMO OFDM link

% Output data: channel_matrix - the original channel matrix 
%                                      (num_rx_antenna,num_tx_antenna,channel_length,num_frame)
%              channel_matrix_equivalent - the equivalent channel matrix including the effect of the pluse shaping
%         if pulse shaping is disabled, (num_rx_antenna,num_tx_antenna,channel_length,num_frame)
%         if pulse shaping is enabled, (num_rx_antenna,num_tx_antenna,channel_length_equivalent,num_frame)
%              channel_matrix_upsampling - the upsampling channel matrix
%         if pulse shaping is disabled, 0
%         if pulse shaping is enabled, (num_rx_antenna,num_tx_antenna,channel_length*pulse_oversampling,num_frame)
%              MIMO_OFDM_link_par - parameters of MIMO OFDM link

% Begin:  Feb, 2006
% Last edited: Feb 14, 2007
%   Feng Wan    

num_frame=MIMO_OFDM_system_par.num_frame;
num_tx_antenna=MIMO_OFDM_system_par.num_tx_antenna;
num_rx_antenna=MIMO_OFDM_system_par.num_rx_antenna;

pulse_type=MIMO_OFDM_system_par.pulse_type;%%%%% 5/9/2008 cz pulse_type = 0;  0, the pulse shaping is disabled;

if(pulse_type==1)                    %%%%%%%%%if 0, the pulse shaping is disabled; if 1, the pulse shaping is enabled
    pulse_oversampling=MIMO_OFDM_system_par.pulse_oversampling;
    pulse_filt_order=MIMO_OFDM_system_par.pulse_filt_order;
    pulse_rolloff_factor=MIMO_OFDM_system_par.pulse_rolloff_factor;
    pulse_delay=MIMO_OFDM_system_par.pulse_delay;
    channel_length_equivalent=MIMO_OFDM_system_par.channel_length_equivalent;
end

channel_model_type=MIMO_OFDM_system_par.channel_model_type;
channel_length=MIMO_OFDM_system_par.channel_length;

if(MIMO_OFDM_system_par.channel_time_varying==1)  %%% if 1, channel is time varying; otherwise it is constant. If channel_mode_type=1, Jakes model is used.
    norm_sampling_ch_freq=MIMO_OFDM_system_par.norm_sampling_ch_freq;
    doppler_shift=MIMO_OFDM_system_par.doppler_shift;
    num_harmoniic_funs=MIMO_OFDM_system_par.num_harmoniic_funs;
end

%+++++++++++++++++++++++++++++++++++++++++++ 
% Make Channel Matrix
%+++++++++++++++++++++++++++++++++++++++++++ 

if(channel_model_type==1)    % 3GPP channel model  

        
elseif(channel_model_type==2)  % random matrix model    
    
    LOOP_Sig=0;
    
    while(LOOP_Sig~=1)

    
    channelcoef=randn(channel_length,1);   % ISI channel coefficients    
    abs_channelcoef=sqrt(sum(abs(channelcoef).^2));
    stdvar_channelcoef=(channelcoef/abs_channelcoef);
    channel_matrix_original=zeros(num_rx_antenna,num_tx_antenna,channel_length,num_frame);

    channel_gain_time=zeros(channel_length,num_frame);
    
    if(MIMO_OFDM_system_par.channel_time_varying==1)

        [channel_gain_time,ar_coef_ave,ar_noise_variance_ave]=multi_jakes(channel_length,doppler_shift,...
            norm_sampling_ch_freq,num_frame,num_harmoniic_funs);
        % set the coefs of first order AR model for the estimation of Kalman filter
        MIMO_OFDM_link_par.ar_coef=ar_coef_ave;
        MIMO_OFDM_link_par.ar_noise_variance=ar_noise_variance_ave;
    else
        for n1=1:num_frame
            for n2=1:channel_length
                channel_gain_time(n2,n1)=1;
            end
        end
    end

       tempm=randn(num_rx_antenna,num_tx_antenna,channel_length)/sqrt(2)+...
          sqrt(-1)*randn(num_rx_antenna,num_tx_antenna,channel_length)/sqrt(2);
       for ntime=1:num_frame        
           for n2=1:channel_length
%                channel_matrix_original(:,:,n2,ntime)=stdvar_channelcoef(n2)*tempm(:,:,n2)*...
%                    channel_gain_time(n2,ntime); % Transpose
               channel_matrix_original(:,:,n2,ntime)=tempm(:,:,n2)*...
                   channel_gain_time(n2,ntime)/sqrt(channel_length); % Transpose               
           end
       end
       sum_norm_temp=0;
       for n1=1:channel_length
           sum_norm_temp=sum_norm_temp+norm(channel_matrix_original(:,:,n1,1),'fro').^2;
       end
       norm_temp=norm(channel_matrix_original(:,:,1,1),'fro').^2/sum_norm_temp;
       if((norm_temp>channel_0_range_min)&&(norm_temp<channel_0_range_max))
           LOOP_Sig=1;
       end

    
    
    if(pulse_type==1)          % upsampling the channel by adding zeros in frquency domain
        channel_length_upsampling=channel_length*pulse_oversampling; 

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% Assumption Here. March, 2006: make the problem easy
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        channel_no_upsampling=1;   % Don't upsampling the channel coefficients
        if(channel_no_upsampling==1)
            channel_matrix_upsampling=zeros(num_rx_antenna,num_tx_antenna,channel_length_upsampling,num_frame);
            for n1=1:num_frame
                for n2=1:num_rx_antenna
                    for n3=1:num_tx_antenna
                        for n4=1:channel_length
                            channel_matrix_upsampling(n2,n3,(n4-1)*pulse_oversampling+1,n1)=...
                                channel_matrix_original(n2,n3,n4,n1);
                        end
                    end
                end
            end
        end
        
        
        channel_matrix=channel_matrix_upsampling;   
  
        % generate the equivalent channel matrix combing the effect of pluse shaping: Method2
        channel_matrix_equivalent=zeros(num_rx_antenna,num_tx_antenna,channel_length_equivalent,num_frame);
        pulse_rrcfilter=rcosine(1,pulse_oversampling,'fir/sqrt',pulse_rolloff_factor,pulse_delay);% Create a square root raised cosine filter.
        channel_vec_temp11=1;
        channel_vec_temp12=rcosflt(channel_vec_temp11,1,pulse_oversampling,'filter',pulse_rrcfilter);
        channel_vec_temp13=zeros(channel_length_upsampling,1);
        
        for n1=1:num_frame
            for n2=1:num_rx_antenna
                for n3=1:num_tx_antenna
                    for n4=1:channel_length_upsampling
                        channel_vec_temp13(n4)=channel_matrix_upsampling(n2,n3,n4,n1);
                    end
                    channel_vec_temp14=conv(channel_vec_temp13,channel_vec_temp12);
                    channel_vec_temp15=rcosflt(channel_vec_temp14(1:length(channel_vec_temp14)-...
                        pulse_oversampling+1),1,pulse_oversampling,'Fs/filter',pulse_rrcfilter);
                    channel_vec_temp16=downsample(channel_vec_temp15,pulse_oversampling); % Downsample.
%                     channel_vec_temp17=channel_vec_temp16(2*pulse_delay+1:end-2*pulse_delay);
                    for n4=1:channel_length_equivalent
                        channel_matrix_equivalent(n2,n3,n4,n1)=channel_vec_temp16(n4);
                    end
                end
            end
        end  
        % Let the channel_matrix as the original channel matrix
        channel_matrix=channel_matrix_original;
    else    
        channel_matrix=channel_matrix_original;
        channel_matrix_equivalent=channel_matrix;
        MIMO_OFDM_link_par=MIMO_OFDM_link_par;
        channel_matrix_upsampling=0;  % no channel_matrix_upsampling
    end   
    
    end
else    
    error('input channel_model_type is not wrong!');
end