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

New to MATLAB?

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

Asked by Gianfranco

Gianfranco (view profile)

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

Gianfranco (view profile)

Tags

Products

No products are associated with this question.

2 Answers

Answer by Youssef Khmou

Youssef Khmou (view profile)

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

Youssef Khmou (view profile)

Answer by Youssef Khmou

Youssef Khmou (view profile)

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

Youssef Khmou (view profile)

Contact us