# simple fourier analysis and synthesis of a real function f(t)
# f[] --> a[],b[] coefficients of cos and sin
# Wolfgang.Urban@schule.at, 7/2004

from Tkinter import *

########################################################################
#
def getmax(f): return float(max( 0.001,abs(max(f)),abs(min(f)) ))
    
def getmax2(a,b): return max(getmax(a),getmax(b))

########################################################################
#
class Raster:

    def __init__(self,n=128,width=650):
        root = Tk()
        root.title("Fourier: function of time")

        self.sx = width/n
        self.sy = int(self.sx*n/6.28)
        self.mid = self.sy+4
        
        self.cv = Canvas(root,width=n*self.sx+4, height=2*self.sy+4)
        self.cv.pack()

        self.cv.create_line(2,self.mid,2+self.sx*n,self.mid,fill="#000000")

        self.lines = [0]*(n+1)
        for i in range(n+1):
            self.lines[i] = self.cv.create_line(0,0,0,0,fill="#0000ff")

    def show(self,l,maxvalue=0):
        if abs(maxvalue)<0.0001: maxvalue = getmax(l)
        scale = self.sy/float(maxvalue)
        n = len(l)
        h = [0]*(n+1)
        for i in range(n): h[i] = int(round(self.mid - scale*l[i]))
        h[n] = h[0]
        
        for i in range(n):
            self.cv.coords(self.lines[i],2+i*self.sx,h[i],2+(i+1)*self.sx,h[i+1])
        

########################################################################
#
class Bars:
    def __init__(self,n=64,width=650):
        root = Tk()
        root.title("Fourier: spectre cos/sin")
        
        height = 420
        self.mida = height/4
        self.midb = height*3/4
        self.scale = float(height/4-1)
        sx=self.sx = width/n
        bx=self.bx = self.sx-4
        
        self.cv = Canvas(root,width=sx*n+2, height=height)
        self.cv.pack()

        self.cv.create_line(2,self.mida+1,0+sx*n,self.mida+1,fill="#000000",width=2)
        self.cv.create_line(2,self.midb+1,0+sx*n,self.midb+1,fill="#000000",width=2)

        self.barsa = [0]*n
        for i in range(n):
            self.barsa[i] = self.cv.create_rectangle(
                2+i*sx,self.mida,2+i*sx+bx,self.mida,fill="#0000ff")#,outline="#0000ff")
        self.barsb = [0]*n
        for i in range(n):
            self.barsb[i] = self.cv.create_rectangle(
                2+i*sx,self.midb,2+i*sx+bx,self.midb,fill="#ff0000")#,outline="#ff0000")

    def show(self,a,b,maxvalue=0):
        if abs(maxvalue)<0.0001: maxvalue = getmax2(a,b)
        sc = self.scale/float(maxvalue)
        ha = [0]*len(a)
        for i in range(len(a)): ha[i] = self.mida - int(round(sc*a[i]))
        hb = [0]*len(b)
        for i in range(len(b)): hb[i] = self.midb - int(round(sc*b[i]))
        for i in range(len(b)):
            self.cv.coords(self.barsa[i],2+i*self.sx,ha[i],2+i*self.sx+self.bx,self.mida)            
        for i in range(len(b)):
            self.cv.coords(self.barsb[i],2+i*self.sx,hb[i],2+i*self.sx+self.bx,self.midb)            
        self.cv.update()

    # don't use a[0]
    def show1(self,a,b,maxvalue=0):
        if abs(maxvalue)<0.0001: maxvalue = getmax2(a[1:],b)
        sc = self.scale/float(maxvalue)
        ha = [0]*len(a)
        for i in range(len(a)): ha[i] = self.mida - int(round(sc*a[i]))
        hb = [0]*len(b)
        for i in range(len(b)): hb[i] = self.midb - int(round(sc*b[i]))
        for i in range(1,len(b)):
            self.cv.coords(self.barsa[i],2+i*self.sx,ha[i],2+i*self.sx+self.bx,self.mida)            
        for i in range(1,len(b)):
            self.cv.coords(self.barsb[i],2+i*self.sx,hb[i],2+i*self.sx+self.bx,self.midb)            
        self.cv.update()

########################################################################
########################################################################

def analysis(f):
    from math import pi,sin,cos
    num = len(f)
    a,b = [0]*num,[0]*num

    for n in range(num):
        phi = 2.0*pi*n/num
        sum1 = sum2 = 0.0
        for k in range(num):
            sum1 += f[k]*cos(phi*k)
            sum2 += f[k]*sin(phi*k)
        a[n] = sum1/float(num)
        b[n] = sum2/float(num)
    return a,b


def synthesis(a,b):
    from math import pi,sin,cos
    num = len(a)
    f = [0]*num

    for n in range(num):
        phi = 2.0*pi*n/num
        sum = a[0] 
        for k in range(1,num):
            sum += a[k]*cos(phi*k) + b[k]*sin(phi*k)
        f[n] = sum 
        
    return f

########################################################################
########################################################################

def normalize(ff):
    ff = f[:]
    m = -1
    for y in f:
        if abs(y)>m: m=abs(y)
    if m<0.0001: m = 0.0001
    for i in range(len(f)):
        f[i] = f[i]/m
    return f

def power(a,b):
    from math import hypot
    powr = [0]*len(a)
    for i in range(len(a)):
        powr[i] = hypot(a[i],b[i])
    normalize(powr)
    return powr
        
#################################################################

# test synthesis: given a,b
# show a,b and f
# large cos with smaller sinusses
def test():
    r = Raster(128,800)
    spectre = Bars(128,800)
    a = [0]*128
    b = [0]*128
    b[1] = 1
    a[5] = 0.2
    a[10] = -0.2
    f = synthesis(a,b)
    r.show(f)
    spectre.show(a,b)

# test analysis: given f
# show a,b and f
# f has pulse at 0
def test1():
    r = Raster(64,800)
    spectre = Bars(64,800)
    f = [0]*64
    f[0] = f[1] = 1.0
    a,b = analysis(f)
    r.show(f)
    spectre.show(a,b)


# test analysis: given f
# show half of a,b (no redundant data) and f
# f is a step function
def test2():
    r = Raster(128,800)
    spectre = Bars(64,800)
    f = [0]*128
    for i in range(128):
        if i<32: f[i]=-1.0
        else: f[i] = 1.0
    a,b = analysis(f)
    r.show(f)
    spectre.show(a[:64],b[:64])


# kind of synthesis animation
# given series of f
# show half of a,b and f
def test3():
    r = Raster(64,600)
    spectre = Bars(32,600)
    f = [0]*64
    for k in range(64):
        for i in range(64):
            if i<k: f[i]=1.0
            else: f[i] = -1.0
        a,b = analysis(f)
        r.show(f)
        spectre.show(a[:32],b[:32],1)

# same as test4()
# but without constant a[0], we can use magnification
def test4():
    r = Raster(64,600)
    spectre = Bars(32,600)
    f = [0]*64
    for k in range(64):
        for i in range(64):
            if i<k: f[i]=1.0
            else: f[i] = -1.0
        a,b = analysis(f)
        r.show(f)
        spectre.show1(a[:32],b[:32],0.5)



test4()