function [X] = Rayleigh_sub(N, f, T, M)
% Returns a simulated Rayleigh distribution [2, N] representing the time
% variable with complex value of a distribution with N samples.
% f is the max Doppler frequency.
% T is the time period of sampling.
% M is a constant, usually 8 or 16.
%
% 03-08-07 D.Bozarth, Engineering Science Dept, Sonoma State University. 

% Initialize the return value array.
X = zeros(2, N); % Col 1: time value, Col 2: complex value

% Generate a uniform RV representing a constant angle for each sample.
theta = rand(1, 1) * 2 * pi;

    % Generate uniform RV's representing constant angles for each summation point.
    alpha = rand(1, M) * 2 * pi;
    beta  = rand(1, M) * 2 * pi;

% Obtain each time value, and each corresponding complex value.
for n = 1:N % Each sample
    
    % Calculate the time value corresponding to this sample.
    Tn = T * (n - 1);
    X(1, n) = Tn;
    
    % Calculate a single sample complex value.
    m = 1:M;
    z = 2 * pi * f * cos(( (2 * m - 1) * pi + theta ) / (4 * M) ) * Tn;
    Y(m) = cos(z + alpha(m)) + i * sin(z + beta(m));
    X(2, n) = sum(Y) / sqrt(M);
end

% end of program
