compute frequencies

fft.py

import matplotlib.pyplot as plt
import numpy as np
import scipy
import scipy.fftpack
import scipy.io.wavfile
import heapq
from datetime import datetime
# from itertools import izip

Fs, y = scipy.io.wavfile.read("sine.wav")
# Fs, y = scipy.io.wavfile.read("noise.wav")
# Fs, y = scipy.io.wavfile.read("pluck.wav")
# Fs, y_orig = scipy.io.wavfile.read("sine_8h.wav")
# Fs, y = scipy.io.wavfile.read("output/0.wav")

Ts = 1.0/Fs  # sampling interval
n = len(y)  # length of the signal
t = np.arange(0, len(y)/Fs, Ts)  # time vector

k = np.arange(n)
T = n/Fs
frq = k/T  # two sides frequency range
frq = frq[range(n/2)]  # one side frequency range

Y = np.fft.fft(y)/n  # fft computing and normalization
Y = Y[range(n/2)]
yy = abs(Y)
yy = yy/yy.max()

result = set()
starttime = datetime.now()

ind = np.argpartition(yy, -15)[-15:]
for i in ind:
    if yy[i] > 0.5:
#        print(frq[i], yy[i])
        result.add(frq[i])

# for f, v in izip(frq[ind], yy[ind]):
#     if v > 0.5:
#         print(f, v)
#         result.add(f)

# largests = heapq.nlargest(15, xrange(len(yy)), yy.take)
# for l in largests:
#     if yy[l] > 0.5:
#         print(frq[l])
#         result.add(frq[l])

# for i in range(len(frq)):
#     if yy[i] > 5000:
#         print(frq[i], yy[i])
#         result.add(frq[i])

print("result:", list(result))
endtime = datetime.now()
# print("took {0} time".format(endtime - starttime))

# fig, ax = plt.subplots(2, 1)
# ax[0].plot(t, y)
# ax[0].set_xlabel('Time2')
# ax[0].set_ylabel('Amplitude')
# ax[1].plot(frq, abs(Y), 'r')  # plotting the spectrum
# ax[1].set_xlabel('Freq (Hz)')
# ax[1].set_ylabel('|Y(freq)|')

plt.plot(frq, yy, 'r')  # plotting the spectrum
plt.xlabel('Freq (Hz)')
plt.ylabel('|Y(freq)|')

plt.show()