Derivative with noise

gistfile1.py

import numpy as np
import matplotlib.pyplot as plt
import scipy.signal

# generate some example data
t = np.linspace(0, 1, 1000)
# interesting temperature profile
temperature = 10*(t-0.4)*(t-0.2)*(t-0.8) + 2.0
# with exact derivative
dtempdt = 10*((t-0.2)*(t-0.8) + (t-0.4)*(t-0.8) + (t-0.4)*(t-0.2))
# our true signal is a function of temp and its derivative
signal = 0.8*temperature+0.5*dtempdt+0.5
# noisy measurements
noisy_signal = signal + np.random.normal(0, 0.2, len(t))
noisy_temperature = temperature + np.random.normal(0, 0.2, len(t))

# design the low-pass fiter.
# 50 taps means convolving with a length-100 window
# cutoff is nyquist*0.01 (sampling*0.005)
numtaps=50
filter = scipy.signal.fir_filter_design.firwin(numtaps, 0.01)
filtered_temperature = scipy.signal.fftconvolve(noisy_temperature, filter,
        mode='same')
# plot temperature filtering results
plt.figure(1)
plt.clf()
ax=plt.subplot(211)
plt.plot(t, temperature, label="true")
plt.plot(t, noisy_temperature, '+', label="measured")
plt.plot(t[numtaps:-numtaps], filtered_temperature[numtaps:-numtaps], label="filtered")
ax.legend(loc="best")
ax.set_title("Temperature")
# can do the derivative and filtering all at once by differentiating the filter
derivative_filter = np.gradient(filter)/(t[1]-t[0])
derivative_temperature = scipy.signal.fftconvolve(noisy_temperature,
        derivative_filter, mode='same')
# plot the temperature derivative
ax=plt.subplot(212)
plt.plot(t, dtempdt, label="true")
plt.plot(t, np.gradient(noisy_temperature)/np.gradient(t), '+', label="unfiltered")
plt.plot(t[numtaps:-numtaps], derivative_temperature[numtaps:-numtaps],
        label="filtered")
ax.set_ylim(-12,12)
ax.legend(loc="best")
ax.set_title("Derivative Temperature")
ax.set_xlabel("time")
plt.draw()
plt.figure(2)
# plot the signal filtering result
plt.clf()
filtered_signal = scipy.signal.fftconvolve(noisy_signal, filter, mode='same')
ax=plt.subplot(111)
plt.plot(t, signal, label="true")
plt.plot(t, noisy_signal, '+', label="measured")
plt.plot(t[numtaps:-numtaps], filtered_signal[numtaps:-numtaps], label="filtered")
ax.legend(loc="best")
ax.set_title("Signal")
plt.draw()
plt.figure(3)
# plot the interesting things
plt.clf()
ax=plt.subplot(211)
# first, signal vs temperature
plt.plot(filtered_temperature[numtaps:-numtaps],
        filtered_signal[numtaps:-numtaps], '+', label="filtered")
plt.plot(temperature, signal, label="true")
ax.set_xlabel("temperature")
ax.set_ylabel("signal")
ax.legend(loc="best")
ax=plt.subplot(212)
# then, signal vs dtemp/dt
plt.plot(derivative_temperature[numtaps:-numtaps],
        filtered_signal[numtaps:-numtaps], '+', label="filtered")
plt.plot(dtempdt, signal, label="true")
ax.legend(loc="best")
ax.set_xlabel("temperature derivative")
ax.set_ylabel("signal")
plt.show()
plt.draw()