Thursday, September 22, 2011

The sampling theorem explained with numpy

The sampling theorem states that a continuous signal x(t) bandlimited to B Hz can be recovered from its samples x[n] = x(n*T), where n is an integer, if T is greater than or equal to 1/(2B) without loss of any information. And we call 2B the Nyquist rate.
Sampling at a rate below the Nyquist rate is called undersampling, it leads to the aliasing effect. Let's observe the aliasing effect with the following script:
from numpy import linspace,cos,pi,ceil,floor,arange
from pylab import plot,show,axis

# sampling a signal badlimited to 40 Hz 
# with a sampling rate of 800 Hz
f = 40;  # Hz
tmin = -0.3;
tmax = 0.3;
t = linspace(tmin, tmax, 400);
x = cos(2*pi*t) + cos(2*pi*f*t); # signal sampling
plot(t, x)

# sampling the signal with a sampling rate of 80 Hz
# in this case, we are using the Nyquist rate.
T = 1/80.0;
nmin = ceil(tmin / T);
nmax = floor(tmax / T);
n = arange(nmin,nmax);
x1 = cos(2*pi*n*T) + cos(2*pi*f*n*T);
plot(n*T, x1, 'bo')

# sampling the signal with a sampling rate of 35 Hz
# note that 35 Hz is under the Nyquist rate.
T = 1/35.0;
nmin = ceil(tmin / T);
nmax = floor(tmax / T);
n = arange(nmin,nmax);
x2 = cos(2*pi*n*T) + cos(2*pi*f*n*T);
plot(n*T, x2, '-r.',markersize=8)

axis([-0.3, 0.3, -1.5, 2.3])
show()
The following figure is the result:
The blue curve is the original signal, the blue dots are the samples obtained with the Nyquist rate and the red dots are the samples obtainde with 35 Hz. It's easy to see that the blue samples are enough to recover the blue curve, while the red ones are not enough to capture the oscillations of the signal.

5 comments:

  1. Another wonderful post, full of new stuff to learn.

    Thanks!

    ReplyDelete
  2. Hello, first of all, thank you for your nice posts, they are very informative.

    On this one, I think it would be better to show the sampling effect on the frequency domain. You word it correctly, but it is easy to get the wrong message by that picture. For instance, try sampling with 85Hz... in the time domain you will not get a perfect fit, but in the frequency domain you will get a good approximation.

    Btw, 35Hz is even bellow the resolution of the signal itself. 60Hz would be a more appropriate example.

    Anyway, keep up the good work! :)

    ReplyDelete
  3. I like very much your work. Thanks for sharing!

    ReplyDelete
  4. I demonstrated the same effect in an interactive manner in matplotlib https://github.com/jlugao/aliasing_demo
    tell me what you think of it. It is still on beta testing as I still want to move it to a QT interface.

    ReplyDelete
    Replies
    1. Hey Joao, that great. I'll post it on twitter.

      Delete