def taper(flag, m, samples):
    # m = samples for taper
    # samples = total samples
    window = kaiser(2*m+1, beta=14)

    if flag == 'front':
        # cut and replace the second half of window with 1s
        ones = np.ones(samples-m-1)
        window = window[0:(m+1)]
        window = np.concatenate([window, ones])

    elif flag == 'end':
        # cut and replace the first half of window with 1s
        ones = np.ones(samples-m-1)
        window = window[(m+1):]
        window = np.concatenate([ones, window])

    elif flag == 'all':
        ones = np.ones(samples-2*m-1)
        window = np.concatenate([window[0:(m+1)], ones, window[(m+1):]])

    # avoid concatenate error
    if window.size < samples:
        window = np.append(window, 1)

    if window.size != samples:
        print(window.size)
        print(samples)
        print("[ERROR]: taper and data do not have the same number of samples.")
        window = np.ones(samples)

    return window

def apply_kaiserbessel_window(X, alpha=6.5):
    """
    Apply a Kaiser-Bessel window to X.

    Parameters
    ----------
    X : ndarray, shape=(n_samples, n_features)
        Input array of samples

    alpha : float, optional (default=6.5)
        Tuning parameter for Kaiser-Bessel function. alpha=6.5 should make
        perfect reconstruction possible for MDCT.

    Returns
    -------
    X_windowed : ndarray, shape=(n_samples, n_features)
        Windowed version of X.
    """
    beta = np.pi * alpha
    win = sg.kaiser(X.shape[1], beta)
    row_stride = 0
    col_stride = win.itemsize
    strided_win = as_strided(win, shape=X.shape,
                             strides=(row_stride, col_stride))
    return X * strided_win

def firwin_kaiser_lpf(f_pass, f_stop, d_stop, fs = 1.0, N_bump=0):
    """
    Design an FIR lowpass filter using the sinc() kernel and
    a Kaiser window. The filter order is determined based on 
    f_pass Hz, f_stop Hz, and the desired stopband attenuation
    d_stop in dB, all relative to a sampling rate of fs Hz.
    Note: the passband ripple cannot be set independent of the
    stopband attenuation.

    Mark Wickert October 2016
    """
    wc = 2*np.pi*(f_pass + f_stop)/2/fs
    delta_w = 2*np.pi*(f_stop - f_pass)/fs
    # Find the filter order
    M = np.ceil((d_stop - 8)/(2.285*delta_w))
    # Adjust filter order up or down as needed
    M += N_bump
    N_taps = M + 1
    # Obtain the Kaiser window
    beta = signal.kaiser_beta(d_stop)
    w_k = signal.kaiser(N_taps,beta)
    n = np.arange(N_taps)
    b_k = wc/np.pi*np.sinc(wc/np.pi*(n-M/2)) * w_k
    b_k /= np.sum(b_k)
    print('Kaiser Win filter taps = %d.' % N_taps)
    return b_k

 def test_spectrogram_energy_conservation(self):
     """Test the energy conservation property of the spectrogram."""
     signal, _ = fmlin(128, 0.1, 0.4)
     window = kaiser(17, 3 * np.pi)
     tfr, ts, freqs = cohen.Spectrogram(signal, n_fbins=64, fwindow=window).run()
     e_sig = (np.abs(signal) ** 2).sum()
     self.assertAlmostEqual(tfr.sum().sum() / 64, e_sig)

def THD(signal, sample_rate):
    """Measure the THD for a signal
    
    This function is not yet trustworthy.
    
    Returns the estimated fundamental frequency and the measured THD,
    calculated by finding peaks in the spectrum.
    
    There are two definitions for THD, a power ratio or an amplitude ratio
    When finished, this will list both
    
    """
    # Get rid of DC and window the signal
    signal -= mean(signal) # TODO: Do this in the frequency domain, and take any skirts with it?
    windowed = signal * kaiser(len(signal), 100)
    
    # Find the peak of the frequency spectrum (fundamental frequency)
    f = rfft(windowed)
    i = argmax(abs(f))
    true_i = parabolic(log(abs(f)), i)[0]
    print 'Frequency: %f Hz' % (sample_rate * (true_i / len(windowed)))
    
    print 'fundamental amplitude: %.3f' % abs(f[i])
    
    # Find the values for the first 15 harmonics.  Includes harmonic peaks only, by definition
    # TODO: Should peak-find near each one, not just assume that fundamental was perfectly estimated.
    # Instead of limited to 15, figure out how many fit based on f0 and sampling rate and report this "4 harmonics" and list the strength of each
    for x in range(2, 15):
        print '%.3f' % abs(f[i * x]),
   
    THD = sum([abs(f[i*x]) for x in range(2,15)]) / abs(f[i])
    print '\nTHD: %f%%' % (THD * 100)
    return

def freq_from_hps(signal, fs):
    """
    Estimate frequency using harmonic product spectrum

    Low frequency noise piles up and overwhelms the desired peaks

    Doesn't work well if signal doesn't have harmonics
    """
    signal = asarray(signal) + 0.0

    N = len(signal)
    signal -= mean(signal)  # Remove DC offset

    # Compute Fourier transform of windowed signal
    windowed = signal * kaiser(N, 100)

    # Get spectrum
    X = log(abs(rfft(windowed)))

    # Remove mean of spectrum (so sum is not increasingly offset
    # only in overlap region)
    X -= mean(X)

    # Downsample sum logs of spectra instead of multiplying
    hps = copy(X)
    for h in range(2, 9):  # TODO: choose a smarter upper limit
        dec = decimate(X, h, zero_phase=True)
        hps[:len(dec)] += dec

    # Find the peak and interpolate to get a more accurate peak
    i_peak = argmax(hps[:len(dec)])
    i_interp = parabolic(hps, i_peak)[0]

    # Convert to equivalent frequency
    return fs * i_interp / N  # Hz

def freq_from_fft(signal, fs):
    """
    Estimate frequency from peak of FFT

    Pros: Accurate, usually even more so than zero crossing counter
    (1000.000004 Hz for 1000 Hz, for instance).  Due to parabolic
    interpolation being a very good fit for windowed log FFT peaks?
    https://ccrma.stanford.edu/~jos/sasp/Quadratic_Interpolation_Spectral_Peaks.html
    Accuracy also increases with signal length

    Cons: Doesn't find the right value if harmonics are stronger than
    fundamental, which is common.
    """
    signal = asarray(signal)

    N = len(signal)

    # Compute Fourier transform of windowed signal
    windowed = signal * kaiser(N, 100)
    f = rfft(windowed)

    # Find the peak and interpolate to get a more accurate peak
    i_peak = argmax(abs(f))  # Just use this value for less-accurate result
    i_interp = parabolic(log(abs(f)), i_peak)[0]

    # Convert to equivalent frequency
    return fs * i_interp / N  # Hz

def firwin_kaiser_hpf(f_stop, f_pass, d_stop, fs = 1.0, N_bump=0):
    """
    Design an FIR highpass filter using the sinc() kernel and
    a Kaiser window. The filter order is determined based on 
    f_pass Hz, f_stop Hz, and the desired stopband attenuation
    d_stop in dB, all relative to a sampling rate of fs Hz.
    Note: the passband ripple cannot be set independent of the
    stopband attenuation.

    Mark Wickert October 2016
    """
    # Transform HPF critical frequencies to lowpass equivalent
    f_pass_eq = fs/2. - f_pass
    f_stop_eq = fs/2. - f_stop
    # Design LPF equivalent
    wc = 2*np.pi*(f_pass_eq + f_stop_eq)/2/fs
    delta_w = 2*np.pi*(f_stop_eq - f_pass_eq)/fs
    # Find the filter order
    M = np.ceil((d_stop - 8)/(2.285*delta_w))
    # Adjust filter order up or down as needed
    M += N_bump
    N_taps = M + 1
    # Obtain the Kaiser window
    beta = signal.kaiser_beta(d_stop)
    w_k = signal.kaiser(N_taps,beta)
    n = np.arange(N_taps)
    b_k = wc/np.pi*np.sinc(wc/np.pi*(n-M/2)) * w_k
    b_k /= np.sum(b_k)
    # Transform LPF equivalent to HPF
    n = np.arange(len(b_k))
    b_k *= (-1)**n
    print('Kaiser Win filter taps = %d.' % N_taps)
    return b_k

def freq_from_hps(signal, fs):
    """Estimate frequency using harmonic product spectrum
    
    Low frequency noise piles up and overwhelms the desired peaks
    """
    N = len(signal)
    signal -= mean(signal) # Remove DC offset
    
    # Compute Fourier transform of windowed signal
    windowed = signal * kaiser(N, 100)
    
    # Get spectrum
    X = log(abs(rfft(windowed)))
    
    # Downsample sum logs of spectra instead of multiplying
    hps = copy(X)
    for h in arange(2, 9): # TODO: choose a smarter upper limit
        dec = decimate(X, h)
        hps[:len(dec)] += dec
    
    # Find the peak and interpolate to get a more accurate peak
    i_peak = argmax(hps[:len(dec)])
    i_interp = parabolic(hps, i_peak)[0]
    
    # Convert to equivalent frequency
    return fs * i_interp / N # Hz

 def test_spectrogram_linearity(self):
     signal, _ = fmlin(128, 0.1, 0.4)
     window = kaiser(17, 3 * np.pi)
     tfr1, _, _ = cohen.spectrogram(signal, n_fbins=64, window=window)
     tfr2, _, _ = cohen.spectrogram(signal * 2, n_fbins=64, window=window)
     x = np.sum(np.sum(tfr2))
     y = np.sum(np.sum(tfr1))
     self.assertEqual(x / y, 4)

def calc_freq_fft(signal):
  N = len(signal)
  # Compute Fourier transform of windowed signal
  windowed = signal * ss.kaiser(N, 100)
  f = np.fft.rfft(windowed)
  # Find the peak and interpolate to get a more accurate peak
  i_peak = np.argmax(abs(f)) # Just use this value for less-accurate result
  i_interp = parabolic(np.log(abs(f)), i_peak)[0]
  # Convert to equivalent frequency
  return sample_rate * i_interp / N # Hz

 def test_spectrogram_linearity(self):
     """Test the linearity property of the spectrogram."""
     signal, _ = fmlin(128, 0.1, 0.4)
     window = kaiser(17, 3 * np.pi)
     tfr1, _, _ = cohen.Spectrogram(signal, n_fbins=64,
                                    fwindow=window).run()
     tfr2, _, _ = cohen.Spectrogram(signal * 2, n_fbins=64,
                                    fwindow=window).run()
     x = np.sum(np.sum(tfr2))
     y = np.sum(np.sum(tfr1))
     self.assertEqual(x / y, 4)

 def test_basic(self):
     assert_allclose(signal.kaiser(6, 0.5),
                     [0.9403061933191572, 0.9782962393705389,
                      0.9975765035372042, 0.9975765035372042,
                      0.9782962393705389, 0.9403061933191572])
     assert_allclose(signal.kaiser(7, 0.5),
                     [0.9403061933191572, 0.9732402256999829,
                      0.9932754654413773, 1.0, 0.9932754654413773,
                      0.9732402256999829, 0.9403061933191572])
     assert_allclose(signal.kaiser(6, 2.7),
                     [0.2603047507678832, 0.6648106293528054,
                      0.9582099802511439, 0.9582099802511439,
                      0.6648106293528054, 0.2603047507678832])
     assert_allclose(signal.kaiser(7, 2.7),
                     [0.2603047507678832, 0.5985765418119844,
                      0.8868495172060835, 1.0, 0.8868495172060835,
                      0.5985765418119844, 0.2603047507678832])
     assert_allclose(signal.kaiser(6, 2.7, False),
                     [0.2603047507678832, 0.5985765418119844,
                      0.8868495172060835, 1.0, 0.8868495172060835,
                      0.5985765418119844])

def CircularDeconvolution(x,y,dt,fmax,beta=None):
    
    from numpy.fft import rfft, irfft, fftshift, fft, ifft, fftfreq, ifftshift
    from numpy import floor,zeros,hstack,conj,pi,exp,linspace,ones,real,pi
    from scipy.signal import blackmanharris,kaiser,hann
    from matplotlib.pyplot import plot, show
    
    # mx = abs(x).argmax()
    # my = abs(y).argmax()
    # Nx = len(x)
   #  Ny = len(y)
    Nmin = len(x)+len(y)-1
    
   
    N = FFTLengthPower2((len(x)+len(y)-1))
    
    

    # X = rfft(hstack((x[mx::],zeros(abs(Ny-Nx)),x[0:mx-1])))
    X = fft(x,n=N)
    Sxx = fftshift(X*conj(X))
    # Y = rfft(y[m::],n=Ny)
    # Y = rfft(hstack((y[mx::],y[0:mx-1])))
    Y = fft(y,n=N)
    Syx = fftshift(Y*conj(X))
    
    f = fftshift(fftfreq(N,dt))
    
    fpass = [(f>=-fmax)&(f<=fmax)]
    
    H = Syx[fpass]/Sxx[fpass]
    
    
    
    if beta is None:
    
        H = hann(len(H))*H
        
    else:
        
        H = kaiser(len(H),beta)*H
        

    HH = zeros(N)+0.0*1j
    
    HH[fpass] = H.copy()
    
    H = ifftshift(HH)
  
    h = real(ifft(H))

    return h[0:Nmin], H

def kaiser_derived(M, beta):
    """ Return a Kaiser-Bessel derived window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero or less, an empty
        array is returned.
    beta : float
        Kaiser-Bessel window shape parameter.

    Returns
    -------
    w : ndarray
        The window, normalized to fulfil the Princen-Bradley condition.

    Notes
    -----
    This window is only defined for an even number of taps.

    References
    ----------
    .. [1] Wikipedia, "Kaiser window",
           https://en.wikipedia.org/wiki/Kaiser_window

    """
    M = int(M)
    try:
        from scipy.signal import kaiser_derived as scipy_kd
        return scipy_kd(M, beta)
    except ImportError:
        pass

    if M < 1:
        return np.array([])

    if M % 2:
        raise ValueError(
            "Kaiser Bessel Derived windows are only defined for even number "
            "of taps"
        )

    w = np.zeros(M)
    kaiserw = kaiser(M // 2 + 1, beta)
    csum = np.cumsum(kaiserw)
    halfw = np.sqrt(csum[:-1] / csum[-1])
    w[:M//2] = halfw
    w[-M//2:] = halfw[::-1]

    return w

def firwin_kaiser_bsf(f_stop1, f_pass1, f_pass2, f_stop2, d_stop, 
                  fs = 1.0, N_bump=0):
    """
    Design an FIR bandstop filter using the sinc() kernel and
    a Kaiser window. The filter order is determined based on 
    f_stop1 Hz, f_pass1 Hz, f_pass2 Hz, f_stop2 Hz, and the 
    desired stopband attenuation d_stop in dB for both stopbands,
    all relative to a sampling rate of fs Hz.
    Note: The passband ripple cannot be set independent of the
    stopband attenuation.
    Note: The filter order is forced to be even (odd number of taps)
    so there is a center tap that can be used to form 1 - H_BPF.

    Mark Wickert October 2016    
    """
    # First design a BPF starting from simple LPF equivalent
    # The upper and lower stopbands are assumed to have 
    # the same attenuation level. The LPF equivalent critical
    # frequencies:
    f_pass = (f_pass2 - f_pass1)/2
    f_stop = (f_stop2 - f_stop1)/2
    # Continue to design equivalent LPF
    wc = 2*np.pi*(f_pass + f_stop)/2/fs
    delta_w = 2*np.pi*(f_stop - f_pass)/fs
    # Find the filter order
    M = np.ceil((d_stop - 8)/(2.285*delta_w))
    # Adjust filter order up or down as needed
    M += N_bump
    # Make filter order even (odd number of taps)
    if ((M+1)/2.0-int((M+1)/2.0)) == 0:
        M += 1
    N_taps = M + 1
    # Obtain the Kaiser window
    beta = signal.kaiser_beta(d_stop)
    w_k = signal.kaiser(N_taps,beta)
    n = np.arange(N_taps)
    b_k = wc/np.pi*np.sinc(wc/np.pi*(n-M/2)) * w_k
    b_k /= np.sum(b_k)
    # Transform LPF to BPF
    f0 = (f_pass2 + f_pass1)/2
    w0 = 2*np.pi*f0/fs
    n = np.arange(len(b_k))
    b_k_bs = 2*b_k*np.cos(w0*(n-M/2))
    # Transform BPF to BSF via 1 - BPF for odd N_taps
    b_k_bs = -b_k_bs
    b_k_bs[int(M/2)] += 1 
    print('Kaiser Win filter taps = %d.' % N_taps)
    return b_k_bs

def analyze(filename):
  global correct_count

  print('analyzing ' + filename + ' ...')
  try:
    sampling_rate, signal = scipy.io.wavfile.read(filename)
  except:
    print("unable to analyze " + filename)
    print("")
  else:
    samples_count = len(signal)
    duration = float(samples_count) / sampling_rate
    if not isinstance(signal[0], np.int16):
      signal = [s[0] for s in signal]
    signal = signal * kaiser(samples_count, 100)

    spectrum = np.log(abs(np.fft.rfft(signal)))
    hps = copy(spectrum)
    for h in np.arange(2, 6):
      dec = decimate(spectrum, h)
      hps[:len(dec)] += dec
    peak_start = 50 * duration
    peak = np.argmax(hps[peak_start:])
    fundamental = (peak_start + peak) / duration

    if fundamental < 165:
      verdict = 'M'
    elif 180 < fundamental:
      verdict = 'F'
    else:
      verdict = 'U'

    correct = re.search('([KM])\.wav', filename).group(1)
    if correct == 'K':
      correct = 'F'

    if correct == verdict:
      correct_count += 1

    print("fundamendal frequency: " + "%.5f" % fundamental)
    print("verdict: " + verdict)
    if correct != verdict:
      print("WRONG")
    print("")

def wndow(win_size_N, out_vector_N = -1, **kwargs):
	if win_size_N > out_vector_N:
		out_vector_N = win_size_N
	
	wnd_size_div = out_vector_N / win_size_N
	
	wndw_kind = kwargs.get('kind', 'hamming')
	
	wndw = np.arange(1)
	
	if wndw_kind == 'kaiser':
		shaper = kwargs.get('shaper', 1)
		wndw_s = si.kaiser(win_size_N, shaper)
		wndw = np.concatenate([np.zeros((wnd_size_div -
		1)*out_vector_N/(2*wnd_size_div)), wndw_s])
		wndw = np.concatenate([wndw, np.zeros((wnd_size_div -
		1)*out_vector_N/(2*wnd_size_div))])

	return wndw

 def test_spectrogram_time_invariance(self):
     """Test the time invariance property of the spectrogram."""
     signal, _ = fmlin(128, 0.1, 0.4)
     window = kaiser(17, 3 * np.pi)
     tfr, ts, freqs = cohen.Spectrogram(signal, n_fbins=64, fwindow=window).run()
     shift = 64
     timeshifted_signal = np.roll(signal, shift)
     timeshifted_tfr, _, _ = cohen.Spectrogram(timeshifted_signal, n_fbins=64,
                                         fwindow=window).run()
     rolled_tfr = np.roll(tfr, shift, axis=1)
     # the time invariance property holds mostly qualitatively. The shifted
     # TFR is not numerically indentical to the rolled TFR, having
     # differences at the edges; so clip with two TFRs where there are
     # discontinuities in the TFR.
     edge = 10
     xx = np.c_[timeshifted_tfr[:, edge:(shift - edge)],
                timeshifted_tfr[:, (shift + edge):-edge]]
     yy = np.c_[rolled_tfr[:, edge:(shift - edge)],
                rolled_tfr[:, (shift + edge):-edge]]
     np.testing.assert_allclose(xx, yy)

def genderRecognition(filename):
  global countrights
  try:
    rate, signal = sc.read(filename)
  except:
    print("")
  else:
    frames = len(signal)
    dur = float(frames) / rate
    if not isinstance(signal[0], np.int16):
      signal = [s[0] for s in signal]
    try:
      signal = signal * kaiser(frames, 30)
    except:
      print("")
    else:
      
      signal1 = abs(np.fft.rfft(signal))#np.log(abs(np.fft.rfft(signal)))
      signal2 = copy(signal1)
      
      for h in np.arange(2, 6):
        decsignal = decimate(signal1, h)
        signal2[:len(decsignal)] += decsignal
      
      speak = 50 * dur
      peak = np.argmax(signal2[speak:])
      freq = (speak + peak) / dur

      if freq < 180:
        result = 'M'
      elif 180 <= freq:
        result = 'K'

      truth = re.search('([KM])\.wav', filename).group(1)

      if truth == result:
        countrights += 1