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:
- Obtain the frequency spectrum of your time series data using the Fast Fourier Transform (FFT). MATLAB provides the “fft” function for this purpose.
- 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.
- Calculate the power of the target frequency component by squaring the magnitude of the complex value at the identified frequency bin.
- 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.
- 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!
