function [rt,zt,rt1,rt2]=sosfadeli(nc,fd,T,K,N,M)
%function [rt,zt,rt1,rt2]=sosfadeli(nc,fd,T,K,N,M)
% rt  fading (valore complesso: rt=rt1+jrt2)
% zt  potenza (in dB)
% rt1 componente in fase
% rt2 componente in quadratura
% nc numero campioni
%    i risultati sono matrici di dimensione M*nc
% fd massima frequenza Doppler (Hz)   fd*T<=1
% T periodo di campionamento (s)
% K fattore di Rice (dB); K=-inf for Rayleigh
% N numero oscillatori del simulatore per canale
% M numero canali
%
% questa funzione genera M sequenze di fading incorrelate fra loro
% l'ampiezza segue il modello di Rice
% la correlazione temporale segue il modello di Clarke
% il coefficiente di autocorrelazione di rt1, rt2 e rt è J0(2*pi*fd*n*T)
% il coefficiente di autocorrelazione di 10^(zt/10) è (J0(2*pi*fd*n*T))^2
% il coefficiente di autocorrelazione di abs(rt) è appr. (J0(2*pi*fd*n*T))^2
% rt inviluppo complesso (segnale complesso passa-basso) ricevuto: rt=rt1+j*rt2
% zt inviluppo-quadrato normalizzato rispetto alla potenza media espresso in dB  zt=(rt1^2+rt2^2)

% Simulatore di canale affetto da Fading con il metodo SOS, Li TR-COM 2002, p.1503 (sum-of-sinusoid)
% Realizzazione: Fulvio Babich

% se il segnale trasmesso e' s(t)=Re(u(t)*exp(2*pi*fc*t))
% il segnale ricevuto e' x(t)=Re(r(t)*exp(2*pi*fc*t))
% con: u(t) inviluppo complesso (segnale complesso passa-basso) trasmesso
%      r(t) inviluppo complesso (segnale complesso passa-basso) ricevuto
% rt e' ottenuto ponendo u(t)=1 ossia considerando una trasmissione con
% portante non modulata

% converto il fattore di Rice da dB in ampiezza

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

    % creazione della base dei tempi per gli oscillatori
    % T periodo di campinamento

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

    % calcolo delle fasi casuali

    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


    % calcolo delle componenti in fase e in quadratura

    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 e rt2 risultano variabili aleatorie gaussiane con media nulla e varianza unitaria => N(0,1)

    % modifica di rt1 e rt2 per tenere conto del fattore di Rice (Babich)
    var(rt1);
    var(rt2);
    rt1=rt1/(K+1)^.5;
    rt2=(rt2+(2*K)^.5)/(K+1)^.5;

    % inviluppo complesso (segnale complesso passa-basso) ricevuto

    rt=rt1+j*rt2;

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

    % converto l'inviluppo-quadrato in dB
    zt=10*log10(zt);
end