MATLAB SNR QUESTION: need help using awgn() function to make noisy signal and then assessing that signal SNR with custom code

Question

I am using awgn() to make a signal 'x' with a given amount of SNR, in this case 10 [line 19].

I am trying to accurately assess the SNR of the signal it creates with my own code, but I am calculating way more signal gain from the noise floor than there should be [line 44].

f = 1000;
t = 0:.0001:1;
x = cos(2*pi*f*t);
len = length(x);
fs = (0:len-1)/t;


plot(t,x)
xlim([0 .01])
title('Signal "x"')

%fft of x 
X = fftshift(abs(fft(x)));                              
TF = (-len/2:len/2-1)*fs/len;
plot(TF,X)
title('FFT of "x"')

%add noise
y = awgn(x,10);                                          

plot(t,y)
xlim([0 .01])
title('Signal "y"')

%fft of y
Y = 10*log(abs(fftshift(fft(y)/length(y))));                
plot(TF,Y)
grid
title('FFT of "y"')
xlim([500 1500])

%Finding power of signal
MinH = max(Y) -10;
[pks,locs] = findpeaks(Y,"MinPeakHeight",MinH)
Peaks = tf(locs);
fcSig1 = 0.5*(Peaks(1,1) + Peaks(1,2))
S = (pks(1,1)+pks(1,2))/2

%Finding power of white noise
[negpks,neglocs] = findpeaks(-(abs(Y)));
avgN = sum(negpks)/length(negpks);
N = avgN

SNR = S-N %in dB

I tried sampling the average power of the Signal, the average power of the Noise floor, turning those powers into dB units, and the subtracting the noise dBm from the signal dBm.

My results are here: In case you can't open: [S = -7.1313, N = -55.5702, SNR = 48.4389 (38dB higher than my input of 10)]

My thoughts are that I am sampling my noise floor inaccurately, or I needed to account for bandwidth, but I am not sure how to do that.

Any suggestions or comments are greatly appreciated!

Answer 1

This may not be the answer you are looking for, but Matlab does have a function called snr. The definition is here: https://au.mathworks.com/help/signal/ref/snr.html. The way it is computed is using the root sum square of the signal so it is computed using the time domain rather than the frequency domain. I'm not sure how to derive the root sum square from a magnitude spectrum.

I re-wrote the signal generation code a bit and kept the plotting in to show that it is roughly the same signal as with awgn (which I don't have).

%%
f = 100;
fs = 1000;
t = 0:1/fs:1-1/fs;
x = cos(2*pi*f*t);

% fft of x 
X = 10*log(abs(fft(x)/length(x)));                              
TF = linspace(0,fs,length(x));

%add noise
%y = awgn(x,10);     
% I don't have 2023, so no awgn, but 10dB is about 3*amplitude
n = rand(size(x));
n = n*2;
n = n-1;
n = n / 3;
y = x+n;

%fft of y
figure;
subplot(2,1,1)
plot(TF,X)
grid
title('magnitude spectrum of "x"')
xlim([0 fs/2])
subplot(2,1,2)
Y = 10*log(abs(fft(y)/length(y)));                
plot(TF,Y)
grid
title('magnitude spectrum of "y"')
xlim([0 fs/2])

Using the time domain you can do this:

s = snr(x,y-x) % equivalent to snr(x,n)
% or:
r = mag2db(rssq(x)/rssq(y-x)) 

Outputs:

s =

   11.3894


r =

   11.3894

So, not exactly 10, but I just used a rough estimate (6dB ~= double amplitude, so *3 is about 10dB) to figure the noise scale. Also, the output will vary slightly depending on the actual content of the noise.

I tidied things up a bit. E.g. I define the sampling rate first and use that to determine vector lengths etc. And you don't really need fftshift for a 1D DFT. That way, the frequency bins are defined much more simply.