%public domain

function [x, F] = gen_rayleigh(fd, fs, N, m)
%Generate Rayleigh fading coefficients. 
%
%[x, F] = gen_rayleigh(fd, fs, N, m)
%   fd is the Doppler spread to use
%   fs is the filter's sampling frequency
%   N is the number of coefficents
%   m is the number of paths
%
%   x is the coefficients' vector
%   F is the filter used
%
%The coefficient are generated using the white-noise filtering approach, which
%gives you Jakes' spectrum. The algorithm used is described in
%
%   The Generation of correlated Rayleigh Random Variates by Inverse Discrete
%   Fourier Transform
%   David. J Young, Norman C. Beaulieu
%   IEEE Tran. Comm. vol. 48 no 7 Jul 2000


    %XXX: Si volem una resolució en freqüència de 100Hz necessito 20*10^6 punts
    %per a la FFT.

    debug_on_warning=0;

    %sanity checks
    if (fd<1) x=ones(m,N); return; end;
    if (fs<2*fd) error('fs < 2*fd'); end;

    fdn=fd/fs; km=floor(fdn*N);
    if (km == 0) disp('uh-oh'); end;
    F=zeros(m,N);
    F(:,1)=0;                   %this could be Riciean instead
    for k = 2:1:km
        F(:,k) = sqrt(1 / (2 * sqrt(1 - ((k - 1) / fdn / N)^2)));
    end
    k=km+1; F(:,k)=sqrt(km/2*(pi/2-atan((km-1) / (sqrt(2*km-1)))));
    % això sobra, ja està inicialitzat abans
    %for k=km+1:1:N-km-1
    %   F(1,k)=0;
    %end
    k=N-km+1; F(:,k)=F(:,km+1);
    for k=N-km+2:1:N;
        F(:,k)=sqrt(1/(2*sqrt(1-((N-k+1)/(N*fdn)).^2)));
    end
    F=F/sqrt(sum(F(1,:).^2)/N);
    if (sum(F.^2)/N > abs(1.0 + 100*eps))
        warning ('Power of filter is not 1');
        disp(sum(F.^2)/N);
    end
    A = wgn(m, N, 1); B = wgn(m, N, 1);
    X = A.*F-j*B.*F;
    x = ifft(X, N, 2);
    return
end