Intelligent background subtraction algorithm using wavelets.
Also contains implementations of other wavelets and peak-finding code
Runs 2017/1/9
In [1]:
%matplotlib inline
from __future__ import division
import matplotlib.pyplot as plt
import numpy as np
import os
import sys
from scipy import signal
In [2]:
data1 = np.genfromtxt(os.path.join('..', 'tests', 'data', 'raman-785nm.txt'))
x = data1[:, 0]
y = data1[:, 1]
plt.plot(x, y)
Out[2]:
In [3]:
widths = np.arange(1,71)
cwtmat = signal.cwt(y, signal.ricker, widths)
plt.imshow(cwtmat, aspect='auto', cmap='PRGn')
Out[3]:
In [4]:
# Find local maxima
# make a binary array containing local maximum of transform, with same shape
lmax = np.zeros(cwtmat.shape)
for i in range(cwtmat.shape[0]):
lmax[i, signal.argrelextrema(cwtmat[i, :], np.greater)] = 1
fig, ax = plt.subplots(figsize=(15, 4))
ax.imshow(lmax, aspect='auto', cmap='gray_r')
Out[4]:
In [5]:
# allocate memory
# intial location assigned to peak from the first row
peak_pos_start = np.where(lmax[0,:]==1)[0]
# current position of the ridge
peak_ridge = np.copy(peak_pos_start) # full copy
n_peaks = peak_pos_start.size
# length of the ridge
peak_len = np.ones(n_peaks)
# use the max of the ridge line to find the width of the peaks
peak_pos = np.zeros(n_peaks, dtype='int')
peak_width = np.ones(n_peaks)
peak_width_max = np.zeros(n_peaks)
In [6]:
# Link local maxima (find ridges)
w = 3
# for each row starting at the second
for i in range(1, lmax.shape[0]):
# for each peak
for j in range(n_peaks):
# assume it doesn't extend, and then check
extends = False
p = peak_ridge[j]
if lmax[i, p] == 1:
# if there is one below, it is part of the same ridge
extends = True
else:
# if not search around peak
for k in range(1, w):
if lmax[i, p-k] == 1:
extends = True
peak_ridge[j] -= k
break
elif lmax[i, p+k] == 1:
extends = True
peak_ridge[j] += k
break
# if it extends
if extends:
# it it longer
peak_len[j] += 1
# find width by comparing max vs. previous
if cwtmat[i, p] > peak_width_max[j]:
peak_width_max[j] = cwtmat[i, p]
peak_width[j] = i
peak_pos[j] = p
In [7]:
print peak_pos[:20]
print peak_width[:20]
In [8]:
# generate a simulated spectrum of sorts, with peak positions and the length of the ridge lines
ypeaks = np.zeros(y.shape)
ypeaks[peak_pos] = peak_len*peak_width
fig, ax = plt.subplots(figsize=(15, 4))
ax.plot(x, ypeaks)
Out[8]:
In [9]:
# find peaks using the first ridge position, last ridge position as well using find_peaks
peaks = signal.find_peaks_cwt(y, wavelet=signal.ricker, widths=widths)
peaks_2 = peak_pos[np.all(((peak_width > 0), (peak_len > 5)), axis=0)]
In [10]:
print peaks, peaks_2
For now use scipy.signal.find_peaks_cwt(), compare with my own implementation
In [11]:
fig, ax = plt.subplots(24, figsize=(10,10))
for w in range(3):
for l in range(2, 10):
a = ax[w*8 + (l-2)]
peaks = peak_pos[np.all(((peak_width > w), (peak_len > l)), axis=0)]
a.plot(x,y)
a.plot(x[peaks], y[peaks], 'rx', label='w%i, l%i' % (w,l))
#a.legend()
In [12]:
# find peaks using the first ridge position, last ridge position as well using find_peaks
peaks = signal.find_peaks_cwt(y, wavelet=signal.ricker, widths=widths)
peaks_2 = peak_pos[np.all(((peak_width > 1), (peak_len > 5)), axis=0)]
fig, ax = plt.subplots(figsize=(15,5))
ax.semilogy(x,y)
ax.semilogy(x[peaks], y[peaks], 'kv', alpha=0.8)
ax.semilogy(x[peaks_2], y[peaks_2], 'rd', alpha=0.8, label='filterd width')
#ax.plot(x[peaks_3], y[peaks_3], 'bx', label='filterd length')
ax.set_ylim(200000,600000)
ax.legend()
Out[12]:
In [13]:
# find peaks using the first ridge position, last ridge position as well using find_peaks
peaks = signal.find_peaks_cwt(y, wavelet=signal.ricker, widths=widths)
peaks_2 = peak_pos[np.all(((peak_width > 5), (peak_len > 20)), axis=0)]
fig, ax = plt.subplots(figsize=(15,5))
ax.plot(x,y)
ax.plot(x[peaks], y[peaks], 'kv', alpha=0.8, label='scipy')
ax.plot(x[peaks_2], y[peaks_2], 'rd', alpha=0.8, label='filterd length and width')
#ax.plot(x[peaks_3], y[peaks_3], 'bx', label='filterd length')
ax.set_ylim(200000,520000)
ax.legend()
Out[13]:
Procedure from Zhang et al.
In [14]:
# analyze the ricker wavelet to help build the ricker wavelet
points = 100
for a in range(2, 11, 2):
wave = signal.ricker(points, a)
plt.plot(wave)
In [15]:
# note, all integrate to 0
In [16]:
# make a haar mother wavelet
def haar2(points, a):
"""
Returns a haar wavelet mother wavelet
1 if 0 <= t < 1/2
h(t) = -1 if 1/2 <= t < 1
0 otherwise`
Numpy version, not accurate right now
"""
x = np.arange(0, points) - (points - 1.0) / 2
wave = np.zeros(x.shape)
amp = 2/a
wave[np.where(np.logical_and(0 <= x, x < 0.5*a))[0]] = 1
wave[np.where(np.logical_and(-0.5*a <= x, x < 1))[0]] = -1
return wave*amp
In [17]:
# make a haar mother wavelet
def haar(points, a):
"""
Returns a haar wavelet mother wavelet
1 if 0 <= t < 1/2
h(t) = -1 if 1/2 <= t < 1
0 otherwise`
"""
vec = np.arange(0, points) - (points - 1.0) / 2
wave = np.zeros(vec.shape)
amp = 2/a
for i, x in enumerate(vec):
if 0 <= x < 0.5*a:
wave[i] = 1
elif -0.5*a <= x < 1:
wave[i] = -1
return wave*amp
In [18]:
points = 100
for a in range(2, 11, 2):
wave = haar(points, a)
plt.step(np.arange(points), wave)
In [19]:
hw = signal.cwt(y, haar, widths=widths)
plt.imshow(hw, aspect='auto', cmap='PRGn')
Out[19]:
In [20]:
ahw = np.abs(hw)
plt.imshow(ahw, aspect='auto', cmap='PRGn')
Out[20]:
Search for local minima in in the row corresponding to the peak's scale, within 3x peak scale or peak index
In [21]:
for p in peak_pos:
print p