function fa_comp= rayleigh(fd, fs, Ns)
%*************************************************************************************
%function fade= rayleigh(fd, fs, Ns)
% Copyright 2002 The Mobile and Portable Radio Research Group
%
% function to generate Rayleigh faded waveform
% The input to the function are
%    Doppler frequency
%    Sampling Frequency
%    Number of samples needed
% The output of the funtion is an array of complex number whose amplitude has Rayleigh 
% distribution and the phase is uniformly distributed from -pi to +pi.
% ~~~~~~~~~~~~~~~~~~~~~~~~

% Written by James Hicks
% Modified by Fakhrul Alam

% ~~~~~~~~~~~~~~~~~~~~~~~~

%function fade= rayleigh(fd, fs, Ns)
% constants:
order= 3;  % order of polynomial extrapolation 
		   % along verticle assymptote of Doppler 
		   % Spectra, SEZ.

% Create longer duration of fading envelope than needed in order
% to have a spectrum that is a power of two.
% Set the frequency spacing for the required sampling rate.
% Adjust Doppler shift to assure that assymptote is sampled.

N= 2^ceil(log2(Ns)); 
delta_f= fs/N;   % (Hz) frequency spacing
kd= ceil(fd/delta_f);

% If Doppler shift is very small, envelope will be roughly constant: no
% apparent fade.
if(kd <= order) % assure at least one point outside interpolation keys
   %fade= ones(Ns,1);
   fa_comp=ones(Ns,1)*(1+j);
else

% compute actual adjusted Doppler Shift. 
fm= delta_f*kd;

% Compute the Doppler Power Spectral Density 
SEZ= zeros(kd+1,1);
f= (0:kd-1)*delta_f;
SEZ(1:kd)= 1.5./(pi*fm*sqrt(1-(f/fm).^2));
% Use Polynomial fit to get the component at f= fm
% replace infite value at f=fm with a polynomial extrapolation.
p= polyfit( f(kd-order:kd), SEZ(kd-order:kd).', order);    
SEZ(kd+1)= polyval(p, fm);


I_RAND= [sqrt(2)*randn(1); randn(kd-1,1)+j*randn(kd-1,1); sqrt(2)*randn(1)];
Q_RAND= [sqrt(2)*randn(1); randn(kd-1,1)+j*randn(kd-1,1); sqrt(2)*randn(1)];
Io= sqrt(SEZ).*I_RAND;
Qo= sqrt(SEZ).*Q_RAND;

pad = 2^(nextpow2(2*Ns-1));

If= [Io; zeros(pad-2*kd-1,1); conj(Io(kd+1:-1:2))];
Qf= [Qo; zeros(pad-2*kd-1,1); conj(Qo(kd+1:-1:2))];
Ienv= real(ifft(If));
Qenv= real(ifft(Qf));

fade= sqrt(Ienv(Ns:2*Ns-1).^2 + Qenv(Ns:2*Ns-1).^2);
fa_comp=(Ienv(Ns:2*Ns-1)+j*Qenv(Ns:2*Ns-1))/ sqrt(mean(fade.^2));
%fade= fade/ sqrt(mean(fade.^2));

end;