#   w a v e . p y
#
import sys
from Tkinter import *
from math    import *

class wave :
	def __init__ (self, points=400, formula=None) :
		self.data = [0.0]*points
		self.points= points
		if formula :
			for p in range(points) :
				x = p*pi*2/points
				self.data[p] = eval(formula)

	def __add__ (self, other) :
		target = wave(points=self.points)
		for i in range(self.points) :
			target.data[i] = self.data[i] + other.data[i]
		return target

	def __mul__ (self, other) :
		target = wave(points=self.points)
		if type(other) == type(5) or type(other) == type(5.0) :
			for i in range(self.points) :
				target.data[i] = self.data[i] * other
		else :
			for i in range(self.points) :
				target.data[i] = self.data[i] * other.data[i]
		return target

	def __sub__ (self, other) :
		target = wave(points=self.points)
		for i in range(self.points) :
			target.data[i] = self.data[i] - other.data[i]
		return target

	def integral(self) :
		ans = 0.0
		for pt in self.data : ans = ans+pt
		return ans*2*pi/self.points

	def plot (self, title="??", pixHeight=None, maxY=None, others=[]) :
		if not pixHeight : pixHeight = self.points*2/3   # Pleasant ratio
		pixWidth = self.points
		# find max and min data to scale
		if not maxY :
			maxY = max (max(self.data), -min(self.data))
		offset = pixHeight/2
		scale = offset/maxY

		win = Tk()
		win.title (title)
		canvas = Canvas(win,width=pixWidth,height=pixHeight)
		# create zero line
		canvas.create_line(0,offset,pixWidth,offset)
		canvas.pack()

		self.plotOne (canvas, pixWidth, scale, offset)
		for i in range(len(others)) :
			others[i].plotOne (canvas, pixWidth, scale, offset)
		if sys.platform == "win32" : win.mainloop()

	def plotOne (self, canvas, pixWidth, scale, offset) :
		for x in range(pixWidth) :
			y = offset - self.data[x] * scale
			if x : canvas.create_line(x-1,yprev,x,y)
			yprev = y

	def fft (self) :
		work = self * 1   # Harmonics will be stripped from this
		for harm in range(1,10) :
			formula="-sin(%d*x)" % harm
			area = (wave(formula=formula)*work).integral()
			amplitude = area/-pi
			if amplitude > .000001 :
				 print "Harmonic=",harm, " Amplitude=%.04f" % amplitude
			takeAway = wave(formula="sin(%d*x) * %f" % (harm,amplitude))
			work = work-takeAway

def test() :
	p1 = wave(formula="sin(x)/1")
	p3 = wave(formula="sin(3*x)/3")
	p5 = wave(formula="sin(5*x)/5")
	mys = p1+p3+p5
	mys.fft()

if __name__ == "__main__" : test()