Determining Signal to Noise Ratio

Hi everyone,

I would like to determine the Signal to Noise Ration for an electrical signal, which was already amplified. I have multiple input channels, which experience noise and sometimes a real signal, which is easily visible with the naked eye. My Input arguments are x (Voltage) and y (time). My goal is to determine the SNR accurately for all channels, to compare their SNR and performance.

I have tried using the snr-function in MATLAB in two ways:

1. r = snr(x): I'm not sure if this is correct for my application since this function is meant for a "real-valued sinusoidal input signal x", which I´m not certain my signal is. Also, channels which just experience noise somehow have a higher SNR than channels with visible peaks when I use this formula.
2. r = snr(xi,y): I manually estimated the noise y (not to be confused with time, I just copied the definition from MATLAB), for the upper example about 0.004. The problem is the exact same as before, channels with just noise have a higher SNR than others. Also for some reason it does not matter what I input as noise, the result from the snr-function always stays the same.

One last question: Can a channel which does not record a valuable signal, just noise, even posses a SNR? Maybe thats part of my mistake.

I know I´m doing something wrong or I´m missing something, any help would be greatly appreciated!

[Paul]

I don't fully understand snr. Example straight from that doc page:

rng(100);
fs = 48e3;
t = 0:1/fs:1-1/fs;
A = 1.0;
powfund = A^2/2;
a = 0.4;
powharm = a^2/2;
s = 0.1;
varnoise = s^2;
x = A*cos(2*pi*1000*t) + ...
    a*sin(2*pi*2000*t) + s*randn(size(t));
defSNR = 10*log10(powfund/varnoise)

defSNR = 16.9897

figure
snr(x)

ans = 17.0130

copyobj(gca,figure)

xlim([0 50])

If the sampling frequency is 48 kHz, how can the frequency axis go out to 500 megaHz?

Why is the fundamental shown at 21 mHz, as opposed to 1 kHz?

Why is the noise power shown at around -60 dB? Shouldn't it be around

10*log10(varnoise)

ans = -20

which would then result in the expected (and shown) snr of 17 dB, because the power in the fundamental is

10*log10(powfund)

ans = -3.0103

Can anyone clarify ...

Answer 1

[Star Strider]

Calculating SNR can be difficult if the underlying signal is unknown. An extended discussion is in the Analyzing Harmonic Distortion documentation section.

A sort of ‘hack’ I sometimes use involves first computing the fft of the signal to see what its frequency content is. After that, since broadband noise is not completely eliminated by frequency-selective filters, I use the Savitzky-Golay filter (the sgolayfilt function) to eliminate as much of the broadband noise as I can. (It operates as a sort of FIR filter with multiple attenuation frequencies.) I then subtract that result from the original signal, calculate the noise power, and the ‘filtered’ signal poweor, and be done with it.

With respect to the ‘sinusoidal signal’ requirement, as a general rule, any periodic signal can be expressed as a sum or sinusoids of varying frequency and phase (this underlies the assumptions of the Fourier transform), so that may be enough.

I usually work with physiologic signals that are nonlinear and nonstationary, so always challenging. If there is some sort of magickal procedure that will always produce a robust and correct result for SNR, I have thus far not discovered it.

I am not certain that this answers your question, so if it does not, I will delete it.

[M]

Thanks for your reply!

I´ll definetly read the discussion you linked to help my understanding of the topic and thanks for clarifying that there is no perfect result for SNR.

I also tried to implement your suggestion and ran into the problem that my Savitzky-Golay filtered signal doesn't really look filtered, I used order = 3 and framelength = 501.

C8fft = fft(C8.y);

sgf = sgolayfilt(C8.y,3,501);

plot(sgf,C8.y)

noise = abs(C8.y - sgf);

noisePower = mean(noise);

SNR = mean(C8.y/noisePower)

This is the function I used to calculate SNR based on your recommendation, did I do something wrong or is there room for improvement?

[Star Strider]

My pleasure!

I don’t beleive that you did anything ‘wrong’, just that there is no perfect way to do it.

First, examine the fft of your signal to see if there is actually a signal in all the noise. (I have my own function to do that, and also use the pspectrum function.)

You may have to increase the ‘framelen’ even further. I usually begin using an order 3 polynomial (as you’re doing) with:

L = size(signal,1); % Use The Correct Dimension For Your Signals

framelen = fix(L/10) + rem(fix(L/10),2);

and see what that result is. (The MATLAB signal processing functions are column-major, so getting the column length is important. I use the numel function for vectors.) I then reduce it experimentally until I get the result I want (that of course varies, depending on what the signal itself is). I don’t have your signal (at least I didn’t see it uploaded here), so I can’t experiment with it.

.

[M]

I implemented the changes and it works very well by now, my SNR is as between 5.5-6.5 on average, and behaves the way I hoped it would across channels.

Thank you for your help Star Strider!

[Star Strider]

As always, my pleasure!

Answer 2

[Image Analyst]

You need to smooth the data somehow, like with movmean, movmedian, or sgolayfilt. Then once you have the "noise free signal, you can subtract and get some stats on it.

smoothSignal = movmean(signal, 9);
noiseOnly = signal - smoothSignal; % A vector
snrSignal = signal / abs(noiseOnly);
meanNoise = mean(noiseOnly); % Or use rms function.
meanSNR = mean(snrSignal);

or something similar.

[M]

Thanks for your reply!

I tried implementing your suggestion, the only problem I ran into: snrSignal creates a 50002x50002 array, since I have 50002 data points for each channel. This also means that meanSNR returns a vector instead of an absolute number. Did I do something wrong?

Also I just tried using snr with the calculated noiseOnly:

smoothSignal = movmean(C8.y, 501);

noiseOnly = C8.y - smoothSignal; % A vector

snr8 = snr(C8.y,noiseOnly)

With still faces the same problems of channels with only noise having a high SNR.

Is there anything I should do differently?

[Image Analyst]

Just split it up into each channel individually and process each individually. Attach your data in a .mat file if it's less than 5 MB, otherwise clip out a snippet of it to make it smaller.

[M]

I'm already calculating each channel individually, each channel has 50002 data points per event. But I can handle that by turning off the array size limiter. But my main question was how to fix the meanSNR output, since snrSignal is a matrix, meanSNR outputs a 50002 colunm long vector, instead of a single value for SNR of the channel.

[Voss]

snrSignal = signal / abs(noiseOnly);

should be

snrSignal = signal ./ abs(noiseOnly);

[Image Analyst]

Make sure both signal and smoothSignal are both either row vectors or column vectors. They're probably one of each and automatic expansion caused it to create a huge 2-D matrix.

[M]

Thanks, I managed to get rid of the problem, but my SNR was still all over the place using this method, so I`ll stick to the other one for now.

Nonetheless thanks for your help Image Analyst!

[Image Analyst]

OK, no problem. Though the "other method" was also my method. There are several ways to smooth the data and Star and I both suggested sgolayfilt to smooth the data to get an estimate of the "noise-free" data, and then subtracting that from the signal to get an estimate of only the noise. Then compute the SNR from the ratio of the noise-free signal to the noise-only signal.