The transmitting signal is $$ s(t) = \sin(20\pi t)+\sin(40\pi t)+\sin(60\pi t) $$

The received signal after AWGN channel is $$ r(t) = s(t)+n(t) $$ where the variance of noise $n(t)$ is $\sigma^2$.

The theoretical single spectrum of the $r(t)$ should be

$$ R(f) = \pi\delta(f-10)+ \pi\delta(f-20)+ \pi\delta(f-30)+\sigma^2 $$

where $\delta(.)$ is Delta function. However, the spectrum obtained through MATLAB seems not like the theroetical one.

enter image description here

When the SNR is setted as lower than -25dB, the spikes will be submerged by noise. And how to calculate the threshold SNR where sinusoid part is totally covered by noise?

enter image description here

The relative code is here

fs_init = 10; % basis frequency
T = 2*pi/om; 
tall = 10;
fs_init= 10;
ut = zeros(1,length(t));

%%% Added constants
NSyboml = 3;

for i = 1:NSyboml
    ut = ut+a*sin(2*pi*i*fs_init*t);

SNR = 0;
rt = awgn(ut,SNR,'measured');   % add noise

% FFT 
y = rt;
L = length(ut);    % Length of signal
Ts_v = t./(0:L-1); % sample time
Fs = ceil(1/Ts_v(2));    % sample frequency
deltaf = 0.1; % resolution
NFFT = floor(Fs/deltaf);
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,floor(NFFT/2+1));

% Plot single-sided amplitude spectrum.
 figure;set(gcf,'Color','w','Position',[0 0 500 350]);
xlabel('Frequency (Hz)');ylabel('|R(f)|');
grid on;
title(['SNR = ',num2str(SNR),'dB'])
%%%enter preformatted text here
  • 1
    $\begingroup$ I slightly changed your code to make it work. I hope it is fine $\endgroup$ – Laurent Duval Nov 5 '18 at 12:56
  • 2
    $\begingroup$ your theoretical formula is wrong. The impulses are expected at $f = 5, 10, 15$ Hz locations... $\endgroup$ – Fat32 Nov 5 '18 at 14:18
  • 1
    $\begingroup$ I would have slept more this week-end. Busy at trying the code, I did not notice that... $\endgroup$ – Laurent Duval Nov 5 '18 at 14:24

Notice that in your first image, there are peaks where they would be expected to be, so that's fine. I guess you think that "the spectrum obtained through MATLAB seems not like the theroetical one" because of the noisy spectrum you get, instead of a straight line at $\sigma^2$ on the vertical axis.

The random signal added to $s(t)$ is not random actually, it's just a realization you get when calling awgn(). So the FFT will return just a noisy spectrum for that signal. If you generate lots of realizations of $r(t)$ and plot the average of all their FFTs, you'll obtain what you were expecting for.

For your second question on how to calculate the threshold SNR, some criteria is needed. If you average lots of FFTs the variance of the noisy part of the spectrum will drastically decrease, so the SNR needed to see the peaks will be much smaller. If you use just one realization, then the SNR will be higher. There is not a fixed number for the needed SNR, it depends on what you refer to "sinusoid covered by noise". If the spectrum observed was just like the theoretical one (i.e. you averaged a lot of FFTs), a sinusoid would be "completely covered by noise" when its amplitude is such that the corresponding peak in the frequency domain is below $\sigma^2$.


Your question is related to a central topic: given the theoretical power spectral density (DSP) of a class of signals (your formula), how well can it be obtained in practice from a given practical realization?

This can be tricky, as one the one hand you have the DSP of a continuous infinite stochastic process, and one the other hand, a finite number of samples in a limited window frame. Likewise, the Fourier transform of a random noise does not make sense in general, so a mere FFT may not yield what you expected. For instance, the FFT of a vector of generated white Gaussian is almost never flat.

So, the two main options to improve the estimation is to:

  • use models, and perform parametric spectral analysis
  • use some statistics, and study the appropriateness of non-parametric spectral analysis, and how you can improve the results with windowing, averaging, etc. for instance with periodograms.

I would suggest you to follow the slides "Spectrum estimation using Periodogram, Bartlett and Welch" by Guido Schuster. They "follow closely chapter 8 in the book "Statistical Digital Signal Processing and Modeling" by Monson H. Hayes", and I use them directly in a lecture to give a quick overview of the limits of "just plot the FFT of my signal".

Another support is: Lecture 3 - Spectrum Estimation, Danilo Mandic, Adaptive Signal Processing & Machine Intelligence.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.