Determine snr for specific frequencies

I have noisy time series data with sinusodial target signals (500, 1000, 2000 Hz) recorded at 20-40kHz sampling rate using different bio amplifiers in a controlled setup. Calculating the snr, I try to evaluate the recording quality of the different devices.

Some recordings show massive undesired artefacts (environmental crosstalk) with neighbouring frequencies, which often get mistaken as the fundamentals when applying the snr matlab function. The resulting snr does obviously not reflect the snr of the target frequency, although there is a visible if minor peak in the frequency domain.

I would like to specify the snr of my target frequencies or specify some adequate fundamental encompassing the target frequencies, independet from the fundamental component, that is automatically determined by the matlab function snr.

Any ideas?

Can anyone help me to step by step reproduce the result of the snr matlab function?

Or how would you in detail calculate SNR from time series or PSD by hand?

What unit would you prefer for snr and for what reason?

Many thanks.

This is an example output of snr(signal) where signal is a response to a 2kHz sinusodial tone, amplitude is too low to get detected as the fundamental frequency, but I need to know SNR of the 2kHz component in this recording.

Answer 1

[Karan Singh]

Hi Florian,

From what I understand, the goal is to calculate the SNR (Signal-to-Noise Ratio) of the target frequencies to assess the recording quality. However, there is a challenge due to the presence of undesired artifacts caused by environmental crosstalk with neighbouring frequencies. These artifacts can be mistaken for target frequencies when applying the MATLAB “snr” function, resulting in inaccurate SNR values.

To calculate the SNR (Signal-to-Noise Ratio) of a specific frequency component in your time series data, you can follow these steps:

1. Obtain the frequency spectrum of your time series data using the Fast Fourier Transform (FFT). MATLAB provides the “fft” function for this purpose.
2. Identify the frequency bin corresponding to your target frequency (2kHz in your case) in the frequency spectrum. You can calculate the bin index using the formula “bin_index = round(target_frequency * N / fs)”, where “N” is the number of samples in your time series data and “fs” is the sampling rate.
3. Calculate the power of the target frequency component by squaring the magnitude of the complex value at the identified frequency bin.
4. Calculate the total noise power by summing the power of all the other frequency components except the target component. You can exclude a range of bins around the target frequency to exclude any neighboring frequencies that might be causing artifacts.
5. Finally, calculate the SNR using the formula “SNR = 10 * log10(target_power / noise_power)”, where “target_power” is the power of the target frequency component and “noise_power” is the total noise power.

Here is an example code snippet that demonstrates the steps described above:

% Example time series data with a 2kHz sinusoidal target component 
fs = 40000; % Sampling rate in Hz 
duration = 1; % Duration of the time series data in seconds 
time = 0:(1/fs):duration-(1/fs); 
target_frequency = 2000; % Target frequency in Hz 
target_amplitude = 0.1; % Amplitude of the target component 

% Generate the time series data with noise 
noise_amplitude = 0.05; % Amplitude of the noise 
noise = noise_amplitude * randn(size(time)); 
signal = target_amplitude * sin(2*pi*target_frequency*time) + noise; 

% Calculate the SNR of the target frequency component 
N = length(signal); % Number of samples 
frequencies = (0:N-1) * (fs/N); % Frequency axis for FFT 
spectrum = fft(signal); % Compute the FFT 
power_spectrum = abs(spectrum).^2 / N; % Compute the power spectrum 

% Find the frequency bin index corresponding to the target frequency 
target_bin_index = round(target_frequency * N / fs); 

% Calculate the power of the target frequency component 
target_power = power_spectrum(target_bin_index); 

% Exclude neighboring frequency bins to calculate the noise power 
excluded_bins = [target_bin_index-10 : target_bin_index+10]; % Exclude 10 bins around the target frequency 
excluded_bins(excluded_bins < 1 | excluded_bins > N) = []; % Remove out-of-range bins 
noise_power = sum(power_spectrum(setdiff(1:N, excluded_bins))); 

% Calculate the SNR in dB 
SNR = 10 * log10(target_power / noise_power); 

disp(['SNR of the target frequency component: ' num2str(SNR) ' dB']); 

In this example, the code generates a time series data with a 2kHz sinusoidal target component and additive Gaussian noise. It then calculates the SNR of the target frequency component using the steps outlined above.

The SNR is typically expressed in decibels (dB) as it provides a logarithmic scale that is more suitable for comparing signal power levels. A higher SNR value indicates a better signal quality relative to the noise level.

Attached below are some documentation links that you may find helpful:

• Signal-to-noise ratio - MATLAB snr (mathworks.com)
• Fast Fourier transform - MATLAB fft (mathworks.com)
• Signal Processing Toolbox Documentation (mathworks.com) (For Power Spectral Density (PSD))

Hope this helps!