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 0x7ff7edfbc8f0>]
../_images/964ca843707dd84aace96969b81791e8e68c7457e19ed1378d9cd2b8ce0048bd.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 0x7ff7eddc5190>]
../_images/b30d7b1d65f07b65cc54f223b8738df8cd4a3ff5b8762598823b676a2988feb4.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/a24fab013ba27e77d8708a28923ba61a6742011aa9eb0e804015bbd0c686ec30.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/6da9275b39ca34daffd8f84630f393f65e7663fa8601541701646695fd702ea0.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/29a8441e41e969b8044442ee49635971311dd353eae8abfc9e09594c432ebd88.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/668cdb3fd3179e9524c3d979ea979654f2d5a8b6bb9e891a0ef0a6260caf98df.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.))
np.float64(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/1276cab135bf372a61fe45b2056ec80e2c2d445178f4f78fbb9c5ea8b4063e39.png