% Simulation Of Jakes Model
clear all;
f_max = 26;
M = 7;                                 % # of low frequency oscillators - 1
N = 4*M+2;
Ts=1.024e-04;
sq = 2/sqrt(N);
sigma = 1/sqrt(2);
theta = 0;                              % Fixed Phase
count = 0;
t0=0.005;
    for t = 0:Ts:0.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
  
    len=length(lag_i);
   [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-2001)./2;%mid ... points for drawing figures
      cc=bb+1:1:bb+2001;%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(1);
   plot(tdd,auto_i(cc),'-k',tdd,T_bessel,'.m');%in-phase
   xlabel('t(Second)');
   ylabel('Auto-correlation');
   legend('In-component(in-phase)','T-bessel')
   
   figure(2);
   plot(tdd,auto_q(cc),'-k',tdd,T_bessel,'.m');%quadrature
   xlabel('t(Second)');
   ylabel('Auto-correlation');
   legend('Q-component(quadrature)','T-bessel')
 
   %figure;
   %plot(tdd,corrx2(cc));
   %xlabel('t(Second)');
   %ylabel('cross-correlation');  
figure(3);
co1=1:1000;
semilogy(co1*Ts,envelope(1:1000));
 xlabel('t(Second)');
 ylabel('Rayleigh Coef.'); 
 legend('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(4);
plot(x,pdf_env,'-k',x,pdf_theory,'.m');
legend('Simulated','Theoretic');
xlabel('r');
ylabel('PDF of r');