1% estsnr.m 2% David Rowe May 2017 3% 4% estimate SNR of a sinewave in noise 5 6function snr_dB = estsnr(x, Fs=8000, Nbw = 3000) 7 8 [nr nc] = size(x); 9 if nr == 1 10 x = x'; 11 end 12 13 % find peak in +ve side of spectrum, ignoring DC 14 15 L = length(x); 16 X = abs(fft(x)); 17 st = floor(0.05*L); en = floor(0.45*L); 18 [A mx_ind]= max(X(st:en)); 19 mx_ind += st; 20 21 % signal energy might be spread by doppler, so sum energy 22 % in frequencies +/- 1% 23 24 s_st = floor(mx_ind*0.99); s_en = floor(mx_ind*1.01); 25 S = sum(X(s_st:s_en).^2); 26 27 % real signal, so -ve power is the same 28 29 S = 2*S; 30 SdB = 10*log10(S); 31 32 printf("Signal Power S: %3.2f dB\n", SdB); 33 34 % locate a band of noise next to it and find power in band 35 36 st = floor(mx_ind+0.05*(L/2)); 37 en = st + floor(0.1*(L/2)); 38 39 N = sum(X(st:en).^2); 40 41 % scale this to obtain total noise power across total bandwidth 42 43 N *= L/(en-st); 44 NdB = 10*log10(N); 45 printf("Noise Power N: %3.2f dB\n", NdB); 46 47 % scale noise to designed noise bandwidth /2 fudge factor as its a 48 % real signal, wish I had a better way to explain that! 49 50 NodB = NdB - 10*log10(Fs/2); 51 NscaleddB = NodB + 10*log10(Nbw); 52 snr_dB = SdB - NscaleddB; 53 54 figure(1); clf; 55 plot(20*log10(X(1:L/2)),'b'); 56 hold on; 57 plot([s_st s_en], [NdB NdB]- 10*log10(L), 'r'); 58 plot([st en], [NdB NdB]- 10*log10(L), 'r'); 59 hold off; 60 top = 10*ceil(SdB/10); 61 bot = NodB - 20; 62 axis([1 L/2 bot top]); 63 grid 64 grid("minor") 65endfunction 66