SNR calculation of FFT Spectrum

[Chiel]

Hello,

I am working on a project where I need to calculate the SNR of a signal using a FFT plot.

This plot is of a sampled sine wave from a vector signal generator, sampled by an ADC converter. A plot of the fft is shown in the next figure:

http://img269.imageshack.us/img269/9894/screenshot2id.png

I used a flattopwindow when calculating the FFT.

The sine wave has a magnitude of -35 dBm. Normally the SNR is calculated by the signal power devided by the mean of the noisepower. Because of spectral leakage, this formula does not hold. How can i calculate the SNR from this FFT?


Here is the code that I used to make the FFT

% Load the data into matlab and remove datafile from disk
data=read_complex_binary('data.dat');
system('rm data.dat');

% Create Flat Top weighted window and multiply data with the window
window=flattopwin(length(data),'periodic');
window=window/(mean(window));
data=window.*data;

% Make a fourrier transform of the data and determine the value of the peak
fft_data = abs(fft(data)/(samples-1));
fft_power = ((fft_data.*fft_data)/50)*1000; % Convert voltage to mW
fft_power = 10*log10(fft_power);

plot(fft_power);

Thanks in advance!

[dbd]

On Nov 9, 4:03 am, "Chiel " <c.hakkenb...@student.utwente.nl> wrote:
> Hello,
> ...

> Normally the SNR is calculated by the signal power devided by the mean of the noisepower. Because of spectral leakage, this formula does not hold. How can i calculate the SNR from this FFT?

> ...


Take a look at Fred Harris IEEE Proceedings paper at:

http://www.utdallas.edu/~cpb021000/EE%204361/Great%20DSP%20Papers/Harris%20on%20Windows.pdf

In particular, section IV subsection B. Processing Gain, for
correction due to the window.

Now, what are you going to do about the power in the peaks that may be
at DC and the signal's harmonics?

Dale B. Dalrymple

[Paul Nelson]

dbd <d...@ieee.org> wrote in message <bb89ed8b-4873-4510...@u10g2000prl.googlegroups.com>...

Will do, thanks for the link.

I don't mind calculating the SNR in the time-domain, either, if that is eaiser. I just want to make sure that the signal actually has the SNR defined by awgn (and hopefully understand why).

[Paul Nelson]

OK, so I've read the posted IEEE paper, as well as this reference (page 6) http://www.analog.com/static/imported-files/tutorials/MT-001.pdf

As I understand it, the FFT creates a processing gain, thereby increasing the signal to noise ratio as a function of the number of points in the FFT. I understand this concept now, but am still having a bit of difficulty in implementation.

I calculate the 'correct' SNR in frequency and time (see variables SNR_freq and SNR_time),
but when I attempt to calculate directly from the FFT output, the results are incorrect.

It must be a small error for SNR_freq1, because it is within 3 dB. Perhaps I need to do something like scale the amplitude of the maximum by 3dB? Do I need to do anything special since I'm dealing from -Fs/2 to Fs/2 instead of 0 to Fs/2?

Thanks!

BTW, here is my new code:

% SNR test
clc; clear all; close all;

% Units
MHz = 1e6;

% Variables
Fs = 10*MHz; % Sample rate
Ts = 1/Fs; % Sample period
L = 4096; % Sample record length
m = 0:L-1; % Sample vector
f_if = 1*MHz;

omega_if = 2*pi*f_if; % Analog IF frequency [rad/s]
W_if = 2*pi*f_if/Fs; % Digital IF frequency [rad/s]
theta_n = deg2rad(30); % Phase shift [radians]
SNR = 40; % Signal to noise ratio of message waveform x[m] [dB]

g0 = cos(W_if.*m - theta_n); % Generate ideal signal and phase shift
g1 = awgn(g0,SNR,'measured'); % Add white Gaussian noise to signal

% SNR estimate of signal to verify awgn is properly used
% Signal power (Integrate over samples of original signal, divide by record length)
p_signal = sum(abs(g0).^2,2)./length(g0)
% Signal noise (Calculate variance of noise)
p_noise = std(g1 - g0,0,2).^2
% SNR calculation in time
SNR_time = 10*log10(p_signal./p_noise)


% Perform FFT of g0
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
%NFFT = 1024;
G0 = fft(g0,NFFT,2)/L; % FFT
G0 = fftshift(G0,2); % Shift spectrum to [-Fs/2,Fs/2]
G0_dB = 20*log10(abs(G0)); % Convert to dB

% Perform FFT of g1
G1 = fft(g1,NFFT,2)/L; % FFT
G1 = fftshift(G1,2); % Shift spectrum to [-Fs/2,Fs/2]
G1_dB = 20*log10(abs(G1)); % Convert to dB

f = Fs/2*linspace(-1,1,NFFT); % Generate frequency list

% SNR estimate using spectrum of signal to verify awgn is properly used
% Signal power (Integrate over samples of original signal, divide by record length)
p_signal_f = sum(abs(G0).^2,2)./length(G0)
% Signal noise (Calculate variance of noise)
p_noise_f = std(G1 - G0,0,2).^2
p_noise_f1 = sum(abs(G1 - G0).^2,2)
% SNR calculation
SNR_freq = 10*log10(p_signal_f./p_noise_f)

% Attempt to calculate SNR from spectrum, subtracting processing gain
SNR_freq1 = max(G1_dB)- 10*log10(p_noise_f) - 10*log10(L/2)
SNR_freq2 = max(G1_dB)- 10*log10(p_noise_f1) - 10*log10(L/2)

% Plot time-domain signal
figure;
plot(m,g0,m,g1)
hold on;
%stem(m,g1,'r')
xlabel('Sample Number')
ylabel('Value')
x0 = Fs/(omega_if/(2*pi))*4;
xlim([0 x0])
legend('Signal','Signal + Noise')
hold off;

% Plot spectrum magnitude.
figure;
plot(f./MHz,G0_dB,f./MHz,G1_dB)
title('Spectrum Magnitude of g0(t) and g1(t)')

xlabel('Frequency (MHz)')
ylabel('Magnitude')
legend('Signal','Signal + White Noise')

20*log10(mean(abs(G1))) + 10*log10(L/2) + SNR
max(G1_dB) - (20*log10(mean(abs(G1))) + 10*log10(L/2))