clear
%script per determinare la correlazione e lo spettro di potenza
%del segnale numerico e del processo aleatorio sinusoidale

tc=.01;     %intervallo campionamento segnale
symbols=20; %il tempo va da 0 a symbols con passo tc
realizations=300;

pam=1;    %1 per il segnale numerico, ò per la sinusoide

%parametri segnale numerico binario antipodale
asynch=0; %1 per rendere asincrono il segnale numerico
pn=1/2;   %probabilità 0 nella sequenza numerica, associato a -1

%parametri segnale sinusoidale: vettori di dimensione realizations
amplitudes=2*(rand(realizations,1))-1;
phases=pi*rand(realizations,1);

%parametri filtraggio
W=4; %W/2=banda filtro ideale
t0=5;%ritardo filtro per la realizzazibilità

displayplots=1; %1 per visualizzare la trasformate

samples=round(1/tc);
t=tc:tc:symbols;
bits=rand(realizations,symbols)>pn;
an=2*bits-1;

offset=floor(rand(1,realizations)*samples);
x=zeros(realizations,length(t));
for h=1:realizations
    if pam
        for i=1:symbols 
            for j=1:samples 
                x(h,(i-1)*samples+j+asynch*offset(h))=an(h,i); 
            end 
        end
    else
        x(h,:)=amplitudes(h)*cos(2*pi*t-phases(h));
    end
    x=x(1:realizations,1:length(t));
    X(h,:)=fft(x(h,:))*tc;
end
[corr,rho] = correlazione(x,samples);
[corrt,rhot] = correlazionetemporale(x,2*samples);

hfilter=W*sinc(W*(t-t0))*tc;
t1=tc:tc:2*symbols-tc;
for h=1:realizations
    y(h,:)=conv(hfilter,x(h,:));
    y1(h,:)=y(h,t1>=t0 & t1<symbols+t0);
    Y(h,:)=fft(y1(h,:))*tc;
end
[corry,rhoy] = correlazione(y1,samples);
[corrty,rhoty] = correlazionetemporale(y1,2*samples);

if (displayplots)
    figure
    subplot(1,2,1)
    plot(t(samples: 10*samples),x(1,samples:10*samples),'k')
    hold on
    plot(t(samples: 10*samples),y1(1,samples:10*samples),'r')
    title('segnale numerico (nero) segnale filtrato (rosso)');
    t1=tc:tc:2*symbols-tc;
    fs=1/tc;
    ns=length(t);
    fshift=(-ns/2:ns/2-1)*(fs/ns);
    Xs=fftshift(X);
    subplot(1,2,2)
    plot(fshift,mean(abs(Xs)),'k')
    Ys=fftshift(Y);
    hold on
    plot(fshift,mean(abs(Ys)),'r')
    set(gca,'Xlim',[-3*W,3*W]);
    title('spettro segnale numerico (nero) segnale filtrato (rosso)');
    Hfilter=fft(hfilter);
    Hfs=fftshift(Hfilter);
    plot(fshift,abs(Hfs),'b')
    figure
    subplot(1,2,1)
    plot(t(1*samples:4*samples),corr(1*samples:4*samples,round(samples/4)),'k');
    title('correlazione per \tau=1/4 in funzione di t');
    xlabel('t');
    if pam==1
        set(gca,'Ylim',[-0.1,1.1]);
    else
        set(gca,'Ylim',[-1.1,1.1]);
    end
    subplot(1,2,2)
    if pam==1
        plot(t(1*samples:4*samples),corr(round(5*samples/4),1:3*samples+1),'k')
        set(gca,'Ylim',[-0.1,1.1]);
        title('correlazione per t=1/4 T0 in funzione di \tau');
    else
        plot(t(1*samples:4*samples),corr(10,1:3*samples+1),'k')
        set(gca,'Ylim',[-1.1,1.1]);
        title('correlazione per t=.1 in funzione di \tau');
    end        
    xlabel('\tau');
end