function [H] = jakes(Datalen,chanpf,fm,fs,No,startp)
% Datalen the length of the data
% chanpf path numbers
% fm max doppla
% startp used for continuing calculating scenario
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% for debug;
% clear all;
% clc;
% Datalen= 5000;
% chanpf=1;
% fs=1e5;
% fm=50;
%%%%%%%%%%%%%%%
if (~exist('startp')),
startp = 0;
end
if (~exist('No')),
startp = 0;
end
if startp == 0;
flag = 0;
else
flag = 1;
end
%pathnum = length(chanpf);
pathnum = chanpf;
%Datalen = Datalen + 500;
%clear all
%Datalen = 500000;
%pathnum = 6;
%v = 30 ; % km/h
%fm = 300;
%Fc = 2e9; %carrier frequency
%wc = 2*pi*Fc; %radian frequency
%W = 4e6;
Dopplershift=fm;
%fm = 500;
wm = 2*pi*fm;
%W = 4000*wm;
Ts = 1/fs; %sampling period
%Ts = 0.002;
% the problem which the amplitued is not Rayleigh distribution caused by the value of Ts.
% wm = 2*pi*fm when fm = 10 wm = 62.8 T = 0.1. when fm = 100 wm = 628 T = 0.01
% Ts = 1/384e3 = 2.6e-6 10000*Ts = 2.6e-2.(total time). The number of period is 0.026/0.0159 = 1.6
% the simulation shows when v = 600, fm = 1000, Uder the conditions Ts = 1/384e3 The number of period is 160.
% Then good result can get(rayleigh distribution). So when the v = 6km/h and Ts = 1/384e3, the datalengh should be 1000000
% When Ts = 0.002, 10000*Ts = 20. 20/0.0159 = 1250th period added July 10
% added at Nov 18 2003
% Ts = 1/(8*wm) 10000*Ts/T = 10000/8 = 1250 periods = 10000*fm/f. More peirods, the
% amplitude dirtribution is more like rayleigh distribution and the phase
% distribution is more like united distribution
t = 0:TsDatalen-1)*Ts; %time array
t = t + startp*Ts;
No = 10;
N = (2*No+1)* 2;
% alpha = pi/4;
alpha = 0;
belta = pi*(1:No)/(No+1);
w = wm * cos(2*pi*(1:No)/N);
if flag == 1;
rand('state',0); % test the rand function will output same sequence because the state is set to the same value
end
phase_init =rand(No,1)*2*pi;% [0 0 0 0 0 0 0 0 ]';%rand(No,1)*2*pi;
for k = 1:No
resin = zeros(No,length(t));
resqu = zeros(No,length(t));
for i = 1:No
resin(i,:) = 2*cos(belta(i))*cos(w(i)*t + phase_init(k));
resqu(i,:) = 2*sin(belta(i))*cos(w(i)*t + phase_init(k));
end
res2in = sum(resin) + sqrt(2) * cos(alpha) * cos(wm*t);
res2qu = sum(resqu) + sqrt(2) * sin(alpha) * cos(wm*t);
y(k,:) = sqrt(1/(2*No+1)).*(res2in + sqrt(-1)* res2qu);
%y(k,:) = y(k,:)/exp(k-1); % make it suitable with practice. The power of otherpath is less than the main path.
end
% normalization the channel
for i = 1:No
y(i,:) = y(i,:)/norm(y(i,:))*sqrt(Datalen);
%y(i,:) = y(i,:)/sqrt(exp(i-1));
end
H = y(1:pathnum,1:end);
return;
%mean(abs(H(:,:)))
%%%%%%%%%%%%%%%%%%%%%% plot
%for i = 1:10
% h(:,i) = H(:,(i-1)*50+1);
%end
%figure;
%for i=1:pathnum
%subplot(pathnum,1,i);
% plot(abs(H(i,1:20)));%/norm(H(i,:))));
%end
%hold on;
figure;
channel(1,:)=fft([H(1,1:10) zeros(1,118)]);
%power(1,:)=10*log(abs(channel(1,:)));
%plot(power(1,:));
plot(abs(channel(1,:)));
%figure;
%plot(abs(H(1,11:20)));
%figure;
%plot(abs(H(1,1:20)));
%hold on;
Datalen-1)*Ts; %time array