Filtering signals

After trying out the SciPy signal module and IIR/FIR filter design and all that jazz (and finding it mind-numbing and mostly useless), I've regressed back to just multiplying signals by various frequency domain transfer functions. Here is a set of useful transfer functions for filtering signals and looking at power spectra. The examples are for filtered white gaussian noise.

	# Filtering utilities
	import numpy as np
	from numpy import pi
	from numpy import fft
	from scipy import signal
	import scipy.signal as sig
	import pylab as plt
	#### Functions #################
	def theta(x):
	if x >=0: return 1
	else: return 0
	vtheta = np.vectorize(theta)
	def lowPass(s,tau):
	return 1./(1 + 1j*2*pi*s*tau)
	def highPass(s,tau):
	return (1 + 1j*2*pi*s*tau)
	def bandPass(s,f0=1.,deltaf=1.):
	return deltaf/(2*pi) / ( (deltaf/2)**2 + (s-f0)**2 )
	def allPass(s,tau):
	return 1.
	## Time array
	N = int(1e5)
	dt = 1e-6
	T = N*dt
	t = np.arange(N)*dt
	f = fft.fftfreq(N,d=dt)
	f = fft.fftshift(f)
	## Signal
	#x = np.sin(2*pi*5e4*t)
	x = np.random.normal(0,1.0,(N))
	X = fft.fft(x)/np.sqrt(T)
	X = fft.fftshift(X)
	## Filtered signal
	#Y = X * lowPass(f,1e-5)
	Y = X * bandPass(f,f0=1e5,deltaf=1e3)
	y = fft.ifft( fft.ifftshift(Y) ) * np.sqrt(T)
	#### Plots #####################
	plt.figure(figsize=(8,8))
	plt.subplot(311) # time traces
	plt.plot(t,x, alpha=0.5, label="Raw")
	plt.plot(t,y, label="Filtered")
	plt.ylabel("$x$")
	plt.xlabel("$t$")
	plt.legend(loc='upper right')
	plt.subplot(312)
	plt.xlabel("$f$")
	plt.ylabel("$|Y_T(f)|^2$") # spectra
	#plt.plot(f,np.abs(X)**2, alpha=0.5, label="Raw")
	plt.plot(f,np.abs(Y)**2, 'g', label="Filtered")
	plt.subplot(313)
	plt.xlabel("$f$")
	plt.ylabel("$W_{YY}(f)$") # single sided spectral density of filtered signal
	plt.loglog(f,2 * np.abs(Y)**2,'g')
	plt.grid(axis='both')
	plt.tight_layout()

view raw signalFilter.py hosted with ❤ by GitHub

 

Bandpass

 Lowpass

Cameras and Python

Fun with OpenCV, to acquire video feeds and process them.

	# Display webcam image, plus plasma diagnostics
	# version 2, 2013-09-25
	# Amar
	# Changelog:
	# v2: Very slight modification of http://matplotlib.org/examples/animation/dynamic_image.html
	import numpy as np
	import cv2
	import time
	import matplotlib.animation as animation
	import matplotlib.pyplot as plt
	import matplotlib.cm as cm
	try: import Tkinter as Tk
	except ImportError: import tkinter as Tk
	from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg, NavigationToolbar2TkAgg
	from matplotlib.figure import Figure
	class Camera:
	def __init__(self,channel=0):
	self.capture = cv2.VideoCapture(channel)
	self.success,self.image = self.capture.read()
	print "Beginning acquisition ...",
	def acquire(self):
	self.success,self.image = self.capture.read()
	def close(self):
	if self.capture.isOpened(): self.capture.release()
	print "released camera"
	camera = Camera()
	fig = plt.figure()
	fig.canvas.set_window_title('USB Camera')
	im = plt.imshow(camera.image[:,:,0])
	def updatefig(*args):
	camera.acquire()
	im.set_array(camera.image[:,:,0])
	return im,
	ani = animation.FuncAnimation(fig, updatefig, interval=80, blit=True)
	plt.show()
	camera.close()

view raw camera.py hosted with ❤ by GitHub

Even zetas are sometimes rational

This calculation uses the amazing abilities of SymPy, to evaluate the even-valued zeta functions using Parseval's theorem applied to the functions $t^{s}$.

	# Evaluating zeta functions using Parseval sums
	# Amar
	# version 1, Aug 2013
	from sympy import *
	from sympy import init_printing
	init_printing()
	n = symbols('n',integer=True)
	s = symbols('s',rational=True)
	t = symbols('t')
	def a0(s):
	# DC component of fourier series
	return ((2*pi)**s)/(s+1)
	def b(n,s):
	return integrate(t**s * E**(I*n*t),(t,0,2*pi))/(2*pi)
	def a(n,s):
	# fourier components of the function t**s, defined recursively
	if s!=0 and s > Rational(1,2): return (s)/(I*n) * ( a0(s-1) - a(n,s-1) )
	elif s==Rational(1,2): return b(n,Rational(1,2))
	else: return 0
	def lhs(n,s):
	# fourier side of Parseval sum. Factor of 2 is because only single-sided fourier series is considered
	return apart( 2 * Abs(a(n,s))**2, n ) #,full=True )
	def rhs(s):
	# time side of Parseval sum
	return a0(2*s) - a0(s)**2
	def equation(n,s):
	return lhs(n,s), rhs(s)
	pprint(equation(n,1)) # gives the series sum for zeta(2)
	pprint(equation(n,2)) # gives the series sum for zeta(4)

view raw zeta.py hosted with ❤ by GitHub

The last two lines evaluate to
$$2 \sum_{n=1}^\infty \frac{1}{n^2} = \frac{\pi^2}{3},$$
$$\pi^2 \sum_{n=1}^\infty \frac{1}{n^2} + \sum_{n=1}^\infty \frac{1}{n^4} = \frac{8 \pi^4}{45}.$$

Incidentally, this recursive definition provides an easy proof of why all the even $\zeta$s are rational multiples of $\pi^{2s}$.

 I wonder if the polynomials in $t$ that evaluate to $\zeta(s)$ are special somehow ...