Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Different SNR estimation with FFT and Welch, which is the problem?

Asked by Gianfranco on 1 May 2013

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

0 Comments

Gianfranco

Tags

Products

No products are associated with this question.

2 Answers

Answer by 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)

0 Comments

Youssef  KHMOU
Answer by 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;

0 Comments

Youssef  KHMOU

Contact us