from pylab import *

import os
def play(data,r=1,pad=8000):
    data = array(data*16000.0/amax(abs(data)),'int16')
    if r>1: data = array(list(data)*r)
    data = r_[zeros(pad,data.dtype),data,zeros(pad,data.dtype)]
    data.tofile("_out.raw")
    os.system("play -s2L -r 8000 _out.raw")

x = 2*pi*arange(1024,dtype='d')/1024.0

plot(x,sin(x))
plot(x,cos(x))

data128 = sin(128*x)
data128 /= norm(data128)
plot(data128)
play(data128,r=10)

data64 = sin(64*x)
data64 /= norm(data64)
plot(data64)
play(data64,r=10)

data32 = sin(32*x)
data32 /= norm(data32)
plot(data32)
play(data32,r=10)

data32c = cos(32*x)
data32c /= norm(data32c)
plot(data32c)
play(data32c,r=10)

dataodd = sin(31.3*x)
plot(dataodd)
play(dataodd,r=10)

dot(data32,data32)

dot(data32,data32c)

dot(data32,data64)

dot(data32,dataodd)

data32p = sin(32*x+0.3)
plot(data32p)
play(data32p,r=10)

dot(data32,data32p)

dot(data32c,data32p)

dot(data32,data64)

rec = dot(data32,data32p)*data32 + dot(data32c,data32p)*data32c
print norm(rec-data32p)
plot(rec)
plot(0.9*data32p)

sins = [sin(i*x) for i in range(1,513)]

figsize(12,8)
plot(sins[0])
plot(sins[1])
plot(sins[2])

coss = [cos(2.0*pi*arange(1024)*i/1024.0) for i in range(0,512)]

figsize(12,8)
plot(coss[0])
plot(coss[1])
plot(coss[2])

ylim(-0.1,1.1)
plot(sin(0*x)[-10:])
plot(cos(0*x)[-10:])

plot(sin(512*x)[-10:])
plot(cos(512*x)[-10:])

sins = array(sins)
coss = array(coss)

figsize(24,12)
gray()
subplot(131); imshow(dot(sins,sins.T),vmax=1); print amin(dot(sins,sins.T)),amax(dot(sins,sins.T))
subplot(132); imshow(dot(coss,coss.T),vmax=1); print amin(dot(coss,coss.T)),amax(dot(coss,coss.T))
subplot(133); imshow(dot(sins,coss.T),vmax=1); print amin(dot(sins,coss.T)),amax(dot(sins,coss.T))

figsize(12,12)
subplot(121); imshow(sins[:100,:100])

figsize(10,5)
plot(dot(sins,data32))
plot(dot(coss,data32))

figsize(10,5)
plot(dot(sins,data32p)[:50])
plot(dot(coss,data32p)[:50])

figsize(10,5)
plot(dot(sins,data64)[:-1]**2+dot(coss,data64)[1:]**2)

ft = r_[sins,coss]
ift = inv(ft)
def transform(signal):
    return dot(ft,signal)
def itransform(signal):
    return dot(ift,signal)

def spectrum(signal):
    return dot(sins,signal)[:-1]**2+dot(coss,signal)[1:]**2

from scipy.ndimage import filters
wave = filters.gaussian_filter(randn(1024),5.0)
plot(wave)

plot(spectrum(wave))

plot(transform(wave))

norm(wave-itransform(transform(wave)))

for phase in [0.0,pi/2,0.3,0.7]:
    frequency = 100
    plot(transform(sin(2.0*pi*arange(1024)*frequency/1024.0+phase)))

for phase in [0.0,pi/2,0.3,0.7]:
    frequency = 100
    plot(transform(sin(2.0*pi*arange(1024)*frequency/1024.0+phase))[95:105])

for phase in [0.0,pi/2,0.3,0.7]:
    frequency = 100
    signal = sin(frequency*x+phase)
    plot(spectrum(signal))

for phase in [0.0,pi/2,0.3,0.7]:
    frequency = 100
    signal = sin(frequency*x+phase)
    plot(spectrum(signal)[95:105])

for frequency in [1,10,100,500]:
    signal = sin(frequency*x)
    plot(transform(signal))

def b(k,n,N):
    return exp(-2*pi*1j*k*n/N)
N = 420
n = arange(0,420)

s = b(3,n,N)

ylim(-1.1,1.1)
plot(real(s))
plot(imag(s))
plot(abs(s),color='r')

F = array([b(k,n,N) for k in range(0,N)])

a1 = 1.5
a2 = 2.3
omega1 = 6*2*pi/N
omega2 = 14*2*pi/N
y = a1*sin(n*omega1) + a2*sin(n*omega2)
plot(y)

Y = dot(F,y)
plot(abs(Y)**2)

plot(Y) # ???

[i for i in range(N) if abs(Y[i])>100]

def max1(x): return x*1.0/amax(abs(x))

signal = n
transform = dot(F,signal)
plot(max1(signal))
plot(max1(abs(transform)))

plot(roll(max1(signal),200))

signal = concatenate([n,n])[2*n]
transform = dot(F,signal)
plot(max1(signal))
plot(max1(abs(transform)))

signal = (n==100)
transform = dot(F,signal)
ylim(-0.1,1.1)
plot(max1(signal))
plot(max1(abs(transform)))

ylim(-0.1,1.1)
plot(max1(concatenate([signal for i in range(8)])))

signal = (n>=100) * (n<102)
transform = dot(F,signal)
plot(max1(signal))
plot(max1(abs(transform)))

signal = (n>=100) * (n<110)
transform = dot(F,signal)
plot(max1(signal))
plot(max1(abs(transform)))

signal = exp(-(n-200)**2/100.0)
signal /= amax(signal)
transform = dot(F,signal)
plot(max1(signal))
plot(max1(abs(transform)))

plot(max1(signal))
plot(roll(max1(abs(transform)),200))

signal = (n>=100) * (n<110)
transform = dot(F,signal)
reconstructed = dot(inv(F),transform)

plot(max1(signal))
plot(max1(abs(transform)))
plot(max1(reconstructed))

double = dot(F,transform)
plot(max1(double))

imshow(abs(dot(F,F))); gray()

figsize(10,10); xticks([]); yticks([])
imshow(imread("fft-image.png"))

def myfft(s): 
    N = len(s)
    if N<=1: return s
    Se = fft(s[::2])
    So = fft(s[1::2])
    e = exp(arange(0,N/2)*2j*pi/N)
    return r_[Se+e*So,Se-e*So]
    
myfft(array([1.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0]))

k = arange(256)
signal = 1.0*(k>=100)*(k<110)
transform = myfft(signal)
plot(max1(signal))
plot(max1(abs(transform)))

mix = data32+data128
play(mix,r=10)

plot(abs(fft(mix)))

plot(abs(fft(mix))[:200])

plot(abs(fft(mix))[120:140])

transform = fft(mix)

transform[120:130] = 0
transform[-130:-120] = 0
plot(abs(transform))

plot(real(ifft(transform))[:100])
plot(mix[:100])

play(mix,r=10)

play(real(ifft(transform)),r=10)

test = numpy.fromfile("test.raw","int32")[::2]
test = test*1.0/amax(abs(test))
play(test)
print amax(test)

noise = test+sin(arange(len(test))/10.0)
play(noise)

Noise = fft(noise)

plot(abs(Noise)**.5)

plot((abs(Noise)**.5)[400:500])

Noise[480:500] = 0
Noise[-500:-480] = 0
plot(abs(Noise)**.5)

denoised = ifft(Noise)
play(real(denoised))

play(noise)