Python  - SciPy                                                         Home : www.sharetechnote.com

SciPy - Basics

Importing the package

you can import a specific subpackage you want to use as shown below

 import numpy as np import matplotlib.pyplot as plt from scipy.fftpack import fft   t = np.linspace(0,10*np.pi,100); sig = np.cos(3*t);   sig_freq = fft(sig);   plt.plot(t,sig_freq,'r-'); plt.show();

or you can import the full package as shown below.

 import numpy as np import matplotlib.pyplot as plt import scipy as sp   t = np.linspace(0,10*np.pi,100); sig = np.cos(3*t);   sig_freq = sp.fft(sig);   plt.plot(t,sig_freq,'r-'); plt.show();

Finding Peaks

< Find all the peaks >

 import numpy as np import matplotlib.pyplot as plt import scipy as sp from scipy.signal import find_peaks     t = np.linspace(0,10*np.pi,200); x = np.cos(t);   peaks,_ = find_peaks(x);   plt.plot(t, x,'b-'); plt.scatter(t[peaks],x[peaks],c='red'); plt.show();

 import numpy as np import matplotlib.pyplot as plt import scipy as sp from scipy.signal import find_peaks     t = np.linspace(0,10*np.pi,200); x = np.cos(t) + 0.3*np.sin(5*t);   peaks,_ = find_peaks(x);   plt.plot(t, x,'b-'); plt.scatter(t[peaks],x[peaks],c='red'); plt.show();

< Find the peaks above a threshold >

 import numpy as np import matplotlib.pyplot as plt import scipy as sp from scipy.signal import find_peaks     t = np.linspace(0,10*np.pi,200); x = np.cos(t) + 0.3*np.sin(5*t);   peaks,_ = find_peaks(x, height = 0.5);   plt.plot(t, x,'b-'); plt.scatter(t[peaks],x[peaks],c='red'); plt.show();

< Find all the peaks on ECG >

 import numpy as np import matplotlib.pyplot as plt import scipy as sp from scipy.misc import electrocardiogram from scipy.signal import find_peaks     fs = 360;   x = electrocardiogram()[0:500]; t = (np.linspace(0,x.size,x.size) / fs)[0:500];   peaks,_ = find_peaks(x);   plt.plot(t, x,'b-'); plt.scatter(t[peaks],x[peaks],c='red'); plt.show();

 import numpy as np import matplotlib.pyplot as plt import scipy as sp from scipy.misc import electrocardiogram from scipy.signal import find_peaks     fs = 360;   x = electrocardiogram()[0:500]; t = (np.linspace(0,x.size,x.size) / fs)[0:500];   peaks,_ = find_peaks(x);   plt.plot(t, x,'b-'); plt.scatter(t[peaks],x[peaks],c='red'); plt.show();

ECG (ElectroCardiogram)

 import numpy as np import matplotlib.pyplot as plt import scipy as sp import scipy.misc as spm   x = spm.electrocardiogram(); plt.plot(x,'r-'); plt.show();

 import numpy as np import matplotlib.pyplot as plt import scipy as sp import scipy.misc as spm   fs = 360;   x = spm.electrocardiogram(); t = np.linspace(0,x.size,x.size) / fs;   plt.plot(t, x,'r-'); plt.xlabel('time in s'); plt.ylabel('ECG in mV'); plt.show();

 import numpy as np import matplotlib.pyplot as plt import scipy as sp import scipy.misc as spm   fs = 360;   x = spm.electrocardiogram(); t = np.linspace(0,x.size,x.size) / fs;   plt.plot(t, x,'r-'); plt.xlim(1.0,5.0); plt.ylim(-2,2); plt.xlabel('time in s'); plt.ylabel('ECG in mV'); plt.show();

Low Pass Filter

 import numpy as np import matplotlib.pyplot as plt import scipy as sp from scipy import signal   n = 100; b = signal.firwin(n,cutoff = 0.2, window = "hamming"); plt.stem(b); plt.show();

 import numpy as np import matplotlib.pyplot as plt import scipy as sp from scipy import signal   n = 100; a = 1; b = signal.firwin(n,cutoff = 0.2, window = "hamming"); w,h = signal.freqz(b,a); h_dB = 20 * np.log10(abs(h)) plt.plot(w/np.pi,h_dB); plt.show();

High Pass Filter

 import numpy as np import matplotlib.pyplot as plt import scipy as sp from scipy import signal   n = 101; b = signal.firwin(n,cutoff = 0.8, window = "hamming",pass_zero = False); plt.stem(b,use_line_collection = True); plt.show();

 import numpy as np import matplotlib.pyplot as plt import scipy as sp from scipy import signal   n = 101; a = 1; b = signal.firwin(n,cutoff = 0.8, window = "hamming",pass_zero = False); w,h = signal.freqz(b,a); h_dB = 20 * np.log10(abs(h)) plt.plot(w/np.pi,h_dB); plt.show();

Band Pass Filter

 import numpy as np import matplotlib.pyplot as plt import scipy as sp from scipy import signal   n = 101; b = signal.firwin(n,cutoff = [0.4,0.6], window = "hamming",pass_zero = False);   plt.stem(b,use_line_collection = True); plt.show();

 import numpy as np import matplotlib.pyplot as plt import scipy as sp from scipy import signal   n = 101; a = 1; b = signal.firwin(n,cutoff = [0.4,0.6], window = "hamming",pass_zero = False);   w,h = signal.freqz(b,a); h_dB = 20 * np.log10(abs(h)) plt.plot(w/np.pi,h_dB); plt.show();

 import numpy as np import matplotlib.pyplot as plt import scipy as sp from scipy import signal   n = 101; b = signal.firwin(n,cutoff = [0.4,0.6], window = "hamming");   plt.stem(b,use_line_collection = True); plt.show();

 import numpy as np import matplotlib.pyplot as plt import scipy as sp from scipy import signal   n = 101; a = 1; b = signal.firwin(n,cutoff = [0.4,0.6], window = "hamming");   w,h = signal.freqz(b,a); h_dB = 20 * np.log10(abs(h)) plt.plot(w/np.pi,h_dB); plt.show();

Applying Filter to Signal

 import numpy as np import matplotlib.pyplot as plt import scipy as sp from scipy import signal   t = np.linspace(0,10*np.pi,200); sig = np.sin(t) + 0.2 * np.sin(10*t);   n = 10; filter_taps = signal.firwin(n,cutoff = 0.2, window = "hamming");   sig_filtered = signal.lfilter(filter_taps, 1.0, sig)   plt.subplot(3,1,1); plt.plot(sig);   plt.subplot(3,1,2); plt.stem(filter_taps,use_line_collection = True);   plt.subplot(3,1,3); plt.plot(sig_filtered);   plt.show();

Detrending

< Detrending with detrend() >

 import numpy as np import matplotlib.pyplot as plt import scipy as sp from scipy import signal     t = np.linspace(0,10*np.pi,200); x = 0.2*t + np.random.normal(size=len(t));   x_detrended = signal.detrend(x);   plt.plot(t, x,'b-'); plt.plot(t, x_detrended,'r-'); plt.show();

 import numpy as np import matplotlib.pyplot as plt import scipy as sp from scipy import signal     t = np.linspace(0,10*np.pi,200); x = -0.5*pow(t + 2,2) + np.random.normal(size=len(t));   x_detrended = signal.detrend(x);   plt.plot(t, x,'b-'); plt.plot(t, x_detrended,'r-'); plt.show();

< Detrending with convolution() >

 import numpy as np import matplotlib.pyplot as plt import scipy as sp from scipy import signal     x = -0.5*pow(t + 2,2) + np.random.normal(size=len(t));   w = [1,-1]; x_detrended = np.convolve(x,w,mode='same');   plt.plot(t, x,'b-'); plt.plot(t, x_detrended,'r-'); plt.show();

 import numpy as np import matplotlib.pyplot as plt import scipy as sp from scipy import signal     x =  np.cos(t) + 0.3*np.sin(5*t);   w = [1,-1]; x_detrended = np.convolve(x,w,mode='same');   plt.plot(t, x,'b-'); plt.plot(t, x_detrended,'r-'); plt.show();

 import numpy as np import matplotlib.pyplot as plt import scipy as sp from scipy import signal     x =  np.cos(t) + 0.3*np.sin(5*t); x += 0.05*np.random.normal(size=len(t));   w = [1,-1]; x_detrended = np.convolve(x,w,mode='same');   plt.plot(t, x,'b-'); plt.plot(t, x_detrended,'r-'); plt.show();

Welch

 import numpy as np import matplotlib.pyplot as plt from scipy import signal, fft   N = 1000000; fs = 100000; dt = 1.0/fs; fc = 1234.0; amp = 2*np.sqrt(2); t = dt * np.arange(N); sig = np.cos(2*np.pi*fc*t) + 0.5*np.cos(3*2*np.pi*fc*t) + 0.2*np.cos(5*2*np.pi*fc*t); sig = amp * sig;   freq,Pd = signal.welch(sig,fs);   plt.subplot(2,1,1); plt.plot(freq,Pd);   plt.subplot(2,1,2); plt.plot(freq,20*np.log10(Pd));   plt.show();

 import numpy as np import matplotlib.pyplot as plt from scipy import signal, fft   N = 1000000; fs = 100000; dt = 1.0/fs; fc = 1234.0; amp = 2*np.sqrt(2); t = dt * np.arange(N); sig = np.cos(2*np.pi*fc*t) + 0.5*np.cos(3*2*np.pi*fc*t) + 0.2*np.cos(5*2*np.pi*fc*t); sig = amp * sig;   freq,Pd = signal.welch(sig,fs,nperseg=1024);  # change resolution with nperseg option   plt.subplot(2,1,1); plt.plot(freq,Pd);   plt.subplot(2,1,2); plt.plot(freq,20*np.log10(Pd));   plt.show();

 import numpy as np import matplotlib.pyplot as plt from scipy import signal, fft   N = 1000000; fs = 100000; dt = 1.0/fs; fc = 1234.0; amp = 2*np.sqrt(2); n_amp = 0.0001 * fs/2; t = dt * np.arange(N); sig = np.cos(2*np.pi*fc*t) + 0.5*np.cos(3*2*np.pi*fc*t) + 0.2*np.cos(5*2*np.pi*fc*t); sig = amp * sig; sig = sig + np.random.normal(scale=n_amp, size=t.shape); # add noise for noise floor   freq,Pd = signal.welch(sig,fs,nperseg=1024);  # change resolution with nperseg option   plt.subplot(2,1,1); plt.plot(freq,Pd);   plt.subplot(2,1,2); plt.plot(freq,20*np.log10(Pd));   plt.show();

Reference :

[1]