temp_script_event_DFoF.m

% compute dF/F
    n_events = numel(event_peakframes);
    window = round(1/data_rec(i).metadata.KineticCycleTime * 15); % 1500 at 100 Hz % median filter (fast implementation compiled in C)
    if rem(window,2) == 0
       window = window +1; 
    end
    f = 1/data_rec(i).metadata.KineticCycleTime; % Hz, frequency data aquisition
    f1 = 1000; % 1 frame every 1 ms, upsampling frequency
    event_duration = round(1/data_rec(i).metadata.KineticCycleTime * 5); % window of 5 sec centred around the peaktime of an event
    if rem(event_duration,2) == 0
       event_duration = event_duration +1; 
    end
    DF = zeros([event_duration*n_events,size(pixels_select)]);
    DFoF = zeros([event_duration*n_events,size(pixels_select)]);
%     sDFoF = zeros([size(data,3),size(pixels_select)]);
%     int_dFoF = zeros([length([1:f/f1:size(data,3)]),size(pixels_select)]);
%     tot_length = size(data,3);
    data1 = reshape(data,[size(data,1)*size(data,2),size(data,3)]);
    for e = 1:n_events
        frame_vector = (event_peakframes(e)-floor(event_duration/2)) : (event_peakframes(e)+floor(event_duration/2)); 
        time_vector = frame_vector(1)-floor(window):frame_vector(end)+floor(window);
        centrepoint = frame_vector(floor(event_duration/2)+1);
        tot_length = length(time_vector);
        epoch = ([1:event_duration] + (event_duration) *(e-1));
        for ii = 1:size(DF,3)
            for jj = 1:size(DF,2)
                x = double(data1(pixels_select(jj,ii),time_vector ));
                xic = x(1:floor(window/2));
                xfc = x((tot_length - floor(window/2)+1):tot_length);
                F0temp = fastmedfilt1d(x, window-1, xic, xfc);
                temp1 = x'-F0temp;
                temp2 = temp1./F0temp;
    %             temp3 = smooth(temp2./F0temp,25,'sgolay',3); 
    %             xpoints = [1:size(temp3,1)];
    %             xq = [xpoints(1):f/f1:xpoints(end)];
    %             temp4 = interp1(xpoints,temp3,xq,'spline');
                DF(epoch,jj,ii) = temp1(frame_vector-time_vector(1)) *sign_peaks; 
                DFoF(epoch,jj,ii) = temp2(frame_vector-time_vector(1)) *sign_peaks; 
    %             sDFoF(:,jj,ii)=temp3; int_dFoF(:,jj,ii) = temp4';
            end
        end
    end