Scipy code examples and exercises

In [2]:
import numpy as np
import scipy.signal as signal
In [114]:
%matplotlib inline
import matplotlib.pyplot as plt

nyq_freq = 1
cutoff1 = 0.5
cutoff2 = np.array([0.3, 0.8])
cutoff_normalized1 = cutoff1*2*np.pi/(2*nyq_freq)
cutoff_normalized2 = cutoff2*2*np.pi/(2*nyq_freq)

b1 = signal.firwin(40, 0.5, nyq = nyq_freq)
b2 = signal.firwin(41, [0.3, 0.8], nyq = nyq_freq)
w1, h1 = signal.freqz(b1)
w2, h2 = signal.freqz(b2)

line_x1 = [cutoff_normalized1, cutoff_normalized1]
line_y = [20, -140]
line_x2 = [cutoff_normalized2[0], cutoff_normalized2[0]]
line_x3 = [cutoff_normalized2[1], cutoff_normalized2[1]]
plt.title('Digital filter frequency response')
plt.plot(w1, 20*np.log10(np.abs(h1)), 'b')
plt.plot(w2, 20*np.log10(np.abs(h2)), 'r')
plt.plot(line_x1, line_y, 'b--')
plt.plot(line_x2, line_y, 'r--')
plt.plot(line_x3, line_y, 'r--')
plt.ylabel('Amplitude Response (dB)')
plt.xlabel('Frequency (rad/sample)')
plt.grid()
plt.show()
In [39]:
sample_rate = 100
nsamples = 400

# signal generation
t = np.arange(nsamples) / sample_rate
x = np.cos(2*np.pi*0.5*t) + 0.2*np.sin(2*np.pi*2.5*t+0.1) + \
    0.2*np.sin(2*np.pi*15.3*t) + 0.1*np.sin(2*np.pi*16.7*t + 0.1) + \
    0.1*np.sin(2*np.pi*23.45*t+.8)

# FIR filter creation
nyq_rate = sample_rate / 2.0
cutoff_hz = 10.0
N = 70
taps = signal.firwin(N, cutoff_hz/nyq_rate)

# filter x with FIR filter
filtered_x = signal.lfilter(taps, 1.0, x)

# FIR filter coefficients plot
plt.figure(1)
plt.plot(taps, 'bo-', linewidth=2)
plt.title('Filter Coefficients (%d taps)' % N)
plt.grid(True)
plt.show()

# plot of frequency response of FIR filter
plt.figure(2)
w, h = signal.freqz(taps, worN=8000)
plt.plot((w/np.pi)*nyq_rate, np.absolute(h), linewidth=2)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Gain')
plt.title('Frequency Response')
plt.ylim(-0.05, 1.05)
plt.grid(True)

# The phase delay of the filtered signal
delay = 0.5 * (N-1) / sample_rate
plt.figure(3)
# Plot the original signal
plt.plot(t, x)
# Plot the filtered signal, shifted to compensate for the phase delay.
plt.plot(t-delay, filtered_x, 'r-')
# Plot just the "good" part of the filtered signal  The first N-1
# samples are "corrupted" by the initial conditions
plt.plot(t[N-1:]-delay, filtered_x[N-1:], 'g', linewidth=4)
plt.xlabel('t')
plt.grid(True)
plt.show()
In [132]:
fs = 10e3
N = 1e5
amp = 2*np.sqrt(2)
freq = 1270.0
noise_power = 0.001 * fs / 2
time = np.arange(N) / fs
x = amp*np.sin(2*np.pi*freq*time)
x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)

f, Pper_spec = signal.periodogram(x, fs, 'flattop', scaling='spectrum')

plt.semilogy(f, Pper_spec)
plt.xlabel('frequency [Hz]')
plt.ylabel('PSD')
plt.grid()
plt.show()
In [133]:
fs = 10e3
N = 1e5
amp = 2*np.sqrt(2)
freq = 1270.0
noise_power = 0.001 * fs / 2
time = np.arange(N) / fs
x = amp*np.sin(2*np.pi*freq*time)
x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)

f, Pwelch_spec = signal.welch(x, fs, scaling='spectrum')

plt.semilogy(f, Pwelch_spec)
plt.xlabel('frequency [Hz]')
plt.ylabel('PSD')
plt.grid()
plt.show()
In [165]:
import scipy.misc as misc
image = misc.face(gray=True)
w = np.zeros((50, 50))
w[0][0] = 1.0
w[49][25] = 1.0
image_new = signal.fftconvolve(image, w)

plt.figure()
plt.imshow(image)
plt.gray()
plt.title('Original image')
plt.show()
In [162]:
plt.figure()
plt.imshow(image_new)
plt.gray()
plt.title('Filtered image')
plt.show()
In [49]:
a = np.array([1, 2, 3, 4])
np.save("array_file.npy", a)
b = np.load("array_file.npy")
print(b)
[1 2 3 4]
In [50]:
b *= 6
np.savez("dict_array.npz", ar1=a, ar2=b, ar3=a+b)
data = np.load("dict_array.npz")
print(data['ar1'])
print(data['ar2'])
print(data['ar3'])
[1 2 3 4]
[ 6 12 18 24]
[ 7 14 21 28]
In [52]:
np.savetxt("text_array.txt", data['ar3'])
c = np.loadtxt("text_array.txt")
print(c)
[  7.  14.  21.  28.]
In [24]:
sig = np.loadtxt("ecg.txt")

sample_rate = 1000.0
nyq_rate = sample_rate / 2.0
cutoff_hz = 50.0
nsamples = len(sig)
t = np.arange(nsamples) / sample_rate
N = 101

taps = signal.firwin(N, [cutoff_hz/nyq_rate*0.7, cutoff_hz/nyq_rate*1.3])
filtered_x = signal.lfilter(taps, 1.0, sig)
delay = 0.5 * (N-1) / sample_rate

plt.figure(1)
w, h = signal.freqz(taps, worN=8000)
plt.plot((w/np.pi)*nyq_rate, np.absolute(h), linewidth=2)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Gain')
plt.title('Frequency Response')
plt.ylim(-0.05, 1.05)
plt.grid(True)

plt.figure(2)
plt.plot(t, sig)
plt.plot(t[N-1:]-delay, filtered_x[N-1:], 'r')
plt.axis([0, 1, min(sig), max(sig)])
plt.grid()
plt.show()
In [25]:
f = open("example", 'w+')
f.write("File content")
f.seek(0, 0)
data1 = f.read(4)
f.seek(9, 0)
data2 = f.read(3)
f.close()
print(data1)
print(data2)
File
ent
In [115]:
nyq_freq = 1
freqs = [0.0, 0.3, 0.6, 1.0]
gains = [1.0, 2.0, 0.5, 0.0]
freqs_normalized = [i*2*np.pi/(2*nyq_freq) for i in freqs]

b = signal.firwin2(150, freqs, gains, nyq = nyq_freq)
w, h = signal.freqz(b)

plt.title('Digital filter frequency response')
plt.plot(w, np.abs(h))
line_y = [min(np.abs(h)), max(np.abs(h))]
for i in freqs_normalized:
    line_x = [i, i]
    plt.plot(line_x, line_y, 'b--')
    
plt.title('Digital filter frequency response')
plt.ylabel('Amplitude Response')
plt.xlabel('Frequency (rad/sample)')
plt.grid()
plt.show()
In [118]:
b, a = signal.iirfilter(4, Wn=0.2, rp=5, rs=60, btype='lowpass', ftype='ellip')
w, h = signal.freqz(b, a)

plt.title('Digital filter frequency response')
plt.plot(w, 20*np.log10(np.abs(h)))
plt.title('Digital filter frequency response')
plt.ylabel('Amplitude Response [dB]')
plt.xlabel('Frequency (rad/sample)')
plt.grid()
plt.show()
In [125]:
b, a = signal.iirdesign(wp=100, ws=200, gpass=2.0, gstop=40., analog=True)
w, h = signal.freqs(b, a)

plt.title('Analog filter frequency response')
plt.plot(w, 20*np.log10(np.abs(h)))
plt.ylabel('Amplitude Response [dB]')
plt.xlabel('Frequency')
plt.grid()
plt.show()