function [chann,Alpha,doppler_rate]=modified_jakes(v,fc,fs,U,M,segma)

% Ali Jamoos
% Equipe Signal & Image, UMR 5131 LAPS, universit? Bordeaux1
% 351 Cours de la lib閞ation, 33405 Talence Cedex, France
% email: ali.jamoos@laps.u-bordeaux1.fr

% ----------------------------------------------------------------------------------------------------------------------
% Generation of U uncorrelated Rayleigh fading channels according to the modified Jakes model:
% [P. Dent, G. Bottomley and T. Croft, "Jakes fading model revisited," IEE Electronics Letters, pp.1162-1163, June 1993]
% ----------------------------------------------------------------------------------------------------------------------

% Version (2)
% Last updated 11/02/2005

% v:  Vehicle speed in kmph
% fc: Central carrier frequency in MHz
% fs: Symbol frequency in kbps
% U:  Number of sub-carriers or sub-channels
% M:  Number of channel coefficients to generate
% segma: Variance of the channel coefficients
% chann: channels coefficients matrix of size U by M
% fm: Maximum doppler frequency in Hz
% doppler_rate: Doppler rate or fading rate


% Example:	
% --------
% [chann,fm,doppler_rate]=modified_jakes(100,2000,10,3,10000,0);

    fm=(1e3/3600)*v*fc/3e2; % maximum doppler frequency in Hz
    doppler_rate=fm/(fs*1000);
    No=24; % number of distinct oscilators
    NN=4*No; % total number of osillators
    omega_m=2*pi*fm; % maximum doppler frequency
    H=hadamard(No); % Walsh_Hadamard codes, to generate uncorrelated channels
    
    for t=1:No
        theta(t)=rand*2*pi;  % Randomly generated initial phases  
    end
    
    for u=1:U
        sum=0;
        t=1/(fs*1000):1/(fs*1000):(1/(fs*1000))*M; % sampling frequency
        for n=1:No
            omega_n=omega_m*cos(2*pi*(n-0.5)/NN);
            sum=sum+H(u,n)*(cos(pi*n/No)+j*sin(pi*n/No))*cos(omega_n*t+theta(n)); % sum of No distinct oscillators
        end
        T=sqrt(2/No)*sum;
                     
        T=repmat(10.^(segma/20),1,M).*T; % Define the variance of the channel (segma) in dB.
       
        chann(u,:)=T*exp(j*2*pi*u*1000/fs); % the u'th sub-channel
    end
    Alpha=besselj(0,2*pi*fm/fs/1e3); 
     
    % Plot the Rayleigh envelope for the first channel