flatironinstitute / ironclust

Spike sorting software being developed at Flatiron Institute, based on JRCLUST (Janelia Rocket Cluster)
Apache License 2.0
28 stars 7 forks source link

What exactly is SNR or Sclu.vrSnr_clu? #140 #48

Open kouichi-c-nakamura opened 4 years ago

kouichi-c-nakamura commented 4 years ago

I wonder what SNR or Sclu.vrSnr_clu in the *_jrc.mat output file means.

Here it says:

"SNR: |Vp/Vrms|; Vp: negative peak amplitude of the peak site; Vrms: SD of the Gaussian noise (estimated from MAD)')"

Or here it says:

Signal to noise ratio at the peak site (SNR = Vpeak / Vrms)

RMS?

S.S_clu.vrVrms_site is likely the Vrms mentioned above. Does RMS stand for "Root-mean-square"?

I guess MAD here stands for "median absolute deviation" (not "mean absolute deviation"). Am I right?

According to this MAD*1.4286 can be approximation of SD.

What is the relationship between RMS and MAD?

Vp?

Then what is the Vp above?

The description "negative peak amplitude of the peak site" implies Vp can be S_clu.vrVmin_clu or S_clu.vrVmin_uv_clu.

However, neither S_clu.vrVmin_clu./S.S_clu.vrVrms_site or S_clu.vrVmin_uv_clu./S.S_clu.vrVrms_site gave me the same values as Sclu.vrSnr_clu.

Also I can't find the lines in the source code where this SNR is computed.

So many things are unclear to me. Can somebody help me clarifying things, please?

jamesjun commented 4 years ago

It measures the absolute negative peak amplitude at the peak channel divided by the RMS noise estimated using MAD (see Quiroga et al 2004). I don’t use SD directly because it’s influenced by the firing rate and spike contamination. I will document this more thoroughly since I am back from the Cosyne conference.

-James

On Mar 6, 2020, at 9:30 AM, Kouichi C. Nakamura notifications@github.com wrote:

 I wonder what SNR or Sclu.vrSnr_clu in the *_jrc.mat output file means.

Here it says:

"SNR: |Vp/Vrms|; Vp: negative peak amplitude of the peak site; Vrms: SD of the Gaussian noise (estimated from MAD)')"

Or here it says:

Signal to noise ratio at the peak site (SNR = Vpeak / Vrms)

RMS?

S.S_clu.vrVrms_site is likely the Vrms mentioned above. Does RMS stand for "Root-mean-square"?

I guess MAD here stands for "median absolute deviation" (not "mean absolute deviation"). Am I right?

According to this MAD*1.4286 can be approximation of SD.

What is the relationship between RMS and MAD?

Vp?

Then what is the Vp above?

The description "negative peak amplitude of the peak site" implies Vp can be S_clu.vrVmin_clu or S_clu.vrVmin_uv_clu.

However, neither S_clu.vrVmin_clu./S.S_clu.vrVrms_site or S_clu.vrVmin_uv_clu./S.S_clu.vrVrms_site gave me the same values as Sclu.vrSnr_clu.

Also I can't find the lines in the source code where this SNR is computed.

So many things are unclear to me. Can somebody help me clarifying things, please?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.

kouichi-c-nakamura commented 4 years ago

Thanks, James! I found the relevant code snippet.

%--------------------------------------------------------------------------
function S_clu = S_clu_quality_(S_clu, P, viClu_update)
% 7/5/17 JJJ: Added isolation distance, L-ratio, 
% TODO: update when deleting, splitting, or merging
if nargin<3, viClu_update = []; end
t1 = tic;
fprintf('Calculating cluster quality...\n');

[vrVpp_clu, vrVmin_clu, vrVpp_uv_clu, vrVmin_uv_clu] = deal([]);
try
    mrVmin_clu = squeeze_(min(S_clu.tmrWav_clu,[],1),1);
    mrVmax_clu = squeeze_(max(S_clu.tmrWav_clu,[],1),1);
    mrVmin_uv_clu = squeeze_(min(S_clu.tmrWav_raw_clu,[],1),1);
    mrVmax_uv_clu = squeeze_(max(S_clu.tmrWav_raw_clu,[],1),1);
    vrVpp_clu = mr2vr_sub2ind_(mrVmax_clu-mrVmin_clu, S_clu.viSite_clu, 1:S_clu.nClu);
    vrVmin_clu = mr2vr_sub2ind_(mrVmin_clu, S_clu.viSite_clu, 1:S_clu.nClu);
    vrVpp_uv_clu = mr2vr_sub2ind_(mrVmax_uv_clu-mrVmin_uv_clu, S_clu.viSite_clu, 1:S_clu.nClu);
    vrVmin_uv_clu = mr2vr_sub2ind_(mrVmin_uv_clu, S_clu.viSite_clu, 1:S_clu.nClu);
catch
    disp('S_clu_quality_: error');
end
try
    S0 = get(0, 'UserData');
    if isempty(S0), S0 = load0_(subsFileExt_(P.vcFile_prm, '_jrc.mat')); end
    vrVrms_site = bit2uV_(single(S0.vrThresh_site(:)) / S0.P.qqFactor, P);
%     vrSnr_clu = vrVpp_clu ./ vrVrms_site(viSite_clu);
    vrSnr_clu = abs(vrVmin_clu) ./ vrVrms_site(S_clu.viSite_clu);
    vnSite_clu = sum(bsxfun(@lt, mrVmin_clu, -vrVrms_site * S0.P.qqFactor),1)';
catch
    [vrVrms_site, vrSnr_clu, vnSite_clu] = deal([]);
    disp('no Sevt in memory.');
end
try
    [vrIsoDist_clu, vrLRatio_clu, vrIsiRatio_clu] = S_clu_quality2_(S_clu, P, viClu_update);
catch
    [vrIsoDist_clu, vrLRatio_clu, vrIsiRatio_clu] = deal(nan(size(vrVpp_clu)));
end

S_clu = struct_add_(S_clu, vrVpp_clu, vrSnr_clu, vrVrms_site, vnSite_clu, ...
    vrIsoDist_clu, vrLRatio_clu, vrIsiRatio_clu, vrVpp_uv_clu, vrVmin_uv_clu);
fprintf('\ttook %0.1fs\n', toc(t1));
end %func

https://raw.githubusercontent.com/flatironinstitute/ironclust/master/matlab/irc.m

In essense:

S_clu.vrSnr_clu == abs(S_clu.vrVmin_clu)./S_clu.vrVrms_site(S_clu.viSite_clu)

I confirmed the equation above is true.