% Simulation Of Jakes Model
clear all;
f_max =10;
M = 8;                                 % # of low frequency oscillators - 1
N = 4*M+2;
Ts=1e-03;
sq = 2/sqrt(N);
sigma = 1/sqrt(2);
theta = 0;                              % Fixed Phase
count = 0;
t0=0.35; %(0.24, 0.45);
    for t = 0:Ts:3            % Varying time
              
        count = count + 1
        g(count) = 0;
        for n = 1 : M+1,
            if n <= M
               c_q(count,n) = 2*sigma*sin(pi*n/M);  % Gain associated with quadrature component
               c_i(count,n) = 2*sigma*cos(pi*n/M);  % Gain associated with inphase component
               f_i(count,n) = f_max*cos(2*pi*n/N);  % Discrete doppler frequencies of inphase component
               f_q(count,n) = f_max*cos(2*pi*n/N);  % Discrete doppler frequencies of quadrature component
           else
               c_i(count,n) = sqrt(2)*cos(pi/4);
               c_q(count,n) = sqrt(2)*sin(pi/4);
               f_i(count,n) = f_max;
               f_q(count,n) = f_max;
           end;  % end if 
               g_i(count,n) =  c_i(count,n)*cos(2*pi*f_i(count,n)*(t-t0) + theta); % Inphase component for one oscillator
               g_q(count,n) =  c_q(count,n)*cos(2*pi*f_q(count,n)*(t-t0) + theta); % Quadrature component for one oscillator       
       end;   %end n
       tp(count) =  sq*sum(g_i(count,1:M+1));% Total Inphase component 
       tp1(count) = sq*sum(g_q(count,1:M+1));% Total quadrature component   
   end;  % end count no n again
   envelope=sqrt(tp.^2+tp1.^2);
   rmsenv=sqrt(sum(envelope.^2)/count);

   [auto_i,lag_i] = xcorr(tp,'coeff') ;   % Auto-correlation associated with inphase component
   [auto_q,lag_q] = xcorr(tp1,'coeff');  % Auto-correlation associated with quadrature component
  %lpa1=lag_i(1:10)
    len=length(lag_i)
    %pause
   [corrx2,lag2] = xcorr(tp,tp1,'coeff');
% Cross Correlation between inphase and quadrature components
      aa=-(len-1)/2:1:(len-1)/2;%total duration for lag
      bb=(len-4001)./2;%mid ... points for drawing figures
      cc=bb+1:1:bb+4001;%for getting the mid-values 
      dd=-2000:1:2000;
  %_______________
  tdd=dd*Ts;
  z=2.*pi.*f_max*tdd;
  sigma0=1;
  T_bessel=sigma0.^2.*besselj(0,z);   
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   figure;
   plot(tdd,auto_i(cc),'-k',tdd,real(T_bessel),'.y');%in-phase
   xlabel('t(Second)');
   ylabel('Auto-correlation');
   legend('I-Simulated value','I-Theoretic value');
      figure;
   plot(tdd,auto_q(cc),'-k',tdd,real(T_bessel),'.y');%quadrature
   xlabel('t(Second)');
   ylabel('Auto-correlation');
   legend('Q-Simulated value','Q-Theoretic value');
  
figure
co1=1:1000;
semilogy(co1*Ts,envelope(1:1000));
 xlabel('t(Second)');
 ylabel('Rayleigh Coef.');
%%__________________________________________
length_r=length(envelope);
%____________________________
pdf_env=zeros(1,501);
count=0;
temp=round(100.*envelope);
for k=1:length_r
    if temp(k)<=500
        count=count+1;
        pdf_env(1,temp(k)+1)=pdf_env(1,temp(k)+1)+1;
    end
end    
count
pdf_env=pdf_env./count./0.01;
sgma2=0.5;
x=[0:0.01:5];
pdf_theory=(x./sgma2).*exp(-1.*x.^2./(2.*sgma2));
    
figure
plot(x,pdf_env,'-k',x,pdf_theory,'.y');
legend('Simulated value','Theoretic value');
xlabel('r');
ylabel('PDF of r');