%pylab inline
Populating the interactive namespace from numpy and matplotlib
phs = linspace(0, 20* 2*pi, 1024, endpoint=False)
g = cos(phs + 0.4)
plot(g)
[<matplotlib.lines.Line2D at 0x7ff3aae05dd0>]
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);
Starting phase for first window:
angle(fft.rfft(g[:512])[10])
0.39999999999999808
Starting phase for second window: (no phase change!)
angle(fft.rfft(g[512:1024])[10])
0.40000000000000219
Phase at the end of the window:
0.4 + (10 * 2 * pi)
63.23185307179586
If we wrap the phase:
(0.4 + (10 * 2 * pi)) %(2* pi)
0.3999999999999986
However for a frequency that has a fractional cycle the ending window phase is:
(10.1 * 2 * pi)
63.46017160251382
And wrapped:
(10.1 * 2 * pi)%(2* pi)
0.6283185307179551
Let's see a concrete example:
phs = linspace(0, 23.4* 2*pi, 1024, endpoint=False)
g = cos(phs)
plot(g)
[<matplotlib.lines.Line2D at 0x7ff39e0b91d0>]
plot(abs(fft.rfft(g[:512], n=512)))
xlim(0, 50);
X = abs(fft.rfft(g[:512], n=512))
argmax(X)
12
The precise center is at 23.4 / 2.
Because the number of cycles is not an integer, the starting phase for the second window is different:
plot(g[:512])
plot(g[512:1024])
[<matplotlib.lines.Line2D at 0x7ff39e087b50>]
Starting phase first window:
angle(fft.rfft(g[:512]))[12]
-0.9519681448511772
Starting phase second window:
angle(fft.rfft(g[512:1024]))[12]
-2.8196061225509639
Phase difference:
angle(fft.rfft(g[512:1024]))[12] - angle(fft.rfft(g[:512]))[12]
-1.8676379776997867
The total elapsed phase for the bin should be:
12 * 2 * pi
75.39822368615503
But what actually elapsed (adding the phase difference measured) was:
12 * 2 * pi + __
73.530585708455249
Which if we count cycles in that phase accumulation we get:
_ / (2 * pi)
11.702756183942927
Must multiply by two because we split the signal into two windows and get the frequency we put in (to the presicion allowed by the sampling rate):
_ * 2
147.0611714169105
Two sine tones ( Second one 4 times the frequency and half the amplitude)
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))
25
(argmax(abs(X)) + 1)
26
The frequency we are looking for is in between. Let's do STFT with overlap:
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]))
25
plot(abs(stft[1]))
plot(abs(stft[5]))
plot(abs(stft[9]))
xlim(10, 40)
(10, 40)
imshow(angle(stft).T, aspect='auto', cmap=cm.gray, interpolation='nearest')
gcf().set_figwidth(10)
colorbar();
plot(angle(stft).T[24])
[<matplotlib.lines.Line2D at 0x7ff39c6c9ad0>]
plot(angle(stft).T[25])
[<matplotlib.lines.Line2D at 0x7ff39c6062d0>]
Let's look at the unwrapped phase:
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')
[<matplotlib.lines.Line2D at 0x7ff39c52bf90>]
A line! The frequency is constant across all the bins!
delta_phs = diff(unwrapped)
plot(delta_phs)
[<matplotlib.lines.Line2D at 0x7ff39c463650>]
Now accumulate phase for the hop period:
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
(6.25, 39.269908169872416)
Adjusted phase:
phs_change = delta_phs[0]
phs_change
2.4488454978294771
Expected change:
total_phs_k % (2 * pi)
1.5707963267948983
Phase divergence:
div = phs_change - total_phs_k % (2 * pi)
div
0.8780491710345788
final_phs = total_phs_k + div
final_phs
40.147957340906991
Estimate number of cycles in hop period from total accumulated phase:
final_phs/(2*pi)
6.3897458658478943
Multiply by overlap:
(winsize/float(hopsize)) *final_phs/(2*pi)
25.558983463391577
Now times the number of windows to match the original full length:
(N/float(winsize))*(winsize/hopsize) *final_phs/(2*pi)
204.47186770713262
Yay! (or not...!)
Original frequency was 203.3, but this still got us much closer to the original frequency than what the FFT bins provide:
(N/float(winsize))*25
200.0
(N/float(winsize))*26
208.0
By: Andrés Cabrera mantaraya36@gmail.com For MAT course MAT 201A at UCSB
This ipython notebook is licensed under the CC-BY-NC-SA license: http://creativecommons.org/licenses/by-nc-sa/4.0/