Commit f41cf10f authored by Giulia Illiano's avatar Giulia Illiano
Browse files

- Dalla riga 31 ho aggiunto il rebin dei tempi originali, ora chiamati t_raw 

- Parlando con Filippo potrebbe essere che non serve fare il ciclo for per calcolare i tempi medi dei vari bin di ogni segmento (ciclo che ora va dalla riga 49 alla 61). Alla riga 38 ho scritto un commento a riguardo, per capire se possiamo assegnare quei tempi medi di ogni bin proprio durante l'operazione di rebin

- alla riga 46 ho calcolato il numero di segmenti, usando la funzione fix. E' altamente improbabile che l'ultimo segmento sia esattamente lungo Tseg secondi, quindi si potrebbe scartare a priori. Altrimenti andrà aggiunto un controllo

- ANCORA DA CAPIRE IL FATTO DEI "nismin/nismax" alle righe 91-92, se sono uguali per tutti i segmenti

- Facendo finta di aver creato già una coherent template bank (MA E' DA FARE, NON SAPEVO COME), dalla riga 119 ho fatto un ciclo for dove per ogni segmento: 1) faccio la TdF 2) Scelgo la Roi (righe 128-130) e faccio la trasformata inversa (riga 131) 3) ANDRA' FATTO IL RESEMPLED DELLA SERIE TEMPORALE, MA NON SAPEVO FARLO NON AVENDO I TEMPLATES 4) faccio lo zero-padding (riga 137) 5) rifaccio la trasfromata di fourier, questa volta usando la normalizzazione di Leahy, ma da capire se va bene 6) ANDRA' CALCOLATA LA DETECTION STATISTICS
parent 44cbb315

scsearch.m

+68 −13
Original line number Diff line number Diff line
@@ -5,8 +5,8 @@ close all; 
pathfi = './';



finame = 'franco.fits';

t = fitsread([pathfi,finame],"binarytable");

t = t{1};

t_raw = fitsread([pathfi,finame],"binarytable");

t_raw = t_raw{1};

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

% for i = 1:length(info)

%     if (isequal(info{i},'TIMEDEL'))
@@ -14,27 +14,40 @@ t = t{1}; 
%     end

% end



% % Controlli sui tempi

% for i = 1:(length(t)-1)

%     if (t(i+1)-t(i)>= res*1.001)



Tseg=256; %segments' length in seconds

%NOTE: t0 used to define the orbital phase γ is equal to the mid point

%of the observation span for dataset; M2015, Sec.7

t0 = (t_raw(end)+t_raw(1))/2;



% Prendi da file il tspan

tspan = t(end)-t(1);

tspan = t_raw(end)-t_raw(1);

param = [f_min,f_max,porb_min,porb_max,a_min,a_max,tasc_min,tasc_max];



%NOTE: t0 used to define the orbital phase γ is equal to the mid point

%of the observation span for dataset; M2015, Sec.7

t0 = (t(end)+t(1))/2;



% Aux par

Omega_max = 2*pi/porb_min;

Omega_min = 2*pi/porb_max;

gamma_min = Omega_min*(t0-tasc_max);

gamma_max = Omega_max*(t0-tasc_min);



% Time matrix with bin midpoints for each segment

%Rebin -------------------------------------------------------------------

Nyq=2000; %cambiare in caso

 

dt_psd=1/(2*Nyq); %risoluzione spettro di potenza, diverso dal dt dei segmenti

nbin=round((t_raw(end)-t_raw(1))/dt_psd);

[C,t]=(histcounts(t_raw,nbin)); 

C=C.'; %conteggi

t=(t(2:end)).'; %vettore tempi rebbinato, assegno al tempo l'edge destro di ogni bin

%se serve il centro del bin devo fare: t=t(1:end-1)+(t(2)-t(1))/2; 

%perchè dt=t(2)-t(1)

%VEDERE SE E' NECESSARIO FARE TUTTO IL CICLO SOTTO PER I tm

%Come fatto qui sopra, tutti i conteggi sono assegnati al centro di ogni

%bin che dura dt



Tseg=256; %segments' length in seconds

M=fix(t(end)-t(1))/Tseg; %number of segments

%con fix prendo la parte intera, scarto l'ultimo segmento che tanto non

%sarà mai di lunghezza Tseg (molto improbabile)



% Time matrix with bin midpoints for each segment----------------------

dt = Tseg/N;

tm = zeros(M,N);

tmid = zeros(M,1);
@@ -45,6 +58,8 @@ for m = 1:M 
    end

    tmid(m) = (tm(m,N)+tm(m,1))/2;

end

%CONTROLLARE SE CON LA FUNZIONE RESHAPE PUOI EVITARE IL CICLO FOR

%PERCHE' tm dovrebbe essere diverso per i vari segmenti?????



% Time series for each segment

x = zeros(M,N);
@@ -94,6 +109,46 @@ while(1) 
    end

end



%PARTI DA FARE GUARDANDO L'OVERLEAF:

% - Create (...)a coherent template bank in the ν(m) space

% - Save this template bank, we'll need it later



%Fourier transform on original time-series --------------------------------

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



%for each segments (lavoro su tm)

for m=1:M

    [C,edges]=(histcounts(tm(m,:),round((tm(m,end)-tm(m,1))/dt_psd))); 

    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 

    Y=fft(C).'; clear C

    F=((0:length(Y)-1)./edges).'; clear edges

    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.

    

    %SCRIVERE CODICE PER RESAMPLED TIME SERIES!!!!

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

    

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

    X1 = [X1;zeros(L-length(F),1)]; %concateno alla timeseries downsapled 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 -----------------------------------------------------

    [C,edges]=(histcounts(X1,round((X1(end)-X1(1))/dt_psd))); 

    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(C).*abs(fft(C)).^2).'; %normalizzazione Leahy, giusto?????

    clear C

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

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

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



    %CALCOLARE LA DETECTION STATISTICS!!!!

end





%COMPLETARE DOPO AVER STUDIATO I RANDOM TEMP. BANKS

%E VEDERE CODICE https://pyfstat.readthedocs.io/en/latest/pyfstat.html#pyfstat.core.SemiCoherentSearch