Different SNR estimation with FFT and Welch, which is the problem?
Show older comments
I wrote this code in which I calculate a noisy signal with a specified SNR and successively I try to estimate the PSD of the signal and of the noise in order to verify the SNR.
Fs=96e6; fc=24e6; B=7.92e6; N=2401; span=48e6; snr=5; deltaf=span/N; Ts=1/Fs; Toss=1/deltaf; t = -Toss/2:Ts:Toss/2; %x=sinc(t*B).*cos(2*pi*fc*t); x=2*cos(2*pi*fc*t); for i=1:100, y=awgn(x,snr,'measured'); n=y-x; nfft = length(x); hp = spectrum.welch('hann',nfft,10); hpopts=psdopts(hp,x); Hpsd = psd(hp,x,hpopts,'NFFT',nfft,'Fs',Fs); hpopts=psdopts(hp,n); Hpsd_noise = psd(hp,n,hpopts,'NFFT',nfft,'Fs',Fs);
signal_power=Hpsd.avgpower; noise_power=Hpsd_noise.avgpower; SNR(i)=10*log10(signal_power/noise_power); Pxx1 = abs(fft(x,nfft)).^2/length(x)/(Fs); Pnoise1 = abs(fft(n,nfft)).^2/length(n)/(Fs); Hpsd1 = dspdata.psd(Pxx1(1:length(Pxx1)/2)*2,'Fs',Fs); Hpsd1_noise = dspdata.psd(Pnoise1(1:length(Pnoise1)/2)*2,'Fs',Fs); signal_power1=Hpsd1.avgpower; noise_power1=Hpsd1_noise.avgpower; SNR1(i)=10*log10(signal_power1/noise_power1); end mean(SNR) mean(SNR1)
When the input signal is a single tone, I have results using the two PSD estimators (FFT and Welch). Unfortunatly it is not the same when I consider a different signal. What I'm not considering or what is wrong in my code. Thank You
Answers (2)
Youssef Khmou
on 1 May 2013
hi, i tried but still can not find the reason of the doubled SNR, i format the code for next time, however the signals below give the correct SNR :
%x=cos(2*pi*fc*t); % the signal in vlots
%x=sin(t*B);
%x=square(2*pi*fc*t);
%x=sawtooth(fc*t,.2);
While these signals give #SNRs :
%x=chirp(t,0,fc,t(end));
%x=sinc(fc*t);
Reformatted code :
clear;
Fs=96e6; % samping frequency in Mhz
fc=24e6; % carrier frequency
N=2401; % Number os samples
span=48e6;
B=7.92e6;
snr=5; % signal to noise ratio in dB
deltaf=span/N;
Ts=1/Fs;
Toss=1/deltaf;
t=-Toss/2:Ts:Toss/2; % time simulation
x=cos(2*pi*fc*t); % the signal in vlots
for i=1:100
y=awgn(x,snr,'measured');
%y=x+sigma*randn(size(x));
n=y-x;
nfft = length(x);
hp = spectrum.welch('hann',nfft,10);
hpopts=psdopts(hp,x);
%set(hpopts,'SpectrumType','TwoSided');
Hpsd = psd(hp,x,hpopts,'NFFT',nfft,'Fs',Fs);
hpopts=psdopts(hp,n);
Hpsd_noise = psd(hp,n,hpopts,'NFFT',nfft,'Fs',Fs);
signal_power=Hpsd.avgpower;
noise_power=Hpsd_noise.avgpower;
SNR(i)=10*log10(signal_power/noise_power);
Pxx1 = abs(fft(x,nfft)).^2/length(x)/(Fs);
Pnoise1 = abs(fft(n,nfft)).^2/length(n)/(Fs);
Hpsd1 = dspdata.psd(Pxx1(1:end/2)*2,'Fs',Fs);
Hpsd1_noise = dspdata.psd(Pnoise1(1:end/2)*2,'Fs',Fs);
signal_power1=Hpsd1.avgpower;
noise_power1=Hpsd1_noise.avgpower;
SNR1(i)=10*log10(signal_power1/noise_power1);
end
mean(SNR)
mean(SNR1)
Youssef Khmou
on 1 May 2013
hi,
So far to have the correct SNR for the sinc(B*t) signal change spectrum method to periodogram :
%hp = spectrum.welch('hann',nfft,10);
hp = spectrum.periodogram;
Categories
Find more on Parametric Spectral Estimation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!