returns the difference matrix

derivativeMatrix.py

import numpy
from scipy.linalg import toeplitz


def pascal_row(n):
  row = [1]
  for col in range(1, n): row.append(row[-1] * (n - col) / col)
  return numpy.array(row)
  


def computeMatrix(signalLength, derivativeOrder, delta = 1):
  if derivativeOrder % 2 == 0:
    coef_signs = numpy.empty((derivativeOrder+1,))
    coef_signs[::2] = 1
    coef_signs[1::2] = -1
    row = pascal_row(derivativeOrder+1)*coef_signs
  else:
    coef_signs = numpy.empty((derivativeOrder+1,))
    coef_signs[::2] = -1
    coef_signs[1::2] = 1
    row = pascal_row(derivativeOrder+1)*coef_signs       
  
  D =  toeplitz(numpy.hstack([row[0], numpy.zeros([signalLength - derivativeOrder - 1,])]), numpy.hstack([row, numpy.zeros([signalLength - derivativeOrder - 1,])]))
  
  return (1.0/delta**derivativeOrder)*D
  

def example():
  from pylab import *
  delta = .01
  x_points = numpy.arange(0,5,delta)
  sample_func = x_points**5
  D4 = computeMatrix(len(sample_func), 4, delta)
  
  res = numpy.dot(D4, sample_func)
  analytic_deriv = 5*4*3*2*x_points
  
  plot(x_points, sample_func, label = 'Original Function')
  plot(x_points[:-4], res, label = "Estimated Derivative using Difference Method")
  plot(x_points, analytic_deriv, label = "Analytic Derivative")
  legend(loc='upper left')
  show()
 
 
if __name__ == '__main__':
  example()