Versione essenzialmente identica a quella girata da Filippo; ~150 ks PER...

a_tru = 0.343356; %lt-s

tasc_tru = 58107.009477; %MJD?

f_gr = zeros(5,1);

%f_gr = zeros(5,1);

porb_gr = zeros(5,1);

a_gr = zeros(5,1);

tasc_gr = zeros(5,1);

tic

f_gr=f_tru+(-2:2).';

for j = 1:5

    f_gr(j) = f_tru*((j^3)/9);

    %f_gr(j) = f_tru*((j^3)/9);

    porb_gr(j) = porb_tru*((j^3)/9);

    a_gr(j) = a_tru*((j^3)/9);

    tasc_gr(j) = tasc_tru*((j^3)/9);

tic

aux1 = M*N;

aux2 = C(1:aux1);

x = reshape(aux2,[M,N]);

aux2 = C(1:aux1); clear aux1

x = reshape(aux2,[M,N]); clear aux2

toc

tic

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Lasciarlo così per circa 13E+6 bin temporali richiede ~88 ore (su questo

%% computer): trovare modo più sensato

%for m = 1:M

%    for j = 1:N

%        pippo = t>=tm(m,j)-dt/2 & t<tm(m,j)+dt/2;

% Try s* and check \nu_s range

g_jj=((pi*Tseg)^2)/3.*[1; (Tseg^2)/60; (Tseg^3)/1344; (Tseg^4)/172800]; %eq. 22 M2015 + calcoli da eq. 21 M2015

mu_s=0.5; %massimo mismatch sulla griglia coerente da scegliere

%% Troppo alto mi ma non abbiamo RAM adesso

%% Troppo alto μ ma non abbiamo RAM adesso

%% Magari spezzare a blocchi i nibank e poi cercare la minima distanza 

%% in ciascuno?

%s_s=uint8(4);

    tic

 %   [Cm,edges]=(histcounts(x(m,:),round((tm(m,end)-tm(m,1))/dt_psd))); % R - convincitene

 %   edges=edges(end)-edges(1); %mi dà il tempo preciso di tutta la TdF, che sarà leggermente diversa da length(C)*dt per come è definito histcounts 

    Y=fft(x(m,:)).'; 

    F=((0:length(Y)-1)./(tm(m,end)-tm(m,1))).'; %freq. da 0 a 2Nyq.

    L=length(F); %lunghezza iniziale, servirà per lo zero-padding

    %Choose freq. region of interest (RoI) and inverse-fourier transf. over RoI

    cond = F>=f_min & F<=f_max;

    F=F(cond);

    Y=Y(cond); 

    X=ifft(Y); %inverse-fourier transf.

    % Y=fft(x(m,:)).'; 

    % F=((0:length(Y)-1)./(tm(m,end)-tm(m,1))).'; %freq. da 0 a 2Nyq.

    % L=length(F); %lunghezza iniziale, servirà per lo zero-padding

    % 

    % %Choose freq. region of interest (RoI) and inverse-fourier transf. over RoI

    % cond = F>=f_min & F<=f_max;

    % F=F(cond);

    % Y=Y(cond); 

    % X=ifft(Y); %inverse-fourier transf.

    toc

    %FARE IL RICAMPIONAMENTO!!!!

    % Dev'essere fatto per ciascun template, ergo mettiamo un bel ciclo 

    % sulle possibili combinazioni di ni_s

    for i = 1:length(nibank)

    tic

    for i = 1:length(nibank)

        % Tento ricampionamento (Eq. 17 MP2015)

        % Per prima cosa, generiamo la nuova coordinata temporale per 

        nizero = nibank(i,s_s+1);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        tau = zeros(N,1);

        for j = 1:N

            % ta(j) = ; No Tau_zero: indifferente per il calcolo della FFT   

            for s = 1:s_s

                tau(j) = tau(j) + (nibank(i,s)/(nizero*factorial(s)))*(tm(m,j)-tmid(m))^s;

            end

        end

        X1 = interp1(tm(m,:),X,tau,'linear',0);

        % tau = zeros(N,1);

        % tic

        % for j = 1:N

        %     % ta(j) = ; No Tau_zero: indifferente per il calcolo della FFT   

        %     for s = 1:s_s

        %         tau(j) = tau(j) + (nibank(i,s)/(nizero*factorial(s)))*(tm(m,j)-tmid(m))^s;

        %     end

        % end

        tau = sum((nibank(i,1:s_s)./(nizero*factorial(1:s_s))).*((tm(m,1:N)-tmid(m)).^((1:s_s).').'),2);

        X1 = interp1(tm(m,:),x(m,:),tau,'linear',0);

        %X1 è la timeseries ricampionata (controllare che sia un vettore colonna)

        %zero-padding (metto gli zeri alla fine) -------------------------------

        X1 = [X1;zeros(L-length(F),1)]; %concateno alla timeseries downsampled and resampled gli zeri

        % X1 = [X1;zeros(L-length(F),1)]; %concateno alla timeseries downsampled and resampled gli zeri

        %il secondo è un vettore colonna di zeri con lunghezza pari al numero di elementi che mi serve per tornare alla risoluzione originale

        %Fourier transform -----------------------------------------------------

        Y1 = fft(X1).';

        clear X1

        F1=((0:length(Y1)-1)./(tm(m,end)-tm(m,1))).';

        F1=F1(1:round(length(F1)/2));

        Y1=Y1(1:round(length(Y1)/2));

        % F1=F1(1:round(length(F1)/2));

        % Y1=Y1(1:round(length(Y1)/2));

        cond = F1>=f_min & F1<=f_max;

        F1=F1(cond);

        Y1=Y1(cond);

        %Calcolo della detection statistic

        %QUESTA SARA' DA INIZIALIZZARE ALL'ESTERNO tipo Lambda = zeros(length(template),1);

        Lambda(i,m)=sum(abs(Y1).^2)/sum(x(m,:)); %CREDO (oppure prendono la potenza massima?)

        toc

    end

    toc

    disp('Ciao');

end