fold_all.m

%fold PSRB0329+54 data from Acre Road
clear;
numfiles = 0;
showplots=false;

secboundary = false; %whether to fold on seconds rather than the pulsar period
shownoise = true; %show the noise plot
sample_interval = 2e-3; % nominally 2ms
samprate = 1/sample_interval;% samples per second, nominally 500
resolution=2e-3; %phase resolution of folded profile in seconds
numbins = 357; %number of folded bins (= period/resolution)


f_ephem = load('0329-freq-2014.txt'); % load the ephemeris

files = dir('data/*preproc.dat'); % just work through the preprocessed data files, here identified by the "preproc" in the name.
for file = files' % look through all the files
    datafile = file.name;
    C= strsplit(datafile,'-');
    gps_start_time = str2num(C{1}); 
    if gps_start_time < 1082269300 % before I sorted out the timestamps
        delay =  0.305-0.04+0.008;
    else
        delay = 0;
    end;
    gps_start_time = gps_start_time+delay; 
    datenum_start_time = (gps_start_time-16)/3600/24+datenum('1980-01-06 00:00');
    startstr =  datestr(datenum_start_time);
    
    fileinfo = dir(strcat('data/',datafile));
    filesize = fileinfo.bytes;
    hours = floor(filesize/4/samprate/3600);
    obssamples = hours*3600*samprate;
    
    profile_file = strcat('profiles/',num2str(gps_start_time+delay, '%11.6f'),'.pro'); % the  name of the corresponding profile file
    if ~exist(profile_file,'file') % generate a folded profile if one doesn't already exist
        fprintf('generating profile %s (%s)\n',profile_file,startstr);
        numfiles = numfiles+1;
        % generate the frequency list based on the start gps time. fr contains the
        % start frequency for each observed hour, interpolated from the reference
        % file.
        fr = zeros(16,1);
        for i=1:16
            fr(i) =  interp1(f_ephem(:,1),f_ephem(:,2),gps_start_time+(i-1)*3600);
        end;
        
        
        pulse = zeros(numbins,1);
        count = zeros(numbins,1);
        %load preprocessed file
        fid = fopen(strcat('data/',datafile),'r');
        phi0=0;
        fprintf('hour ');
        weight = zeros(16,1);
        noise = zeros(16*3600,1);
        for h = 1:4%hours
            fprintf('%d ',h);
            f = fr(h);
            fdot = (fr(h+1)-fr(h))/3600;
            if secboundary
                f=1;fdot=0;
            end;
            pulse_period = 1/f;
            data = single(transpose(fread(fid,obssamples/hours,'float')));
            weight(h) = (1/std(data))^2;
            for m=1:3600
                st = (m-1)*500;
                noise((h-1)*3600+m)=std(data((st+1):(st+500)));
            end;
            le=length(data);
            T =le/samprate;
            
            if ~secboundary
                data(data>10 | data<-10)=0; %identify noisy data by setting it to zero
            end;
            
            t=zeros(3600*samprate,1);
            for i=1:3600*samprate
                t(i) = sample_interval*i; % time now resets every hour, and the offset is handled in phi0
            end;
            
            phase = 1+floor(numbins*( (phi0+f*t+fdot.*t.*t/2) - floor(phi0+f*t+fdot.*t.*t/2)) );
            phi0 = phi0+ f*3600+fdot*3600^2/2; %initial phase for the next hour
            
            for i=1:3600*samprate
                if  (data(i) ~= 0)
                    w = (1/noise(floor(1+i/60/samprate+(h-1)*3600)))^2; %calculate the weight
                    pulse(phase(i))= pulse(phase(i))+w*data(i);
                    count(phase(i))= count(phase(i))+w;
                end;
            end;
        end;
        fprintf('\n');
        fclose(fid);
        
        if shownoise
            sample_interval = 1/(24*60.0*60);
            for i=1:length(noise)
                taxis(i) = datenum_start_time+sample_interval*i;
            end;
            figure;
            plot(taxis,noise,'b');
            dynamicDateTicks();
        end;
        
        % move the pulse to the centre
        profile = pulse./count;
        rolled_profile = zeros(length(profile),1);
        if ~secboundary
            for i=1:numbins
                j =-1+ i-(floor(numbins/2)-find(profile==max(profile)));
                if j<1
                    j=j+numbins;
                end;
                if j>numbins
                    j=j-numbins;
                end;
                rolled_profile(i) = profile(j);
            end;
        else
            rolled_profile = profile;
        end;
        
        
        %estimate the noise floor and subtract it.
        offset = (mean(rolled_profile(1:160))+mean(rolled_profile(195:end)))/2;
        rolled_profile = rolled_profile-offset;
        profile = profile-offset;
        if showplots
            if shownoise
                figure;
            end;
            stairs(linspace(0,pulse_period,numbins), rolled_profile,'b');
            xlim([0,pulse_period]);
            ylim([0.25*min(rolled_profile),1.05*max(rolled_profile)]);
            xlabel('pulse phase (seconds)');
            ylabel('approx. antenna temperature (K)');
            title(['PSR B0329+54 ',' start ',startstr,'  -  ',num2str(h),'h integration (Acre Rd)']);
        end;
        %
        % write out the profiles (unrolled)
        dlmwrite(strcat('profiles/',num2str(gps_start_time+delay, '%11.6f'),'.pro'),profile');
        
    end; % if no profile file
end; % over all preprocessed files
if numfiles==1
    fprintf('finished folding %d dat file\n',numfiles);
else
    fprintf('finished folding %d dat files\n',numfiles);
end;