#
# The code to calculate the Savitzky-Golay filter coefficients
# is a shameless copy of the sg_filter.py module from Uwe Schmitt
# available from http://public.procoders.net/sg_filter
#
# Therefore PyMca author(s) do not claim any ownership of that code
# and are very grateful to Uwe for making his code available to the
# community.
#
# Copyright (C) 2008 Uwe Schmitt
#
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation
# files (the "Software"), to deal in the Software without
# restriction, including without limitation the rights to use,
# copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following
# conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
# OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
# WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
# OTHER DEALINGS IN THE SOFTWARE.
#
__author__ = "Uwe Schitt"
__copyright__ = "Uwe Schmitt"
__license__ = "MIT"
import numpy
from numpy.linalg import solve
ODD_SIGN = 1.0
[docs]def calc_coeff(num_points, pol_degree, diff_order=0):
""" calculates filter coefficients for symmetric savitzky-golay filter.
see: http://www.nrbook.com/a/bookcpdf/c14-8.pdf
num_points means that 2*num_points+1 values contribute to the
smoother.
pol_degree is degree of fitting polynomial
diff_order is degree of implicit differentiation.
0 means that filter results in smoothing of function
1 means that filter results in smoothing the first
derivative of function.
and so on ...
"""
# setup interpolation matrix
# ... you might use other interpolation points
# and maybe other functions than monomials ....
x = numpy.arange(-num_points, num_points+1, dtype=numpy.int)
monom = lambda x, deg : pow(x, deg)
A = numpy.zeros((2*num_points+1, pol_degree+1), numpy.float)
for i in range(2*num_points+1):
for j in range(pol_degree+1):
A[i,j] = monom(x[i], j)
# calculate diff_order-th row of inv(A^T A)
ATA = numpy.dot(A.transpose(), A)
rhs = numpy.zeros((pol_degree+1,), numpy.float)
rhs[diff_order] = 1
wvec = solve(ATA, rhs)
# calculate filter-coefficients
coeff = numpy.dot(A, wvec)
if (ODD_SIGN < 0) and (diff_order %2):
coeff *= ODD_SIGN
return coeff
[docs]def smooth(signal, coeff):
""" applies coefficients calculated by calc_coeff()
to signal """
N = numpy.size(coeff-1)/2
res = numpy.convolve(signal, coeff)
return res[N:-N]
[docs]def getSavitzkyGolay(spectrum, npoints=3, degree=1, order=0):
coeff = calc_coeff(npoints, degree, order)
N = numpy.size(coeff-1)/2
if order < 1:
result = 1.0 * spectrum
else:
result = 0.0 * spectrum
result[N:-N] = numpy.convolve(spectrum, coeff, mode='valid')
return result
[docs]def replaceStackWithSavitzkyGolay(stack, npoints=3, degree=1, order=0):
coeff = calc_coeff(npoints, degree, order)
N = numpy.size(coeff-1)/2
convolve = numpy.convolve
mcaIndex = -1
if hasattr(stack, "info") and hasattr(stack, "data"):
actualData = stack.data
mcaIndex = stack.info.get('McaIndex', -1)
else:
actualData = stack
if not isinstance(actualData, numpy.ndarray):
raise TypeError("This Plugin only supports numpy arrays")
# take a view
data = actualData[:]
oldShape = data.shape
if mcaIndex in [-1, len(data.shape)-1]:
data.shape = -1, oldShape[-1]
for i in range(data.shape[0]):
data[i,N:-N] = convolve(data[i,:],coeff, mode='valid')
if order > 0:
data[i, :N] = data[i, N]
data[i, -N:] = data[i,-(N+1)]
data.shape = oldShape
elif mcaIndex == 0:
data.shape = oldShape[0], -1
for i in range(data.shape[-1]):
data[N:-N, i] = convolve(data[:, i],coeff, mode='valid')
if order > 0:
data[:N, i] = data[N, i]
data[-N:, i] = data[-(N+1), i]
data.shape = oldShape
else:
raise ValueError("Invalid 1D index %d" % mcaIndex)
return
if getSavitzkyGolay(10*numpy.arange(10.), npoints=3, degree=1,order=1)[5] < 0:
ODD_SIGN = -1
if __name__ == "__main__":
x=numpy.arange(100.)
y=100*x
print("Testing first derivative")
yPrime=getSavitzkyGolay(y, npoints=3, degree=1,order=1)
if abs(yPrime[50]-100.) > 1.0e-5:
print("ERROR, got %f instead of 100." % yPrime[50])
else:
print("OK")
print("Testing second derivative")
y=100*x*x
yPrime=getSavitzkyGolay(y, npoints=3, degree=2,order=2)
if abs(yPrime[50]-100.) > 1.0e-5:
print("ERROR, got %f instead of 100." % yPrime[50])
else:
print("OK")
print("Testing third order derivative")
y=100*x*x*x
yPrime=getSavitzkyGolay(y, npoints=5, degree=3,order=3)
if abs(yPrime[50]-100.) > 1.0e-5:
print("ERROR, got %f instead of 100." % yPrime[50])
else:
print("OK")