from __future__ import division
import pyaudio
import wave

import sys
import scipy
import numpy as np
import struct


#from scikits.audiolab import flacread
from numpy.fft import rfft, irfft
from numpy import argmax, sqrt, mean, diff, log
import matplotlib
from scipy.signal import blackmanharris, fftconvolve
from time import time
from parabolic import parabolic

#for frequency estimation using harmonic product spectrum
from pylab import subplot, plot, log, copy, show

chunk = 2**8

def freq_from_fft(sig, fs):
	"""Estimate frequency from peak of FFT
	"""
	# Compute Fourier transform of windowed signal
	windowed = sig * blackmanharris(chunk)
	fftData = rfft(windowed)
	# Find the peak and interpolate to get a more accurate peak
	i = argmax(abs(fftData)) # Just use this for less-accurate, naive version
	if i != len(fftData)-1:
		true_i = parabolic(log(abs(fftData)), i)[0]
		# Convert to equivalent frequency
		return fs * true_i / len(windowed)
	else:
		return fs * i / len(windowed) 

def freq_from_fft_quadratic(sig, rate):
	fftData = np.abs(np.fft.rfft(sig))**2
	#print(fftData[0:4])
	# find the maximum
	fftMax = fftData[1:].argmin() + 1