%gives the complex fade coefficients
%output power is approximately unity
%inputs:max_doppler_freq (Hz); sample_freq (Hz); num_samples, start_time
%(sec)
%outputs:alpha_complex (fades), end_time (sec)

function [alpha_complex, end_time] = jakes(max_doppler_freq, sample_freq, ...
    num_samples, start_time)
   
    num_generators = 32;
    T = 1/sample_freq;  %sample time
    N = 2*(2*num_generators+1.0);
    wm = max_doppler_freq*2*pi;
    
    for k = 1:num_samples
%         t = k*T+10000; %start at some non-zero point ;
        t = (k-1)*T+start_time;
        I = 0;
        Q = 0;
        for nn = 1:num_generators
            
            %calculate input signal
            wn = wm*cos(mod((2*pi*nn)/N,2*pi));
            input_signal = cos(mod((wn*t),2*pi));
                        
            %calculate Q channel and sum
            beta_n= (pi*nn)/(num_generators+1);
            Q = Q + 2*sin(beta_n)*input_signal;
            
            %calculate I channel and sum
            I = I + 2*cos(beta_n)*input_signal;
        end
        
        %added input
        added_input_signal = (1/sqrt(2))*cos(mod((wm*t),2*pi));

        %Q Channel added input signal
        alpha = pi/4;
        added_Q_output = (2*sin(alpha))*added_input_signal;
        Q = Q + added_Q_output;

        %I Channel added input signal
        alpha = pi/4;
        added_I_output = (2*cos(alpha))*added_input_signal;
        I = I + added_I_output;
            
      alpha_complex(k) = (I + j*Q);
      
      %since E{alpha_re^{2}}=No+1 and E{alpha_im^{2}}=No
      %if you normalize it with respect to sqrt(2No+1) as below 
      %E{alpha_re^2}+E{alpha_im^2}=1
      alpha_complex(k) = (I + j*Q)/sqrt(2*num_generators+1);
 
      end_time = t;
      
      end

     %check the mean power level of output
     %real part power 
%      alpha_re=(real(alpha_complex)); 
%      alpha_im=(imag(alpha_complex));
% 
%      mean(alpha_re.^2);
%      mean(alpha_im.^2);