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();

 

 

< Find the peaks above a threshold 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();

 

 

 

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();

 

 

 

Band Stop 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");

 

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]