jakes.m

jakes信真程序

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=1000;
% chanpf=5;
% fs=3.84e6;
% fm=100;
%%%%%%%%%%%%%%%
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:Ts:(Datalen-1)*Ts; %time array
t = t + startp*Ts;
No = 8;
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;
%%%%%%%%%%%%%%%%%%%%%% plot
% figure;
% for i=1:pathnum
%    subplot(pathnum,1,i);
%    plot(abs(h(i,:)/norm(h(i,:))));
%  end