00001 """
00002 cookb_signalsmooth.py
00003
00004 from: http://scipy.org/Cookbook/SignalSmooth
00005 """
00006
00007 import numpy as np
00008 import matplotlib.pyplot as plt
00009
00010 def smooth(x, window_len=10, window='hanning'):
00011 """smooth the data using a window with requested size.
00012
00013 This method is based on the convolution of a scaled window with the signal.
00014 The signal is prepared by introducing reflected copies of the signal
00015 (with the window size) in both ends so that transient parts are minimized
00016 in the begining and end part of the output signal.
00017
00018 input:
00019 x: the input signal
00020 window_len: the dimension of the smoothing window
00021 window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
00022 flat window will produce a moving average smoothing.
00023
00024 output:
00025 the smoothed signal
00026
00027 example:
00028
00029 import numpy as np
00030 t = np.linspace(-2,2,0.1)
00031 x = np.sin(t)+np.random.randn(len(t))*0.1
00032 y = smooth(x)
00033
00034 see also:
00035
00036 numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
00037 scipy.signal.lfilter
00038
00039 TODO: the window parameter could be the window itself if an array instead of a string
00040 """
00041
00042 if x.ndim != 1:
00043 raise ValueError, "smooth only accepts 1 dimension arrays."
00044
00045 if x.size < window_len:
00046 raise ValueError, "Input vector needs to be bigger than window size."
00047
00048 if window_len < 3:
00049 return x
00050
00051 if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
00052 raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
00053
00054
00055
00056
00057
00058 s = x
00059
00060 if window == 'flat':
00061 w = np.ones(window_len,'d')
00062 else:
00063 w = getattr(np, window)(window_len)
00064 y = np.convolve(w/w.sum(), s, mode='same')
00065 return y[window_len-1:-window_len+1]
00066
00067
00068
00069
00070 from scipy import signal
00071
00072 def gauss_kern(size, sizey=None):
00073 """ Returns a normalized 2D gauss kernel array for convolutions """
00074 size = int(size)
00075 if not sizey:
00076 sizey = size
00077 else:
00078 sizey = int(sizey)
00079 x, y = np.mgrid[-size:size+1, -sizey:sizey+1]
00080 g = np.exp(-(x**2/float(size) + y**2/float(sizey)))
00081 return g / g.sum()
00082
00083 def blur_image(im, n, ny=None) :
00084 """ blurs the image by convolving with a gaussian kernel of typical
00085 size n. The optional keyword argument ny allows for a different
00086 size in the y direction.
00087 """
00088 g = gauss_kern(n, sizey=ny)
00089 improc = signal.convolve(im, g, mode='valid')
00090 return(improc)
00091
00092
00093 def smooth_demo():
00094 t = np.linspace(-4,4,100)
00095 x = np.sin(t)
00096 xn = x + np.random.randn(len(t)) * 0.1
00097 y = smooth(x)
00098 ws = 31
00099
00100 plt.subplot(211)
00101 plt.plot(np.ones(ws))
00102
00103 windows=['flat', 'hanning', 'hamming', 'bartlett', 'blackman']
00104
00105 plt.hold(True)
00106 for w in windows[1:]:
00107
00108 plt.plot(getattr(np, w)(ws))
00109
00110 plt.axis([0,30,0,1.1])
00111
00112 plt.legend(windows)
00113 plt.title("The smoothing windows")
00114 plt.subplot(212)
00115 plt.plot(x)
00116 plt.plot(xn)
00117 for w in windows:
00118 plt.plot(smooth(xn, 10, w))
00119 l = ['original signal', 'signal with noise']
00120 l.extend(windows)
00121 plt.legend(l)
00122 plt.title("Smoothing a noisy signal")
00123
00124
00125
00126 if __name__=='__main__':
00127
00128
00129 smooth_demo()
00130
00131
00132 X, Y = np.mgrid[-70:70, -70:70]
00133 Z = np.cos((X**2+Y**2)/200.)+ np.random.normal(size=X.shape)
00134 Z2 = blur_image(Z, 3)
00135 plt.figure()
00136 plt.imshow(Z)
00137 plt.figure()
00138 plt.imshow(Z2)
00139 plt.show()
00140