[SciPy-User] frequency components of a signal buried in a noisy time domain signal

[Ivo Maljevic]

Nils,
I will not answer your questions to Anne then (I could, but you asked her :)
), but I'll give you another advice, given that you are new to signal
processing.
Whatever you do, please do not use the scipy's snr function call to
calculate the SNR in signal processing.

I'll illustrate the problem with the 'statistical SNR' in signal processing
with 3 closely related and simple examples.

1. Infinitely long sinusoidal signal with additive white noise

Let x(t)=s(t)+n(t)=A sin(2\pi f t) + n(t)

Power of the signal component can be found as P_s=1/(2T) \int_T^T s^2(t) dt
= A^2/2, and let the power of noise be P_n=var(n) = 1.

If, we chose A=100, then intuitively, we know that the SNR is large, and
indeed: SNR = P_s/ P_N = A^2/2 = 5000

On the other hand, if you use statistical definition (used in optical
astronomy, and scipy and implemented in SCIPY), you would conclude quite the
opposite:

SNR = \mu_s / signa_n = mean(s) / signa_n= 0/1 = 0, and that would be the
case regardless of what the value for A is.

2. Discrete time signal with finite duration

Let's take this same signal and sample it in a such way that we have 100
cycles of a sinusoid:


>>> f=200
>>> t=np.linspace(0,0.5,2000)
>>> s=100*np.sin(2*np.pi*f*t)

Using sampled signal definitions for powers, you get:

P_s = np.mean(s**2) = 4997.5

and you can calculate the power of the noise either as:

>>> np.var(n) = 0.97726608241925972

or

the same way as the power of the signal:

>>> np.mean(n**2) = 0.977977713

The results will slightly differ, and that has to do with biased/unbiased
estimator.

The SNR=4997.5/0.9772660=5113.756

which is pretty close to analog domain case with infinitely long signal.

Using the scipy approach, you would however conclude the following:

SNR=mean(s)/sigma_n=-5.6495516860933771e-15/0.9885676923808808=
-5.7148860210947362e-15 (which comes close to zero).

3. Same example, but let's take fractional number of cycles.

>>> t=np.linspace(0,0.5+np.pi/800,2000+int(5*np.pi)) # this is 100 cycles +
1/4 of a cycle
>>> s=100*np.sin(2*np.pi*f*t)

>>> P_s = np.mean(s**2)
>>> SNR=P_s/np.var(n)
>>> SNR
5117.8701015651177  # very close to the two results above

But, using the other definition:

>>> SNR=np.mean(s)/np.std(n)
>>> SNR
0.098933308385744489

As you can see, your SNR would be a function of the frequency of the signal
(that is, it would matter how you truncate the cycles) and it wouldn't
matter much what is the amplitude of the signal.

Just my 2 cents,

Ivo



On 27 February 2010 04:41, Nils Wagner <nwagner at iam.uni-stuttgart.de> wrote:

> On Fri, 26 Feb 2010 16:32:01 -0500
>  Anne Archibald <peridot.faceted at gmail.com> wrote:
> > Hi,
> >
> > Looking at a periodic signal buried in noise is a
> >well-studied
> > problem, with many techniques for attacking it. You
> >really need to be
> > a little more specific about what you want to do. For
> >example, is your
> > input signal really a sinusoid, or does it have harmonic
> >content? Are
> > you trying to detect a weak periodic signal or are you
> >trying to
> > extract the features of a strong periodic signal? Is
> >your signal
> > exactly periodic, does it have some (deterministic or
> >random) wander,
> > or are you looking for the power spectrum of a broadband
> >signal?
> >
> > If your input data are non-uniformly sampled, everything
> >becomes more
> > difficult (and computationally expensive), but there are
> >solutions
> > (e.g. the Lomb-Scargle periodogram).
> >
> > Anne
> >
> Hi Anne,
>
> Thank you very much for your hints !
>
> BTW, a BSD licensed code for the Lomb-Scargle periodogram
> is available at
> http://www.mathworks.com/matlabcentral/fileexchange/993-lombscargle-m
>
> http://www.mathworks.com/matlabcentral/fileexchange/20004-lomb-lomb-scargle-periodogram
>
>
> I am newbie to signal processing.
> Is there a good introduction that you can recommend ?
> There are so many books on signal processing. It should
> cover engineering applications.
>
> What makes a signal weak/strong periodic ?
>
> The signals come from real-life application (pressure /
> acceleration data).
> Do I need a filter before I apply FFT ?
>
> What would you do if you know nothing about the origin of
> the signal ?
>
> How can I distinguish between deterministic and random
> wander ?
>
> Cheers,
>                        Nils
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100227/5b3bd3cf/attachment.html>