This is a long answer, but I don't have time to write a shorter one.
My reading of your calc_snr() code is that it assumes the signal power in the ECG is in the frequency band 0.67 to 150 Hz. These frequency values were no doubt chosen because they match the bandpass cutoff frequencies recommended by the 2007 Guidelines for ECG Recording, published by the Amer Heart Assoc and other cardiology societies. (An upper imit of 250 Hz is recommended for ECGs of children. The previous guidline, of 1990, specified 0.05 Hz to 150 Hz. The change from 0.05 Hz to 0.67 Hz is coupled to the fact that the 2007 Guidelines specifies the use of a forward-backward, zero phase filter.) The 2007 Guidelines do not explicitly state that filtering should be digital or analog, but the forward-backward, zero-phase filter guideline for the 0.67 Hz edge implies a digital filter.
The implementation of the 0.67 to 150 Hz band concept, in calc_snr (), has flaws.
Function calc_snr includes
N = fft(n);
spect = 2*abs(N)/length(N);
The factor of 2 on the right indicates the author intends spect to be the normalized one-sided spectrum. The inclusion of indices out to "end" when computing the noise power (PS), below, also suggest the author intends spect to be a one-sided spectrum. But it is not. To get the one-sided spectrum, spect should only include the first (half+1) of the elements of N. So that looks like an error to me.
calc_snr computes the signal power as follows:
Psignal=sum(PS(1+round(0.67/df):1+round(149/df)));
Since df is the frequency spacing of the DFT, in Hz, and bin 1 corresponds to frequency=0, the formula above computes signal power as the total power from 0.67 Hz to 149 Hz, inclusive. The use of round() means the closest frequencies to 0.67 and 149 Hz will be included. Note that this definition of signal power includes the line frequency and its 1st harmonic: 50 and 100 Hz, or 60 and 120 Hz. That will lead to overestimation of signal power, if there is line frequency noise.
calc_snr() computes the noise power as follows:
Pnoise = sum(PS(1:round(0.67/df)))/round(0.67/df)
+ sum(PS(1+round(150/df):end))/(round(500/df)-(1+round(150/df)))
+ PS(round(50/df));
The numerator of the first term is sum(PS(1:round(0.67/df))), which is the total power from 0 up to but not including 0.67 Hz.
The numerator of the second term is sum(PS(1+round(150/df):end)), which is the total power from 150 to the upper end of the spectrum (if PS is a one-sided spectrum - but it's not - see above).
The final term is PS(round(50/df)). This is the power at the bin closest to 50 Hz. The inclusion of this term is probably an attempt to include power at the (European) line frequency as "noise" power. But if that's the case, then it would have made sense to subtract that term from the "signal power" - which was not done. Also, if you are in North America, the line frequency is 60 Hz, not 50. Also, when there is line frequency contamination, it can be large at at harmonics of the line frequency, as well as at the line frequency. It happened to me in one recording. Due to a grounding issue, there was line noise in the ECG, which I did not notice until later. In the power spectrum, it appeared significantly at 60, 120, 180, 240, and 300 Hz.
Note that the noise power term from 0 to 0.67 Hz is divided by the number of frequency bins in that range:
Pnoise = sum(PS(1:round(0.67/df)))/round(0.67/df) + ...
round(0.67/df) is the number of frequency bins from 0 up to, but not including, 0.67 Hz. It is strange to divide the power in the range of frequencies by the number of bins. It was not done for the signal power, from 0.67 to 149 Hz, so why do it for the noise power?
Likewise, the power from 150 Hz to the end is divided by the number of frequency bins from 150 Hz to 500 Hz:
Pnoise = ... + sum(PS(1+round(150/df):end))/(round(500/df)-(1+round(150/df))) + ...
This suggests the sampling rate is probably 1 kHz, and therefore the top frequency in the spectrum is 500 Hz. But Fs is passed to calc_snr(), so it would make sense to use number of bins to the end, or (equivalently) to Fs/2, instead of hard coding it for 500 Hz. In any case, why divide the power in the 150-500 Hz band by the number of frequency bins in the band? It wasn't done for the signal power.
It is possible that the division by the number of bins is part of an attempt to estimate the "noise per bin". If this were the idea, then the noise per bin should be estimated corrctly (which it is not) and should be multipled by the total number of bins, from 0 to the Nyquist frequency, to estimate total noise power, if one assumes that the noise is white across all frequencies in the record, not just the bins below 0.67 Hz and above 150 Hz. But the noise per bin is not estimated correctly, and it is not multiplied by the total number of bins. So I'm not sure what the original author was trying to do.
I have explained the underlying theory of calc_snr(). I have also noted multiple inconsistencies in the calculation. Even if the calculation were done "right", I am not convinced that I like it. How do you deal with the fact that the "signal" from 0.67 to 150 Hz may itself include noise, including line frequency (50 or 60 Hz) noise? One could assume the noise is white, and esitmate the noise power per Hz from the non-noise band, but the power outside 0.67 to 150 Hz may have been suppressed by a digital bandpass filter built in to the ECG machine. That depends on th recording equipment and settings. And I am not convinced that the assumption of whiteness is justified. Noise at the line frequency is certainly not white.
I could fix the SNR calculation in calc_snr(), by dealing with the one-sided versus two-sided problem, and by esitmating noise per Hz in a more correct way from the out-of-band spectrum. I could post the result in the File Exchange. But I hesitate to do so, because I have quesitons about the theoretical justification for such a defintiion of SNR in the ECG, and I do not want to post a peice of code that could reify an unjustified concept.
Matlab's snr() is based on the assumption that the signal is produced by a weeakly nonlinear system driven by a pure sinusoid. The output of such a system will have output signal power at the fundamental and its first few harmonics. This could be a good model for testing an audio amplifier, for example. But it is not a good model for the ECG. Therefore I would not use Matlb's snr() for the ECG. Details: Matlab's snr() finds the frequency (fmax) in the power spectrum where power is maximal. This is the assumed fundamental frequency of an input sinusoid. The power at adjacent frequency bin with power monotonically declining from the peak are also included in the "signal" power. Signal power also includes power at the next five harmonics of the fundamental (i.e. power at fmax, 2*fmax, ..., 6*fmax). Noise power is power estimated from power at the non-signal frequencies. The Help does not make it clear exactly how the noise power is estimated.
