Sampling, clipping and aliasing#

import io
from nbformat import read

def execute_notebook(nbfile):
    
    with io.open(nbfile) as f:
        nb = read(f,3)
    
    ip = get_ipython()
    
    for cell in nb.worksheets[0].cells:
        if cell.cell_type != 'code':
            continue
        ip.run_cell(cell.input)
        
execute_notebook("create_plot_signal.ipynb") # also imports numpy as np and matplotlib.pyplot as plt
[<matplotlib.lines.Line2D at 0x7f6e0c65a830>]
../_images/fd52d5c817af87716c923293954c7abf88c18fe2646b733dba8226149ddc6923.png
# example
t = np.linspace(0,10,10000) # almost continuous
y = 9+np.sin(2*np.pi*0.1*t)

ts,yq,tr,yr = adc(t,y,fs=1,N=4,miny=0,maxy=10,method='soh') # monopolar
plt.figure()
plt.plot(t,y,'k--',lw=0.1)
plt.plot(ts,yq,'ro')
plt.plot(tr, yr,'b-')
[<matplotlib.lines.Line2D at 0x7f6e0a3d0bb0>]
../_images/12a63fc37087f436986c837a221c131ec7eac1cfc4e94338a8143a2e9669bdb9.png
# an example from the A/D lecture
t = np.linspace(0,1,1000) # almost continuous
y = 3 + 3*np.sin(2*np.pi*10*t)# - np.pi/2.)

sample at 15 Hz#

ts,yq,tr,yr = adc(t,y,fs=15.,N=24,miny=0,maxy=10,method=None) # monopolar
plt.figure(figsize=(10,8))
plt.plot(t,y,'m--',lw=0.1)
plt.plot(ts,yq,'ro')
plt.plot(tr, yr,'b-')
plt.xlabel('$t$ [s]',fontsize=18)
plt.ylabel('$y$ [V]',fontsize=18)
plt.xlim([0,1.0])
plt.title('$f_s = 15 $ Hz ',fontsize=22)
Text(0.5, 1.0, '$f_s = 15 $ Hz ')
../_images/664543e4c634aff97ffa427942d610a6464bde9d96cf8e7bf297dc1f8e746f77.png

sample at 11 Hz#

ts,yq,tr,yr = adc(t,y,fs=11.,N=24,miny=0,maxy=10) # monopolar
plt.figure(figsize=(10,8))
plt.plot(t,y,'k--',lw=0.1)
plt.plot(ts,yq,'ro')
plt.plot(tr, yr,'b-',lw=.5)
plt.xlabel('$t$ [s]',fontsize=18)
plt.ylabel('$y$ [V]',fontsize=18)
plt.xlim([0,1.0])
plt.title('$f_s = 11$ Hz',fontsize=22)
Text(0.5, 1.0, '$f_s = 11$ Hz')
../_images/ae7e40b478dcca1d95d5d01af8a89074112e81733c0d6447a888979353b1b363.png
ts,yq,tr,yr = adc(t,y,fs=9.,N=24,miny=0,maxy=10) # monopolar
plt.figure(figsize=(10,8))
plt.plot(t,y,'k--',lw=0.1)
plt.plot(ts,yq,'ro')
plt.plot(tr, yr,'b-',lw=.5)
plt.xlabel('$t$ [s]',fontsize=18)
plt.ylabel('$y$ [V]',fontsize=18)
plt.xlim([0,1.0])
plt.title('$f_s = 9$ Hz',fontsize=22)
Text(0.5, 1.0, '$f_s = 9$ Hz')
../_images/95093c96451a7b51daadaa0e7480664e71322ccf611cecc8d7fb5ee7c0434b10.png
ts,yq,tr,yr = adc(t,y,fs=6.,N=24,miny=0,maxy=10) # monopolar
plt.figure(figsize=(10,8))
plt.plot(t,y,'k--',lw=0.1)
plt.plot(ts,yq,'ro')
plt.plot(tr, yr,'b-',lw=.5)
plt.xlabel('$t$ [s]',fontsize=18)
plt.ylabel('$y$ [V]',fontsize=18)
plt.xlim([0,1.0])
plt.title('$f_s = 6$ Hz',fontsize=22)
Text(0.5, 1.0, '$f_s = 6$ Hz')
../_images/616695adf02b2d1e8d254225186b61845912293f6defa6f6f332d4fc9bd03fb1.png

Estimate the aliasing using the folding diagram or the formula:#

if \(f_s > 2 f\) , no aliasing

if \(2/3 f < f_s < 2 f\), \(f_a = |f_s - f|\)

if \(f_s < 2/3 f\), \(f_a = (f/f_f)f_f\), where \(f_f = f_s/2\)

in most cases: $\(f_a = \left|f- f_s \cdot NINT (f/f_s) \right|\)$

How do we use it? we use trial/error and smart guesses:#

let’s check different options:

if it’s aliasing we could be in one of the two regions, either 6 Hz < 2/3 f, i.e. f > 1.5*6, f > 9 or

it can be that 6Hz is above 2/3f and below 2f, i.e. f is between 3 and 9 Hz

theoretically speaking if f is below 3Hz, then we’re above the Nyquist frequency.

So, we see 2 Hz (2 peaks in 1 second) so it can be that we sampled correctly using 6 Hz or incorrectly and got aliasing.

Therefore we could test few things, but first let’s prepare some options:

  1. if we are fine and 2Hz is there, so we just see not a nice sine signal and we can increase the sampling frequency and get it nice

  2. if we are aliasing, then it can be what? 2Hz = |f - 6Hz * NINT(f/6Hz)|, so f = 6Hz * NINT(f/6Hz) + 2Hz. Let’s see if 10 Hz is reasonable: 10/6 = 1.6667 and the nearest integer is 2 and it means that 10 - 6*2 = 2 Hz.

  3. it can be also lower than 6Hz and we got 2Hz for instance if we sampled f = 4Hz

np.abs(10 - 6 * np.round(10./6.))
2.0
ts,yq,tr,yr = adc(t,y,fs=4.,N=24,miny=0,maxy=10) # monopolar
plt.figure(figsize=(10,8))
plt.plot(t,y,'k--',lw=0.1)
plt.plot(ts,yq,'ro')
plt.plot(tr, yr,'b-',lw=.5)
plt.xlabel('$t$ [s]',fontsize=18)
plt.ylabel('$y$ [V]',fontsize=18)
plt.xlim([0,1.0])
plt.title('$f_s = 4$ Hz',fontsize=22)
Text(0.5, 1.0, '$f_s = 4$ Hz')
../_images/82cd4056285a14f3ff9cb75d968c85502a8a0033d9c4f417156416a97c519940.png