VSD_SNR.m

function [ data_rec ] = VSD_SNR( data_rec, pixels, plotflag, size_image )
% This function computes the SNR of VSD data defined as the amplitude of
% the VSD signal divided by the RMS value of the baseline noise of the dFoF
% traces
%   Detailed explanation goes here:

traces = data_rec.dFoF;

% locate events
avgSignal= squeeze(mean(mean(data_rec.data,1),2));
nFrames = length(avgSignal);
window = round(1/data_rec.metadata.KineticCycleTime * 15); % 1500 at 100 Hz % median filter (fast implementation compiled in C)
minSBWnd=300;
minSBInterval=300;
xic = avgSignal(1:floor(window/2));
xfc = avgSignal(nFrames - floor(window/2)+1:nFrames);
avgF0 = fastmedfilt1d(avgSignal,window,xic,xfc);
avgdFoF=(avgSignal-avgF0)./avgF0;
DataPoints = [1:nFrames] .* (data_rec.metadata.KineticCycleTime*1000);
[BS,BP,BE,BI]=eventLocator(-avgdFoF,DataPoints,[1;1;1;nFrames],minSBWnd,minSBInterval,'stdThresh',10,'medianWindow',2000);

% event windows
for i = 1:length(BS)
    frames = find(DataPoints>= BS(i) & DataPoints<= BE(i));
    if i == 1
        eventframes = frames;
        events = traces(frames,:);
    else
        eventframes = [eventframes, frames];
        events = cat(1,events,traces(frames,:));
    end
end

amplitudes = prctile(events,95);

% determine values of corrected baseline and noise baseline
baseline = traces(setdiff(floor(window/2):nFrames-floor(window/2),eventframes),:);
noise_RMS = sqrt(mean(baseline.^2,1));

% compute SNR

SNR = amplitudes ./ noise_RMS ;

data_rec.SNR_rois = SNR;

% plotting

if plotflag
    figure, imagesc(reshape(SNR,size_image(1), size_image(2))), colorbar, title('map of pixel-wise SNR')
    
    [snr_sorted order] = sort(SNR);
    snr_3 = order(find(snr_sorted>=3 & snr_sorted<4));
    snr_4 = order(find(snr_sorted>=4 & snr_sorted<5));
    snr_5 = order(find(snr_sorted>=5 &snr_sorted<6));
    snr_7 = order(find(snr_sorted>=7 &snr_sorted<8));
    snr_10 = order(find(snr_sorted>=10 &snr_sorted<11));
    snr_11plus =  order(find(snr_sorted>=11));

    figure, 
    a = randi(length(snr_5),1,2);
    b = randi(length(snr_7),1,2);
    c = randi(length(snr_10),1,2);
    d = randi(length(snr_11plus),1,2);
    e = randi(length(snr_3),1,2);
    f = randi(length(snr_4),1,2);
    plot(traces(:, snr_5(a) ),'k','LineWidth',1.6), hold on
    plot(traces(:,snr_7(b)),'b','LineWidth',1.6),
    plot(traces(:,snr_10(c)),'g','LineWidth',1.6), 
    plot(traces(:,snr_11plus(d)),'r','LineWidth',1.6), 
    plot(traces(:, snr_3(e) ),'Color',[0.75 0.75 0.75],'LineWidth',1.6)
    plot(traces(:, snr_4(f) ),'c','LineWidth',1.6)
    title('example dFoF traces for different SNR levels')

    traces_int = data_rec.dFoF_interpolated;
    figure, 
    plot(traces_int(:, snr_5(a) ),'k','LineWidth',1.6), hold on
    plot(traces_int(:,snr_7(b)),'b','LineWidth',1.6),
    plot(traces_int(:,snr_10(c)),'g','LineWidth',1.6), 
    plot(traces_int(:,snr_11plus(d)),'r','LineWidth',1.6), 
    plot(traces_int(:, snr_3(e) ),'Color',[0.75 0.75 0.75],'LineWidth',1.6)
    plot(traces_int(:, snr_4(f) ),'c','LineWidth',1.6)
    title('example upsampled dFoF traces for different SNR levels')

    figure,
    plot(mean(traces(:, snr_4),2),'c','LineWidth',1.8), hold on
    plot(mean(traces(:, snr_5),2),'k','LineWidth',1.8),
    plot(mean(traces(:, snr_7),2),'b','LineWidth',1.8),
    plot(mean(traces(:, snr_10),2),'g','LineWidth',1.8), 
    plot(mean(traces(:, snr_11plus),2),'r','LineWidth',1.8), 
    plot(mean(traces(:, snr_3),2),'Color',[0.75 0.75 0.75],'LineWidth',1.8)
    title('mean trace different SNR levels')
end

end