Fun with signal processing and ... matplotlib

Well, I remember overhearing a conversation one day in a  laboratory once 10 years ago.

Electronic engineers were proud to have helped their team to higher CDD resolution with a trick they published:
they would use a median filter to filter out the photon arriving from the CDD from the thermal noise. They were even very proud of seeing a very signal far smaller than the noise.

I thought to myself. Such a scam, this is to stupid it cannot work.

So I decided to play with matplotlib to test it.

First we generate white noise, then given a probability of 1/n a photon is seen, and its amplitude is 1/A of the signal.

Then, I used #python-fr to have advices on the best way to do a moving average (which is the usual filter used for filtering out signals), then I tried the median filter.

A median filter is just as silly as taking the median in a moving window on a time serie.

It cant work, can't it?

Well, it does.

Is it not hard  to code, but beware your eyes will cry from PEP8 insanity and mapltotlib style mixed with pythonic style (call it art):
(https://gist.github.com/3794506)

	#!/usr/bin/python
	# -*- coding: utf-8 -*-
	import matplotlib.pyplot as plt
	from scipy.signal import medfilt
	import numpy as np
	NOISE_AMPLITUDE = 4.0
	PERIOD = 200
	PRESENCE_PROBABILITY = .4
	OPT_MA_WINDOW = 75
	OPT_MEDIAN_WINDOW = 41
	# size of sample
	S_O_S = 4000
	def canonical_ma( a_signal, size_of_window):
	""" method used canoncically in signal processing to obtain
	moving average the fastest way possible.
	In plain english: it is the average of the nterms around a
	point in a serie
	"""
	window = np.ones( size_of_window )
	return np.roll(
	np.convolve( window / size_of_window , a_signal, 'valid' ) ,
	size_of_window / 2
	)
	#contributed by generic_genus (reddit)
	signal = np.where(np.arange(4000) % PERIOD > PERIOD / 2, 0.5, -0.5)
	noise = np.random.uniform(
	low=-0.5*NOISE_AMPLITUDE,
	high=0.5*NOISE_AMPLITUDE,
	size=S_O_S
	)
	noise = np.where(
	np.random.uniform(size=len(noise)) < PRESENCE_PROBABILITY,
	signal, noise
	)
	fig=plt.figure(figsize=(12,10))
	plt.suptitle(
	"Square signal mixed in white noise filtering with P=%f proba and NS ratio=%f" %
	( PRESENCE_PROBABILITY,NOISE_AMPLITUDE)
	)
	tint = 255
	ax=plt.subplot(211)
	ax.plot( noise, color='blue', label="signal mixed in random noise")
	# amer45 (reddit) made me realize I should use existing available best way to do it
	ax.plot(
	medfilt(noise,[OPT_MEDIAN_WINDOW]),
	color ='#%x0000' % tint, label="median filter"
	)
	ax.plot( signal , color ='yellow', label="initial signal" )
	ax.legend()
	ax=plt.subplot(212)
	ax.plot( noise, color='blue', label="signal mixed in random noise")
	ax.plot( signal , color ='yellow', label="initial signal" )
	ax.plot(
	canonical_ma(noise,OPT_MA_WINDOW),
	color ='#00%x00' % tint, label="moving average"
	)
	ax.legend()
	plt.show()

view raw median_vs_arma.py hosted with ❤ by GitHub

Moral of the story? 

 

Well, I should have to publish the code that shows why 40 for median and 75 for movering average are the optimal settings, and how dazzling fast matplotlib is.

But the real point is python + matplotlib are amazing. And that sometimes simple ideas give amazing results.

Have fun!