SNR calculation via FFT-based crosscorrelation

Summary

A Matlab / Octave example to obtain the signal-to-noise ratio of a signal against a reference signal (correlation)

It works on both real- and complex valued data, for example IQ data streams.

Crosscorrelation

The method uses the normalized crosscorrelation coefficient

between two signals

and

tells what part of

is explained by knowing

, and vice versa.
For example:

implies that and are equal, except for a (possibly complex) scaling factor.
When , and have “half” of their waveform in common.

Samples and power

Samples could represent a physical quantity such as electrical voltage or current. Power is proportional to the square of the amplitude.
The same applies, even if the samples do not represent an electric signal.
For example,

explains one quarter (

) of one signal's power knowing the other signal.

Signal-to-noise ratio

The common part

is fully correlated. The remainder

is fully uncorrelated, and considered as noise.
The signal-to-noise-ratio follows then as:

$SNR_{\mathrm{lin}}=\frac{|c|^2}{1-|c|^2}$

Since it is a ratio of powers, it is converted to dB using

$SNR_{\mathrm{dB}}=10\log_{10}SNR_{\mathrm{lin}}$

Accuracy

Lower SNRs will require a longer signal for a given accuracy, since part of the noise is randomly correlated with the signal.

Example code: Octave (tested) or Matlab

download complete example (right-click, save as)

The program includes the following steps:

loop through a number of signal-to-noise ratios
generate a signal
add noise according to the nominal SNR
find crosscorrelation coefficients for all possible time shifts using cyclic convolution (FFT-based)
find the peak, which is
calculate the SNR from the signal
Plot the resulting SNR versus the nominal SNR

% SNR detection via FFT / correlation
% Markus Nentwig 2007
% this program is in the public domain, provided without any warranty

close all; clear all;
N=65536;

% transmitted signal
a=randn(1, N);

x=[];
y=[];
for SNR_in_dB=-30:5:60
  vRatio=power(10, -SNR_in_dB/20);

  % construct time shifted, noisy "received" signal
  b=circshift(a, [0, -1234])+vRatio*randn(1, N);

  % scale to arbitrary amplitude...
  a=a*12345;
  b=b*23456;

  % ... and normalize
  ea=sum(a .* conj(a)) / N;
  eb=sum(b .* conj(b)) / N;
  a=a/sqrt(ea);
  b=b/sqrt(eb);

  % crosscorrelation
  corr=ifft(fft(a) .* conj(fft(b)))/N;

  % find peak (should be 1234+1)
  maxBin=find(abs(corr)==abs(max(corr)));
  c=abs(corr(maxBin));

  % c is the ratio of "voltage" in a explained by b and vice versa
  % c^2 is the ratio of power => signal
  % 1-c^2 is the remaining part unexplained by signal => noise
  % note, low SNRs require a longer signal, because part of the noise will randomly correlate with the reference signal.

  SNR_out=c^2/(1-c^2);
  SNR_out_dB=10*log10(SNR_out);
  x=[x, SNR_in_dB];
  y=[y, SNR_out_dB];
end

plot(x, y, '+-');
xlabel('SNR in/dB'); ylabel('SNR out/dB');

Result

Output from example program

Timing

A delay mismatch between

and

results in a reduction of reported SNR.
The result isn't particularly sensitive, as long as the signal is narrowband compared to the sampling rate (because the center peak of the signal's autocorrelation function is wide).
For signals using a large percentage of the bandwidth, for example oversampling by 2, even small timing errors may have a significant impact on SNR.
In such a case, it is possible to raise the sampling rate of signal and reference through FFT-interpolation, giving better timing resolution.

Limiting the bandwidth

A meaningful SNR definition requires a filter.
Noise that is filtered out by the (ideal) receiver should be excluded.

Crosscorrelation-based SNR without FFT

The above example uses FFT to calculate the crosscorrelation. As an additional benefit it finds the optimum timing to align signal and reference, with a resolution of one sample.
But frequently,

and

are already time-aligned, and no FFT is needed.
In this case, the crosscorrelation can be calculated directly on the normalized sample streams. The conjugate is only required, when samples are complex-valued.
download complete example

...
% crosscorrelation
c=sum(a .* conj(b) / N);
...

Comments? Questions?

Please send me a mail! mnentwig@elisanet.fi