%pylab inline

phs = linspace(0, 20* 2*pi, 1024, endpoint=False)
g = cos(phs + 0.4)
plot(g)

from scipy.signal import freqz

plot(abs(fft.rfft(g[:512], n=512)))
xlim(0, 400);

plot(abs(fft.rfft(g[512:1024], n=512)))
xlim(0, 400);

angle(fft.rfft(g[:512])[10])

angle(fft.rfft(g[512:1024])[10])

0.4 + (10 * 2 * pi)

(0.4 + (10 * 2 * pi)) %(2* pi)

(10.1 * 2 * pi)

(10.1 * 2 * pi)%(2* pi)

phs = linspace(0, 23.4* 2*pi, 1024, endpoint=False)
g = cos(phs)
plot(g)

plot(abs(fft.rfft(g[:512], n=512)))
xlim(0, 50);

X = abs(fft.rfft(g[:512], n=512))
argmax(X)

plot(g[:512])
plot(g[512:1024])

angle(fft.rfft(g[:512]))[12]

angle(fft.rfft(g[512:1024]))[12]

angle(fft.rfft(g[512:1024]))[12] - angle(fft.rfft(g[:512]))[12]

12 * 2 * pi

12 * 2 * pi + __

_ / (2 * pi)

_ * 2

f0 = 203.3
p0 = 0.3
N = 4096

test_phs = linspace(0, (f0 * 2*pi), N, endpoint=False)
test_sig = sin(test_phs + p0) + 0.5 * cos(p0 + (test_phs * 4))

X = rfft(test_sig[:512])
plot(abs(X))
# xlim(0,1000);

X = rfft(test_sig[:512])
plot(abs(X))
xlim(10, 40);

argmax(abs(X))

(argmax(abs(X)) + 1)

hopsize = 128
winsize = 512
win_start = linspace(0, len(test_sig) - winsize, hopsize)

stft = []

for start in win_start:
    X = rfft(test_sig[start: start + winsize]* hanning(winsize))
    stft.append(X)

argmax(abs(stft[1]))

plot(abs(stft[1]))
plot(abs(stft[5]))
plot(abs(stft[9]))
xlim(10, 40)

imshow(angle(stft).T, aspect='auto', cmap=cm.gray, interpolation='nearest')
gcf().set_figwidth(10)
colorbar();

plot(angle(stft).T[24])

plot(angle(stft).T[25])

unwrapped = [angle(stft).T[25][0]]

for p in angle(stft).T[25][1:]:
    while (p < unwrapped[-1]):
        p += 2 * pi
    unwrapped.append(p)

plot(unwrapped, 'x')

delta_phs = diff(unwrapped)

plot(delta_phs)

k = 25 # bin number
num_cycles = k * (hopsize/ float(winsize)) # how many cycles of bin between STFT hops

total_phs_k = num_cycles * 2* pi # Total accumulated phase for bin
num_cycles, total_phs_k

phs_change = delta_phs[0]
phs_change

total_phs_k % (2 * pi)

div = phs_change - total_phs_k % (2 * pi)
div

final_phs = total_phs_k + div
final_phs

final_phs/(2*pi)

(winsize/float(hopsize)) *final_phs/(2*pi)

(N/float(winsize))*(winsize/hopsize) *final_phs/(2*pi)

(N/float(winsize))*25

(N/float(winsize))*26