#!/usr/bin/env python
'''
-a library of helping functions for spotting force spectroscopy peaks
+peakspot.py
-(c)Fabrizio Benedetti and Massimo Sandal , 2007
+A library of helping functions for spotting force spectroscopy peaks.
+
+Copyright 2007 by Fabrizio Benedetti and Massimo Sandal
+
+This program is released under the GNU General Public License version 2.
'''
-import numpy as np
+
+from numpy import mean
def conv_dx(data,vect):
'''
dim=len(data)
window=len(vect)
temparr=[0.0]*dim
-
+
end=dim-window
for j in range(end):
- for k in range(window):
- temparr[j]+=data[j+k]*vect[k]
+ for k in range(window):
+ temparr[j]+=data[j+k]*vect[k]
+
+ return temparr
- return temparr
-
def absdev(arr):
'''
Calculates the absolute deviation of a vector
'''
med=0.0
absD=0.0
-
- med=np.mean(arr)
+
+ med=mean(arr)
for j in arr:
absD+=abs(j-med)
return absD/len(arr)
-
+
def noise_absdev(data,positive=False,maxcut=0.2,stable=0.005):
'''
Returns the standard deviation of the noise.
The logic is: we cut the most negative (or positive) data points until the absolute deviation
becomes stable (it doesn't vary more than 0.005) or we have cut more than maxcut*len(data) points.
Then calculate the absolute deviation.
-
+
If positive=True we cut the most positive data points, False=we cut the negative ones.
'''
+ #TOD: check if this is necessary
out=[item for item in data] #we copy just to be sure...
out.sort()
if positive:
out.reverse()
-
+
temp_absdev=absdev(out)
-
+
for index in range(len(out)):
cutindex=(index+1)*5
cut_absdev=absdev(out[cutindex:]) #we jump five points after five points...
temp_absdev=cut_absdev
else:
break
-
+
return cut_absdev
-
+
def abovenoise(convoluted,noise_level,cut_index=0,abs_devs=4):
'''
Generates a vector which is 0 where the vector is less than abs_devs*noise_level ; 1 if not (spike).
'''
#calculate absolute noise deviation
#noise_level=noise_absdev(convoluted[cut_index:])
-
+
above=[]
-
+
for index in range(len(convoluted)):
if index<cut_index:
above.append(0)
else:
above.append(0)
return above
-
+
def find_peaks(above, seedouble=10):
'''
Finds individual peak location.
abovenoise() finds all points above a given threshold in the convolution. This point is often not unique
but there is a "cluster" of points above a threshold.
- Here we obtain only the coordinates of the largest *convolution* spike for each cluster.
-
+ Here we obtain only the coordinates of the largest *convolution* spike for each cluster.
+
above=vector obtained from abovenoise()
seedouble=value at which we want to "delete" double peaks. That is, if two peaks have a distance
< than $seedouble points , only the first is kept.
nonzero=[]
peaks_location=[]
peaks_size=[]
-
+
for index in range(len(above)):
if above[index] != 0:
nonzero.append(index)
nonzero=[]
else:
pass
-
+
#recursively eliminate double peaks
#(maybe not the smartest of cycles, but i'm asleep...)
temp_location=None
if peaks_location[index+1]-peaks_location[index] < seedouble:
temp_location=peaks_location[:index]+peaks_location[index+1:]
if temp_location != []:
- peaks_location=temp_location
-
+ peaks_location=temp_location
+
return peaks_location,peaks_size
\ No newline at end of file
projects / hooke.git / blobdiff