Commit 7c25fa50 authored by Riccardo La Placa's avatar Riccardo La Placa
Browse files

Joined forces

parent e72836be

scsearch.m

+17 −7
Original line number Diff line number Diff line
@@ -15,7 +15,7 @@ t_raw = t_raw{1}; 
t_raw=t_raw./86400+50814;

MJDREF=t_raw(1);

t_raw=(t_raw-MJDREF).*86400;

t_raw=t_raw(t_raw>=t_raw(1) & t_raw<t_raw(1)+1200);

% t_raw=t_raw(t_raw>=t_raw(1) & t_raw<t_raw(1)+1200);



% info = fitsinfo(finame).BinaryTable.Keywords;

% for i = 1:length(info)
@@ -121,7 +121,7 @@ aux1 = M*N; 
aux2 = C(1:aux1); clear aux1

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

toc



clear C

% 

% Add step to generate par. extrema

param = [f_min,f_max,porb_min,porb_max,a_min,a_max,tasc_min,tasc_max];
@@ -261,6 +261,8 @@ tic 
% vari ni1,ni2,...,nis_s

%%%[kung,fu,fight] = ndgrid(nis{1},nino{2},nino{3});

% tm=gpuArray(tm);

% nibank = gpuArray(combinations(nis{:}).Variables);

% Lambda = zeros(length(nibank),length(f_gr),M,'gpuArray');

nibank = combinations(nis{:}).Variables;

Lambda = zeros(length(nibank),length(f_gr),M);

% Lambda = zeros(length(nibank),length(f_gr),M,'gpuArray');
@@ -268,6 +270,7 @@ toc 
%Fourier transform on original time-series --------------------------------

%per mantenere l'informazione di fase, non faccio il valore assoluto al quadrato della fft

%for each segment (lavoro su tm)

disp(M);

for m=1:M



    % tic
@@ -282,6 +285,9 @@ for m=1:M 
    % F=F(cond);

    % Y=Y(cond); 

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



    % ttemp = gpuArray(tm(m,:)-tmid(m));

    % xtemp = gpuArray(x(m,:));

    ttemp = tm(m,:)-tmid(m);

    xtemp = x(m,:);
@@ -316,7 +322,7 @@ for m=1:M 
            % toc



            % tic

            X1 = interp1(ttemp(:),xtemp(:),tau,'linear',0);

            % X1 = interp1(ttemp(:),xtemp(:),tau,'linear',0);

            % toc



            %X1 è la timeseries ricampionata (controllare che sia un vettore colonna)
@@ -329,9 +335,10 @@ for m=1:M 
            %edges=edges(end)-edges(2); %mi dà il tempo preciso di tutta la TdF, che sarà leggermente diversa da length(C)*dt per come è definito histcounts

            % Y1=(2./sum(X1).*abs(fft(X1)).^2).'; %normalizzazione Leahy, giusto?????

            % tic

            Y1 = fft(X1).';

            clear X1

            F1=((0:length(Y1)-1)./(ttemp(end)-ttemp(1))).';

            Y1 = fft(interp1(ttemp(:),xtemp(:),tau,'linear',0)).';

            % clear X1

            %F1=((0:length(Y1)-1)./(ttemp(end)-ttemp(1))).'; 

            F1=((0:length(Y1)-1)./(length(Y1)*dt_psd)).';

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

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

            % toc
@@ -342,13 +349,14 @@ for m=1:M 


            %Calcolo della detection statistic

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



            % toc

        end

        % disp('1ni')

    end

    toc



    disp('Ciao');

    disp(m);



end
@@ -357,6 +365,7 @@ tic 
%% The line above will have to be substituted with something along the lines

%% of the ones below to account for the metric g_jj in the phase derivatives

gdistance = @(a,b)sqrt(((a-b).^2)*(g_jj(1:s_s)));

% nisearcher = ExhaustiveSearcher(gather(nibank),'Distance',gdistance);

nisearcher = ExhaustiveSearcher(nibank,'Distance',gdistance);

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -394,6 +403,7 @@ end 
toc



disp(bestpar);

% save('C:\Users\Filippo\Desktop\XMM_Jxxx\risultelli.mat');

save('risultelli.mat');



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