function [rt,zt,rt1,rt2]=sosfadeli(nc,fd,T,KdB,N,M)
%function [rt,zt,rt1,rt2]=sosfadeli(nc,fd,T,KdB,N,M)
% rt  fading (complex value: rt=rt1+jrt2)
% zt  power (in dB)
% rt1 in phase component
% rt2 in quadrature component
% nc samples
% the results are M x nc matrixes
% fd maximum Doppler frequency (Hz)   fd*T<=1
% T sampling time (s)
% KdB Rice factor (dB); KdB=-inf for Rayleigh
% N oscillators per sequence
% M sequences
%
% This function generates M uncorrelated fading sequences
% the amplitude follows the Rice model
% The temporal correlation follows the Clarke model
% the autocorrelation coefficient of rt1, rt2 and rt is J0(2*pi*fd*n*T)
% the autocorrelation coefficient of 10^(zt/10) is (J0(2*pi*fd*n*T))^2
% the autocorrelation coefficient of abs(rt) is appr. (J0(2*pi*fd*n*T))^2
% rt complex envelope (complex low-pass signal) received: rt=rt1+j*rt2
% zt envelope-square normalized with respect to the mean power expressed in dB  zt=(rt1^2+rt2^2)

% SOS method in Li TR-COM 2002, p.1503 (sum-of-sinusoid)
% Fulvio Babich

% if the transmitted signal is s(t)=Re(u(t)*exp(2*pi*fc*t))
% the received signal is x(t)=Re(r(t)*exp(2*pi*fc*t))
% with: u(t) complex envelope (complex low-pass signal) transmitted
%       r(t) complex envelope (complex low-pass signal) received
% rt is obtained by setting u(t)=1 i.e. considering a transmission with
% unmodulated carrier

% I convert the Rice factor from dB to amplitude

    K=10^(KdB/10);
    No=N/4;

    % creating the time base for the oscillators


    t=0:T:(nc-1)*T;

    % Calculating random phases

    a00=pi/2/M/N;
    for k=1:M
       for n=1:No
          a(n,k)=2*pi*(n-1)/N+2*pi*(k-1)/N/M+a00;
          fi(n,k)=2*pi*rand(1,1);
          fip(n,k)=2*pi*rand(1,1);
       end
    end


    % Calculation of in-phase and quadrature components

    rt1=zeros(M,nc);
    rt2=zeros(M,nc);

    for k=1:M
       for n=1:No
          rt1(k,:)=rt1(k,:)+(2/No)^.5*cos(2*pi*fd*cos(a(n,k))*t+fi(n,k));
          rt2(k,:)=rt2(k,:)+(2/No)^.5*sin(2*pi*fd*sin(a(n,k))*t+fip(n,k));
       end
    end


    % rt1 and rt2 are Gaussian random variables with zero mean and unit variance => N(0,1)

    % modification of rt1 and rt2 to account for the Rice factor (Babich)
    var(rt1);
    var(rt2);
    rt1=rt1/(K+1)^.5;
    rt2=(rt2+(2*K)^.5)/(K+1)^.5;

    % complex envelope (complex low-pass signal) received

    rt=rt1+j*rt2;

    zt=(rt1.^2+rt2.^2)/2;

    % convert the square envelope to dB
    zt=10*log10(zt);
end