%-------------------------------------------------------
% PROGETTO DI SEGNALI 2021: STACKING
%-------------------------------------------------------
% INIZIO PARTE ASSEGNATA
%-------------------------------------------------------
clear 
close all
%GENERAZIONE DEL SEGNALE DI RIFERIMENTO
N=300;
x=randn(N,1);
%x=ones(N,1);
% FILTRAGGIO PASSA BASSO E FINESTRATURA GAUSSIANA
h=fir1(20,.15);
y=conv(x,h,'same').*gausswin(N,4);
% Si eleva al quadrato per simulare un'immagine rivelata
y=y.^2;
% SOMMA DI RUMORE GAUSSIANO CON DEVIAZIONE SIGMA TUTTO AL QUADRATO
sigma=.15;
yn=y+(randn(N,1)*sigma).^2;
figure(1)
plot(y)
grid
axis([1 N -0.2 0.8])
hold on
pause
plot(yn-y,'r')
title('segnale originale e rumore additivo')
hold off
pause
figure(2)
plot(yn)
grid
axis([1 N -0.2 0.8])
title('segnale + rumore')
pause
Y=fftshift(fft(y));

% GENERAZIONE DEI K SEGNALI TRASLATI E SOMMATI A RUMORI GAUSSIANI INDIPENDENTI
K=100;
y_shifted=zeros(N,K);
SH_dato=zeros(K,1);

for k=1:K
shift=randn(1,1)*5;
ex=exp(j*2*pi*shift/N*(-N/2:N/2-1));
YS=Y.*transpose(ex);
ys=real(ifft(ifftshift(YS)));
ysn=ys+(randn(N,1)*sigma).^2;
% memorizzazione dei segnali traslati rumorosi
y_shifted(:,k)=ysn;
% memorizzazioni degli shift rispetto al segnale di riferimento
SH_dato(k)=shift;
end
% RAPPRESENTAZIONE DEGLI SHIFT DATI RISPETTO AL RIFERIMENTO
figure(3)
stem(SH_dato,'filled')
grid
title('shift dei segnali rispetto al riferimento')
pause
% RAPPRESENTAZIONE DEI SEGNALI DISALLINEATI RUMOROSI
figure(4)
plot(y_shifted)
axis([1 N -0.2 0.8])
grid
%-------------------------------------------------------
% FINE PARTE ASSEGNATA
%-------------------------------------------------------
pause


% ------------------------------------------------------
% IL PROGETTO DA SVOLGERE
%-------------------------------------------------------
% 1 - Stimare lo shift tra i K segnali generati predentemente e memorizzati
%     nella matrice y_shifted(300,K) e il segnale di riferimento yn.
%     Per stimare lo shift occorre: a) calcolare le CROSSCORRELAZIONI TRA 
%     SEGNALE DI RIFERIMENTO E GLI ALTRI K  R=xcorr(yn,ysn);
%     (si può usare anche la FFT se i segnali sono molto lunghi)

    % b) trovare il massimo Rmax e la posizione del massimo nmax della
    % cross-correlazione (si noti che per segnali perfettamente allineati
    % il massimo della cross-correlazione si trova in posizione N
    % Se si vuole passare dalla FFT il max dell'autocorrelazione si trova
    % in posizione N/2+1
    
    % c) interpolare posizione del massimo con parabola y=ax^2+bx+c noti 
    %    3 punti intorno al massimo come spiegato di seguito:
    %    y(0)=Rmax=R(nmax), 
    %    y(-1)=R(nmax-1)  
    %    y(1)=R(nmax+1) 
    %    stima dei coefficienti della parabolo interpolante:
    %    a=y(-1)/2+y(1)/2-y(0)
    %    b=y(-1)/2-y(1)/2
    %    l'annullamento della derivata prima si ha in xmax=-b/(2a)
    % ...............................................
    % ...............................................
    % memorizzare gli shift stimati rispetto al segnale di riferimento
    % e confrontarli con quelli dati

% 2 - Riallineare i K segnali memorizzati nella matrice y_shifted(1000,K) 
%     con il segnale di riferimento yn e farne lo stack (media punto a
%     punto).  
%     Memorizzazione dei segnali riallineati loro visualizzazione e 
%     visualizzazione dello stack.
