def GenerateBLI(sampling_rate = 48000.0,
                cutoff = 15000.0,
                length = 4,
                ppiv = 2700,
                beta = 8.3,
                apodization_factor = 0.5,
                apodization_beta = 0.5):
    """
    Generates a bandlimited impulse
    """
    kEpsilon = 1e-7
    points_count = ppiv * length
    xaxis = numpy.linspace(0.0, points_count, points_count)
    x2 = length * 2.0 * (xaxis - (points_count) / 2.0 + kEpsilon) / (points_count)
    x3 = numpy.pi * cutoff / sampling_rate * x2

    brickwall_impulse = numpy.sin(x3) / x3

    kaiser_window = numpy.kaiser(points_count, beta)
    bli_data = numpy.transpose(brickwall_impulse) * kaiser_window

    # Apodization
    apodization_window = (1.0 - apodization_factor) \
                            * numpy.kaiser(points_count,apodization_beta)
    bli_data *= apodization_window

    return bli_data

def plot_fourier(t, y, subplotnums=[[2,1,1],[2,1,2]], savedata="none", kaiserwindowparameter=0):
   ''' Function for plotting fourier series of some given arrays t and y
       :param t         time variable (plotted along the x-axis)
       :param y         variable to be plotted along the y-axis
       :subplotnums     subplotnum for plotting multiple plots in one window
   '''
   if len(subplotnums) - 1 != len(np.atleast_1d(kaiserwindowparameter)):
      print "Bad subplotnum vs kaiserwindowparameter"
      return
   # First check the t array whether it has a constant dt
   dt = t[1] - t[0]
   for i in xrange(len(t)-1):
      if dt != t[i+1] - t[i]:
         print "Gave bad timestep to plot_fourier, the time step in array t must be constant (for now)"
   # Use kaiser window on y
   _y = y*np.kaiser(len(y), np.atleast_1d(kaiserwindowparameter)[0])
   # Plot the raw data as well as the fitted data (fourier series)
   pl.subplot(subplotnums[0][0],subplotnums[0][1],subplotnums[0][2])
   pl.plot(t,_y,'.',color='r')
   # Do FFT on the data
   fourier=np.fft.fft(_y)*(1/(float)(len(t)))
   # Get frequencies of the fourier
   freq=np.fft.fftfreq(len(fourier), d=dt)
   # Declare t2 (Note: This is the same as t but we want the steps to be thicker so the data looks smoother
   dt2=dt*0.01
   t2=np.arange(len(t)*100)*dt2
   # Declare y2
   y2=np.array([np.sum(fourier*np.exp(complex(0,1)*2*np.pi*freq*T)) for T in t2])
   pl.plot(t2,y2,'-',color='b')
   pl.legend(["data", "fourier_fit"])
   # Plot the frequency spectrums
   j = 0
   for i in subplotnums[1:]:
      pl.subplot(i[0],i[1],i[2])
      # Use window function on the data
      _y = y*np.kaiser(len(y), np.atleast_1d(kaiserwindowparameter)[j])
      # Do FFT on the data
      fourier=np.fft.fft(_y)*(1/(float)(len(t)))
      # Get frequencies of the fourier
      freq=np.fft.fftfreq(len(fourier), d=dt)
      # Get the indexes:
      toIndex = (int)((len(freq)/2)/2.0 + 1)
      #maxValue = np.max(np.abs(fourier[1:len(fourier)/2]))
      #while True:
      #   if toIndex <= 2 or np.abs(fourier[toIndex-1]) > 0.02*maxValue:
      #      break
      #   toIndex = toIndex - 1
      #pl.plot(freq[1:len(freq)/2],2*np.abs(fourier[1:len(fourier)/2]), marker='.', linestyle='-', linewidth=0.5)
      if savedata != "none":
         saveplotdata( freq[1:toIndex], np.log(2*np.abs(fourier[1:toIndex])), savedata + "_" + str(j) + ".npy" )
      pl.plot(freq[1:toIndex],2*np.abs(fourier[1:toIndex]), marker='.', linestyle='-', linewidth=0.5)
      pl.legend(["frequency_spectrum"])
      pl.ylim([0,1.05*max(2*np.abs(fourier[1:len(fourier)/2]))])
      #pl.xlim([0, 0.13])
      #xTicks = np.arange(14)/100.0
      #pl.xticks(xTicks)
      # Put interactive mode on and show the plot:
      #pl.ion()
      #pl.show()
      j = j + 1

 def make_sense(self,u0):
     st=self.st
     L=numpy.shape(u0)[-1]
     u0dims= numpy.ndim(u0)
     
     rows=numpy.shape(u0)[0]
     cols=numpy.shape(u0)[1]
     # dpss rely on ctypes which are not suitable for cx_freeze 
     #    [dpss_rows, eigens] = spectrum.mtm.dpss(rows, 50, 1)
     #    [dpss_cols, eigens] = spectrum.mtm.dpss(cols, 50, 1)
     dpss_rows = numpy.kaiser(rows, 100)
     dpss_cols = numpy.kaiser(cols, 100)
     
     dpss_rows = numpy.fft.fftshift(dpss_rows)
     dpss_cols = numpy.fft.fftshift(dpss_cols)
     
     dpss_rows = appendmat(dpss_rows,cols)
     dpss_cols  = appendmat(dpss_cols,rows)
     
     
     dpss_fil=dpss_rows*dpss_cols.T # low pass filter
     
     rms=numpy.sqrt(numpy.mean(u0*u0.conj(),-1)) # Root of sum square
     st['sensemap']=numpy.ones(numpy.shape(u0),dtype=numpy.complex64)
     
     #    print('L',L)
     #    print('rms',numpy.shape(rms))
     for ll in numpy.arange(0,L):
         st['sensemap'][:,:,ll]=(u0[:,:,ll]+1e-15)/(rms+1e-15)
         st['sensemap'][:,:,ll]=scipy.fftpack.ifft2(scipy.fftpack.fft2(st['sensemap'][:,:,ll])*dpss_fil)
     return st

    def _kaiser_window(self, width, beta):
        """
        Generates a kaiser window

        beta Window shape
        0    Rectangular
        5    Similar to a Hamming
        6    Similar to a Hann
        8.6  Similar to a Blackman
        """
        return np.kaiser(width[0], beta).reshape((-1,1,1)) * np.kaiser(width[1], beta).reshape((-1,1)) * np.kaiser(width[2], beta)

 def signal_filter(self, channel, freqs=(1, 15)):
     """Filter signal from specific channel. A FFT is performed in
     the signal, and the result is multiplied by a window function
     (kaiser function), which nulls the undesired beating
     frequencies, that are outside of 'freqs'. The signal is then
     recreated with an IFFT."""
     # zero padding size:
     zer_pad_filter = 4
     # alias for sweep_size
     N = int(self.sweep_size)
     # FFT with zero padding
     fft = np.fft.rfft(self.single_sweep_data[channel], zer_pad_filter * N)
     # bp=fft[:] #used for other plot functions
     fmin, fmax = freqs
     # creates the beating frequency axis, and finds the position of
     # the frequency limits in the axis
     fft_freq = np.linspace(0, self.rate  / 2.,
                            num=zer_pad_filter * N / 2)
     cmin = (abs(fft_freq - fmin)).argmin()
     cmax = (abs(fft_freq - fmax)).argmin()
     # creates window function for filter. Second argument of kaiser
     # function must be float
     window_func = np.concatenate((np.zeros(cmin + 1),
                                   np.kaiser(cmax - cmin, 2.), np.zeros(zer_pad_filter * N / 2 - cmax)))
     # multiply the window by the signal's FFT.
     bp = np.multiply(fft, window_func)
     # Since the signal is REAL, we use IRFFT, and takes only the
     # first half, since the other half is symmetric.
     newsig = np.fft.irfft(bp)[:N]
     return newsig

def main():
    import sys
    from scipy import misc
    softness = 2
    r_scale = 3
    windowmaker = lambda x: np.kaiser(x, softness)
    for infilename in sys.argv[1:]:
        print("File: %s" % infilename)
        array = misc.imread(infilename)
        circles = get_circles(array, 5)
        l_pnoise, l_order, c_order, r_order, r_pnoise = circles

        window = get_holed_window(windowmaker, l_order[2] * r_scale, 10)
        mask = get_mask(array.shape, window, l_order[1])
        window = windowmaker(l_pnoise[2] * r_scale)
        mask *= 1 - get_mask(array.shape, window, l_pnoise[1])
        window = windowmaker(c_order[2] * r_scale)
        mask *= 1 - get_mask(array.shape, window, c_order[1])
#        showimage(mask * 255)
#        showimage(logscale(mask * array))

        window = get_holed_window(windowmaker, r_order[2] * r_scale, 10)
        mask = get_mask(array.shape, window, r_order[1])
        window = windowmaker(r_pnoise[2] * r_scale)
        mask *= 1 - get_mask(array.shape, window, r_pnoise[1])
        window = windowmaker(c_order[2] * r_scale)
        mask *= 1 - get_mask(array.shape, window, c_order[1])

def periodogram(x,y,z,t,name):
    vel=x*x+y*y
    components=[vel,z*z]
    F=[]
    A=[]

    for i in range(0,len(components)):
        window=np.kaiser(components[i].shape[-1],5)
        vel=components[i]*window
        f,a=DO_FFT(vel,20)
        F.append(f)
        A.append(a)

    
    plt.figure(1,figsize=(11,7))
    plt.title("Spectral Gap - LiCor Fontanella1 Precampagna")
    plt.subplot(121)
    plt.title("FFT horizontal velocity")
    plt.xlabel("Frequency [Hz]")
    plt.ylabel("Power Spectrum [dB]")
    plt.loglog(F[0],A[0],'k', markersize=0.1)

    plt.subplot(122)
    plt.title("FFT vertical velocity")
    plt.xlabel("Frequency [Hz]")
    plt.ylabel("Power Spectrum [dB]")
    plt.loglog(F[1],A[1],'r',markersize=0.1)

    try:
        plt.savefig("graph/"+name+"_FFT.png", figuresize=(8,6), dpi=320, format="png")
        print_ok("Graph saved in: "+"graph/"+name+"_FFT.png")
    except IOError as IoE:
        print_fail("I/O Error! Erro number = {0}; {1}".format(IoE.errno,IoE.strerror))
        exit(IoE.errno)

def kaiser(active):
    """Kaiser tapering window."""
    idx = _windowidx(active)
    # compute coefficients
    window = np.zeros(active.shape)
    window[idx] = np.kaiser(len(idx), 2)
    return window

def hps(x, fs=44100, lf=255, harmonics=3, precision=2, window=lambda l:_np.kaiser(l, 7.14285)):
    """ Estimates the pitch (fundamental frequency) of the given sample array by a standard HPS implementation. """
    x -= _np.mean(x)
    N = x.size
    w = x*window(N)

    # Append zeros to the end of the window so that each bin has at least the desired precision.
    if fs/N > precision:
        delta = int(fs/precision) - N
        w = _np.append(w, _np.zeros(delta))
        N = w.size

    X = _np.log(_np.abs(_np.fft.rfft(w)))

    # Sequentially decimate 'X' 'harmonics' times and add it to itself.
    # 'precision < fs/N' must hold, lest the decimation loses all the precision we'd gain.
    hps = _np.copy(X)
    for h in range(2, 2 + harmonics):
        dec = _sig.decimate(X, h)
        hps[:dec.size] += dec*(0.8**h)

    # Find the bin corresponding to the lowest detectable frequency.
    lb = lf*N/fs

    # And then the bin with the highest spectral content.
    arg_peak = lb + _np.argmax(hps[lb:dec.size])

    # TODO: Return the full array? A ranked list of identified notes?
    return fs*arg_peak/N

def single_taper_spectrum(data, delta, taper_name=None):
    """
    Returns the spectrum and the corresponding frequencies for data with the
    given taper.
    """
    length = len(data)
    good_length = length // 2 + 1
    # Create the frequencies.
    # XXX: This might be some kind of hack
    freq = abs(np.fft.fftfreq(length, delta)[:good_length])
    # Create the tapers.
    if taper_name == 'bartlett':
        taper = np.bartlett(length)
    elif taper_name == 'blackman':
        taper = np.blackman(length)
    elif taper_name == 'boxcar':
        taper = np.ones(length)
    elif taper_name == 'hamming':
        taper = np.hamming(length)
    elif taper_name == 'hanning':
        taper = np.hanning(length)
    elif 'kaiser' in taper_name:
        taper = np.kaiser(length, beta=14)
    # Should never happen.
    else:
        msg = 'Something went wrong.'
        raise Exception(msg)
    # Detrend the data.
    data = detrend(data)
    # Apply the taper.
    data *= taper
    spec = abs(np.fft.rfft(data)) ** 2
    return spec, freq

def preprocess(X, y=None, box_width=256, overlap=32, pad_width=0):
    X = pre_flow(X)
    
    out_X = []
    out_y = []
    m = int(X.shape[0] / box_width) * box_width
    
    w = np.kaiser(box_width + pad_width, 14)

    for start in range(0, m, overlap):
        end = start + box_width
        if end >= len(X):
            break
        
        x = X[start:end,]

        if pad_width > 0:
            x = np.vstack([x, np.zeros((pad_width, x.shape[1]))])

        x = (x.T * w).T

        psd = abs(np.fft.rfft(x, axis=0))

        #freqs, psd = signal.welch(X[start:end,], axis=0)
    
        out_X.append(psd.flatten())
        if y != None:
            out_y.append(stats.mode(y[start:end])[0][0])

    if y != None:
        return np.array(out_X), np.array(out_y)
    else:
        return np.array(out_X)

def simple_segmentation( wav_data, sample_rate, db_diff = 17, max_silence = 0.8 ): 

    # filter data
    filtered_data = filter_chain(wav_data, sample_rate, 1000, 0.6)

    # create figure and axis
    f, ax = plt.subplots()

    # get spectograma nd spectogram plot (i.e. STFT)
    NFFT = 1024
    noverlap = 256
    window = np.kaiser(NFFT, 8)
    spectogram, freqs, time_array, im = ax.specgram(filtered_data, NFFT=NFFT,
                                                    Fs=sample_rate, noverlap=noverlap,
                                                    window=window, cmap = 'Greys')
    ax.set_xlim(0, len(wav_data)/float(sample_rate)) # remove empty plot region
    ax.set_ylim(0, sample_rate/2.)                   # remove empty plot region

    syllables = find_syllables(spectogram, time_array, db_diff)

    segments = join_in_segments(syllables, max_silence) 

    # add segments to plot as red shadowed areas
    hspans = []
    for segment in segments:
        hspans.append(ax.axvspan(segment[0], segment[1], facecolor='r', alpha = 0.2))
    
    return segments, f, ax

def plot_Temperature(fig_size=None):
	#{{{
	N = 256
	hammingWindow = np.hamming(N)
	hanningWindow = np.hanning(N)
	kaiserWindow = np.kaiser(N,5)
	fig = plt.figure()
	if fig_size != None:
		fig.set_size_inches(fig_size[0], fig_size[1])
	for i in range(0,10):
		Time = Array[i-1][:,0]
		PlData1 = Array[i-1][:, 5]
		PlData2 = Array[i-1][:, 6]
		PlData3 = Array[i-1][:, 7]
		PlData4 = Array[i-1][:, 8]
		PlData5 = Array[i-1][:, 9]
		PlData6 = Array[i-1][:,10]
		PlData7 = Array[i-1][:,16]
		ax = fig.add_subplot(5,2,i)
		#plt.legend((NameIndex[5],NameIndex[6],NameIndex[7],NameIndex[8],NameIndex[9],NameIndex[10]),'upper left')
		#plt.plot(Time,PlData1,label=' 5-OriTemp')
		#plt.plot(Time,PlData2,label=' 6-InjTemp')
		#plt.plot(Time,PlData3,label=' 7-NUWTemp')
		plt.plot(Time,PlData4,label=' 8-NUCTemp')
		plt.plot(Time,PlData5,label=' 9-NDCTemp')
		plt.plot(Time,PlData6,label='10-FrnTemp')
		plt.plot(Time,PlData7,label='16-MixBoil')
		plt.legend(fontsize=15,loc = 'upper right')
		plt.xlabel('Time[s]')
		plt.ylabel('Temperature[K]')
		plt.xlim(0,10)
		plt.ylim(-200,1400)

def kaiser(x, width=None, weights=None, do_fit_edges=False):
    """Smooth the values in `x` with the Kaiser windowed filter.

    See: https://en.wikipedia.org/wiki/Kaiser_window

    Parameters
    ----------
    x : array-like
        1-dimensional numeric data set.
    width : float
        Fraction of x's total length to include in the rolling window (i.e. the
        proportional window width), or the integer size of the window.
    """
    if len(x) < 2:
        return x
    if width is None:
        width = guess_window_size(x, weights)
    x, wing, signal = check_inputs(x, width, False)
    # Apply signal smoothing
    window = np.kaiser(2 * wing + 1, 14)
    if weights is not None:
        y, _w = convolve_weighted(window, signal, weights)
    else:
        y = convolve_unweighted(window, signal, wing)
    if do_fit_edges:
        _fit_edges(x, y, wing)  # In-place
    return y

def gen_map(vals, n=nbins, hue=False):
	# saturation
	fvals = vals.flatten()
	yh, xh, patches = plt.hist(fvals, bins=n, range=(0,1), normed=False , cumulative=False , histtype='step')
	if hue:
		# apply window
		M = 9
		win = np.kaiser(M, 3.0)
		yh = np.insert(yh, 0, np.zeros(M/2))
		yh = np.append(yh, np.zeros(M/2))
		yh = rolling_window(yh.T, M)
		yh = np.dot(yh, win)
	yh /= sum(yh)
	if hue:
		# adapted norm
		#yh = np.minimum(yh, hcut)
		yh[yh<=hcut] = 0
		yh /= sum(yh)
	yh = np.cumsum(yh)

	xhi = np.linspace(0,1,256)
	yhi = np.interp(xhi, yh, xh[1:])
	yhinv = np.interp(xhi, xh[1:], yh)
	#plt.plot(xhi, yhi)
	return (yhi, yhinv)

def plotWindowFunc():
	plt.clf()
	plt.plot(np.arange(20), np.kaiser(20,3.5))
	plt.plot(np.arange(20), np.bartlett(20))
	plt.plot(np.arange(20), np.blackman(20))
	plt.plot(np.arange(20), np.hamming(20))
	plt.plot(np.arange(20), np.hanning(20))

def smooth(x,beta):
    """ kaiser window smoothing """
    window_len=11
    s = numpy.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]]
    w = numpy.kaiser(window_len,beta)
    y = numpy.convolve(w/w.sum(),s,mode='valid')
    return y[5:len(y)-5]

    def lowpass_filter(self,x,width):
        #wndw = np.sinc(np.r_[-15:16]/np.pi)/np.pi
        wndw = np.kaiser(width,6)
        wndw /= np.sum(wndw)
        new_array = signal.fftconvolve(x, wndw)

        return new_array[int(width/2):x.size+int(width/2)]

def test_snr(param, freq):
    """Test the SNR using a test tone at the given frequency.

    This works by computing the FFT of the result, zeroing the FFT
    corresponding to the test tone, and comparing the signal power
    with the signal power of the output signal.  Returns the SNR ratio
    (as a ratio, not dB).
    """
    w = 2 * math.pi * freq / param.RATE_IN
    length = param.TEST_LENGTH
    beta = param.TEST_BETA
    indata = (
        0.5 * numpy.sin(w * numpy.arange(length, dtype='float64')))
    outdata = resample_arr(param, indata)
    window = numpy.kaiser(len(outdata), beta)
    outdata *= window
    fft = numpy.fft.rfft(outdata)

    nbins = NBINS
    fbin = round(freq * len(outdata) / param.RATE_OUT)
    bin0 = min(max(fbin-nbins/2, 0), len(fft))
    bin1 = min(max(fbin-nbins/2 + nbins, 0), len(fft))
    fft[bin0:bin1] = 0
    noise = numpy.std(fft) / math.sqrt(len(outdata))
    signal = numpy.average(window)

    return signal / noise

def resample(s, p, q, h=None):
    """Change sampling rate by rational factor. This implementation is based on
    the Octave implementation of the resample function. It designs the 
    anti-aliasing filter using the window approach applying a Kaiser window with
    the beta term calculated as specified by [2].
    
    Ref [1] J. G. Proakis and D. G. Manolakis,
    Digital Signal Processing: Principles, Algorithms, and Applications,
    4th ed., Prentice Hall, 2007. Chap. 6
    Ref [2] A. V. Oppenheim, R. W. Schafer and J. R. Buck, 
    Discrete-time signal processing, Signal processing series,
    Prentice-Hall, 1999
    """
    gcd = fractions.gcd(p,q)
    if gcd>1:
        p=p/gcd
        q=q/gcd
    
    if h is None: #design filter
        #properties of the antialiasing filter
        log10_rejection = -3.0
        stopband_cutoff_f = 1.0/(2.0 * max(p,q))
        roll_off_width = stopband_cutoff_f / 10.0
    
        #determine filter length
        #use empirical formula from [2] Chap 7, Eq. (7.63) p 476
        rejection_db = -20.0*log10_rejection;
        l = numpy.ceil((rejection_db-8.0) / (28.714 * roll_off_width))
  
        #ideal sinc filter
        t = numpy.arange(-l, l + 1)
        ideal_filter=2*p*stopband_cutoff_f*numpy.sinc(2*stopband_cutoff_f*t)  
  
        #determine parameter of Kaiser window
        #use empirical formula from [2] Chap 7, Eq. (7.62) p 474
        beta = signal.kaiser_beta(rejection_db)
          
        #apodize ideal filter response
        h = numpy.kaiser(2*l+1, beta)*ideal_filter

    ls = len(s)
    lh = len(h)

    l = (lh - 1)/2.0
    ly = numpy.ceil(ls*p/float(q))

    #pre and postpad filter response
    nz_pre = numpy.floor(q - numpy.mod(l,q))
    hpad = h[-lh+nz_pre:]

    offset = numpy.floor((l+nz_pre)/q)
    nz_post = 0;
    while numpy.ceil(((ls-1)*p + nz_pre + lh + nz_post )/q ) - offset < ly:
        nz_post += 1
    hpad = hpad[:lh + nz_pre + nz_post]

    #filtering
    xfilt = upfirdn(s, hpad, p, q)

    return xfilt[offset-1:offset-1+ly]