def zpFFTsizeExpt(x, fs):
    """
    Inputs:
        x (numpy array) = input signal (2*M = 512 samples long)
        fs (float) = sampling frequency in Hz
    Output:
        The function should return a tuple (mX1_80, mX2_80, mX3_80)
        mX1_80 (numpy array): The first 80 samples of the magnitude spectrum output of dftAnal for Case-1
        mX2_80 (numpy array): The first 80 samples of the magnitude spectrum output of dftAnal for Case-2
        mX3_80 (numpy array): The first 80 samples of the magnitude spectrum output of dftAnal for Case-3
        
    The first few lines of the code to generate xseg and the windows have been written for you, 
    please use it and do not modify it. 
    """
    
    M = len(x)/2
    xseg = x[:M]
    w1 = get_window('hamming',M)
    w2 = get_window('hamming',2*M)

    mX1, P = dftAnal(xseg, w1, 256)
    mX2, P = dftAnal(x, w2, 512)
    mX3, P = dftAnal(xseg, w1, 512)

    return (mX1[:80], mX2[:80], mX3[:80])

def B_ell(flat_map, component, bin_size, window=None):
    """
    Return the binned beam function of the given flat map component,
    using ell bins of bin_size. If a window name is given, multiply
    the map component by the window function before taking the 2-D
    DFT.

    The returned beam function is the absolute value of the DFT, so it
    has the same units as the map component. The returned bins array
    contains the left edges of the bins corresponding to the returned
    data: for all i in range(bins.size),
    bins[i] <= data[i] < bins[i+1]
    """
    ell_x, ell_y= np.meshgrid(fftshift(fftfreq(flat_map.x.size, flat_map.dx()/(2*np.pi))),
                              fftshift(fftfreq(flat_map.y.size, flat_map.dy()/(2*np.pi))))
    ell_x = ell_x.T
    ell_y = ell_y.T
    ell_r = np.sqrt(ell_x**2 + ell_y**2)
    ell_bins = np.arange(0, np.max(ell_r), bin_size)
    beam = flat_map[component]
    if window is not None:
        beam *= get_window(window, flat_map.y.size)
        beam *= get_window(window, flat_map.x.size)[:, np.newaxis]
    dft = fftshift(fft2(beam))
    # Shift the indices down by 1 so that they correspond to the data
    # indices. These indices should always be valid indices for the
    # DFT data because r_ell has a lower bound at 0 and max(r_ell)
    # will have index ell_bins.size.
    bin_indices = np.digitize(ell_r.flatten(), ell_bins) - 1
    binned = np.zeros(ell_bins.size) # dtype?
    for i in range(binned.size):
        binned[i] = np.sqrt(np.mean(abs(dft.flatten()[i==bin_indices])**2))
    return ell_bins, binned, dft

def zpFFTsizeExpt(x, fs):
    """
    Inputs:
        x (numpy array) = input signal (2*M samples long)
        fs (float) = sampling frequency in Hz
    Output:
        The function should return a tuple (mX1_80, mX2_80, mX3_80)
        mX1_80 (numpy array): The first 80 samples of the magnitude spectrum output of dftAnal for Case-1
        mX2_80 (numpy array): The first 80 samples of the magnitude spectrum output of dftAnal for Case-2
        mX3_80 (numpy array): The first 80 samples of the magnitude spectrum output of dftAnal for Case-3
        
    The first few lines of the code to generate xseg and the windows have been written for you, 
    please use it and do not modify it. 
    """
    M = len(x)/2
    xseg = x[:M]
    w1 = get_window('hamming',M)
    w2 = get_window('hamming',2*M)
    ## Your code here 
    (xM1, pX1) = dftAnal(xseg, w1, M)
    (xM2, pX2) = dftAnal(x, w2, len(x))
    (xM3, pX3) = dftAnal(xseg, w1, len(x))

    plt.plot(np.arange(0,80,2), xM1[:40], color="red")
    plt.plot(xM2[:80], color="blue")
    plt.plot(xM3[:80], color="green")
    plt.show()

    return (xM1[:80], xM2[:80], xM3[:80])

    def test_boxcar(self):
        w = signal.get_window('boxcar', 12)
        assert_array_equal(w, np.ones_like(w))

        # window is a tuple of len 1
        w = signal.get_window(('boxcar',), 16)
        assert_array_equal(w, np.ones_like(w))

def zpFFTsizeExpt(x, fs):
    """
    Inputs:
        x (numpy array) = input signal (2*M samples long)
        fs (float) = sampling frequency in Hz
    Output:
        The function should return a tuple (mX1_80, mX2_80, mX3_80)
        mX1_80 (numpy array): The first 80 samples of the magnitude spectrum output of dftAnal for Case-1
        mX2_80 (numpy array): The first 80 samples of the magnitude spectrum output of dftAnal for Case-2
        mX3_80 (numpy array): The first 80 samples of the magnitude spectrum output of dftAnal for Case-3
        
    The first few lines of the code to generate xseg and the windows have been written for you, 
    please use it and do not modify it. 
    """
    M = len(x)/2
    xseg = x[:M]
    w1 = get_window('hamming',M)
    w2 = get_window('hamming',2*M)
    ## Your code here 
    # NOTE: The shape of the window is affected by the size parameter M. So if you want a 256 size window
    #you need to specifically create it and cant use w[:256]
    # Case-1: Input signal xseg (256 samples), window w1 (256 samples), and FFT size of 256
    N = 256
    mX1_80, pX1 = dftAnal(xseg, w1, N)
    #Input signal x (512 samples), window w2 (512 samples), and FFT size of 512
    N = 512
    mX2_80, pX2 = dftAnal(x[:N], w2[:N], N)
    #Input signal xseg (256 samples), window w1 (256 samples), and FFT size of 512 (Implicitly does a zero-padding of xseg by 256 samples)
    N = 256
    mX3_80, pX3 = dftAnal(x[:N], w1, 512)
    
    return(mX1_80[:80],mX2_80[:80],mX3_80[:80])

def zpFFTsizeExpt(x, fs):
    """
    Inputs:
        x (numpy array) = input signal (2*M samples long)
        fs (float) = sampling frequency in Hz
    Output:
        The function should return a tuple (mX1_80, mX2_80, mX3_80)
        mX1_80 (numpy array): The first 80 samples of the magnitude spectrum output of dftAnal for Case-1
        mX2_80 (numpy array): The first 80 samples of the magnitude spectrum output of dftAnal for Case-2
        mX3_80 (numpy array): The first 80 samples of the magnitude spectrum output of dftAnal for Case-3
        
    The first few lines of the code to generate xseg and the windows have been written for you, 
    please use it and do not modify it. 
    """
    ## Your code here
    M = len(x)/2
    xseg = x[:M]
    w1 = get_window('hamming',M)
    w2 = get_window('hamming',2*M)
    
    mX1 = dftAnal(xseg, w1, M)[0]
    mX2 = dftAnal(x, w2, 2*M)[0]
    mX3 = dftAnal(xseg, w1, 2*M)[0]
    
    mX1_80 = mX1[:80]
    mX2_80 = mX2[:80]
    mX3_80 = mX3[:80]
    
    plt.plot(mX1_80, label="mX1(half/half)", color="red")
    plt.plot(mX2_80, label="mX1(full/full)", color="blue")
    plt.plot(mX3_80, label="mX1(half/full)", color="green")
    plt.show()
    
    return mX1_80, mX2_80, mX3_80

def window(sw, window = 'barthann'):
    """
    Creates 2d window of size sw
    
    --------------------------------------------------------------------------
    Usage:
    
    Call:  w = window(sw, window = 'barthann')
    
    Input: sw       size of window 
           window   string specifying window type

    Output: Window of size sw 
    --------------------------------------------------------------------------

    Copyright (C) 2010 Michael Hirsch
    """    

    w1  = signal.get_window(window,sw[0])
    w1  = (w1 + w1[::-1])/2
    w1 -= w1.min()
    w2  = signal.get_window(window,sw[1])
    w2  = (w2 + w2[::-1])/2
    w2 -= w2.min()

    www = np.outer(w1,w2)
    www = www/www.max()

    www = np.maximum(www, 1e-16)

    return www

 def test_window_external(self):
     x = np.zeros(16)
     x[0] = 1
     f, p = periodogram(x, 10, 'hann')
     win = signal.get_window('hann', 16)
     fe, pe = periodogram(x, 10, win)
     assert_array_almost_equal_nulp(p, pe)
     assert_array_almost_equal_nulp(f, fe)
     win_err = signal.get_window('hann', 32)
     assert_raises(ValueError, periodogram, x,
                   10, win_err)  # win longer than signal

    def window(self, index1=None, index2=None, is_positional=False, win_fcn='boxcar',
               fftbins=False):
        """
        Applies a window to the signal within a given time range.

        Parameters
        ----------
        index1 : float or int, optional
            The start index/position of the window. Default value is minimum of index.

        index2 : float or int, optional
            The end index/position of the window. Default value is maximum of index.

        is_positional : bool, optional
            Indicates whether the inputs `index1` and `index2` are positional or value
            based. Default is :const:`False`, i.e. value based.

        win_fcn : string/float/tuple, optional
            The type of window to create. See the function
            :func:`scipy.signal.get_window()` for a complete list of
            available windows, and how to pass extra parameters for a
            specific window function.

        fftbins : bool, optional
            If True, then applies a symmetric window with respect to index of value 0.

        Returns
        -------
        Signal:
            The windowed Signal signal.

        Note
        ----
          If the window requires no parameters, then `win_fcn` can be a string.
          If the window requires parameters, then `win_fcn` must be a tuple
          with the first argument the string name of the window, and the next
          arguments the needed parameters. If `win_fcn` is a floating point
          number, it is interpreted as the beta parameter of the kaiser window.
        """
        wind = Signal(0, index=self.index)
        if is_positional:
            if isinstance(index1, float) or isinstance(index2, float):
                raise ValueError('Indices are floats, are you sure you want positional indices?')
            index1 = wind.index.values[index1]
            index2 = wind.index.values[index2]

        wind[index1:index2] = get_window(win_fcn, len(wind[index1:index2]))
        if fftbins:
            if wind[-index2:-index1].size == 0:
                raise IndexError('The signal does not have values at the negative of the indices '
                                 'supplied. Disable fftbins for one-sided windowing.')
            wind[-index2:-index1] = get_window(win_fcn, len(wind[-index2:-index1]))
        return self*wind

def sineModelMultiResTest1(inputFile, M1, M2, M3):
    
    print "\n\n\n###############  RUN MULTIRESOLUTION TEST  ###############\n"
    
    #M1 = 4095
    #M2 = 2047
    #M3 = 1023
    
    print "M1: "
    print M1
    print "M2: "
    print M2
    print "M3: "
    print M3
    
    N1 = int(pow(2, np.ceil(np.log2(M1))))      # FFT Size, power of 2 larger than M
    N2 = int(pow(2, np.ceil(np.log2(M2))))      # FFT Size, power of 2 larger than M
    N3 = int(pow(2, np.ceil(np.log2(M3))))      # FFT Size, power of 2 larger than M
    
    print "N1: "
    print N1
    print "N2: "
    print N2
    print "N3: "
    print N3
    
    t = -80.0                                   # threshold

    fs, x = UF.wavread(inputFile)               # read input sound
    
    #print "Ploting \"x\""
    #plt.plot(x)
    
    window = 'blackman'                         # Window type
    w1 = get_window(window, M1)                 # compute analysis window
    w2 = get_window(window, M2)                 # compute analysis window
    w3 = get_window(window, M3)                 # compute analysis window
    
    #B1 = 1000    # the band from 0 to 999 Hz
    #B2 = 5000    # the band from 1000 to 4999 Hz
    #B3 = 22050   # the band from 5000 to 22049 Hz
    B1 = 100
    B2 = 1000
    B3 = 22050
    
    if (B3 > (44100 / 2)):
        B3 = 22050
    
    return sineModelMultiRes(x, fs, w1, w2, w3, N1, N2, N3, t, B1=1000, B2=5000, B3=22050)

def spectrum(data, sampling, NFFT=256, overlap=0.5,\
             window='hanning', detrender=mlab.detrend_linear,\
             sides='onesided', scale='PSD'):

  numpoints  = len(data)
  numoverlap = int(sampling * (1.0 - overlap))

  if isinstance(window,str):
    window=window.lower()

  win = signal.get_window(window, NFFT)

  # calculate PSD with given parameters
  spec,freq = mlab.psd(data, NFFT=NFFT, Fs=sampling, noverlap=numoverlap,\
                       window=win, sides=sides, detrend=detrender)

  # rescale data to meet user's request
  scale = scale.lower()
  if scale == 'asd':
    spec = numpy.sqrt(spec) * numpy.sqrt(2 / (sampling*sum(win**2)))
  elif scale == 'psd':
    spec *= 2/(sampling*sum(win**2))
  elif scale == 'as':
    spec = nump.sqrt(spec) * numpy.sqrt(2) / sum(win)
  elif scale == 'ps':
    spec = spec * 2 / (sum(win)**2)

  return freq, spec.flatten()

def computeEngEnv(inputFile, window, M, N, H):
    """
    Inputs:
            inputFile (string): input sound file (monophonic with sampling rate of 44100)
            window (string): analysis window type (choice of rectangular, triangular, hanning, 
                hamming, blackman, blackmanharris)
            M (integer): analysis window size (odd positive integer)
            N (integer): FFT size (power of 2, such that N > M)
            H (integer): hop size for the stft computation
    Output:
            The function should return a numpy array engEnv with shape Kx2, K = Number of frames
            containing energy envelop of the signal in decibles (dB) scale
            engEnv[:,0]: Energy envelope in band 0 < f < 3000 Hz (in dB)
            engEnv[:,1]: Energy envelope in band 3000 < f < 10000 Hz (in dB)
    """
    
    (fs,x) = UF.wavread(inputFile)
    w = get_window(window, M)
    (xmX, xpX) = stft.stftAnal(x, fs, w, N, H)

    kLow1 = 0
    
    kLow2 = 0
    while (True):
	kLow2 += 1
	if( (kLow2 < N*(fLow2)/float(fs)) & (kLow2 > N*(fLow2)/float(fs) - 1.0 ) ):
	    break
    
    kHigh1 = 0
    while (True):
	kHigh1 += 1
	if( (kHigh1 < N*(fHigh1)/float(fs)) & (kHigh1 > N*(fHigh1)/float(fs) - 1.0 ) ):
	    break
    
    kHigh2 = 0
    while (True):
	kHigh2 += 1
	if( (kHigh2 < N*(fHigh2)/float(fs)) & (kHigh2 > N*(fHigh2)/float(fs) - 1.0 ) ):
	    break
    
    nHops = int(xmX.shape[0])
    out = np.zeros((nHops,2))
    
    i = 0
    while i < nHops:
        subxmX = xmX[i,:]
    
        subLowxmX = subxmX[kLow1+1:kLow2+1]
        subLowxmX = 10**(subLowxmX/20)
        eSignalLow = sum(subLowxmX**2)
        out[i,0] = 10.0*np.log10(eSignalLow)

        subHighxmX = subxmX[kHigh1+1:kHigh2+1]
        subHighxmX = 10**(subHighxmX/20)
        eSignalHigh = sum(subHighxmX**2)
        out[i,1] = 10.0*np.log10(eSignalHigh)
    
        i += 1 

    return out

def suppressFreqDFTmodel(x, fs, N):
    """
    Inputs:
        x (numpy array) = input signal of length M (odd)
        fs (float) = sampling frequency (Hz)
        N (positive integer) = FFT size
    Outputs:
        The function should return a tuple (y, yfilt)
        y (numpy array) = Output of the dftSynth() without filtering (M samples long)
        yfilt (numpy array) = Output of the dftSynth() with filtering (M samples long)
    The first few lines of the code have been written for you, do not modify it. 
    """
    M = len(x)
    w = get_window("hamming", M)
    outputScaleFactor = sum(w)

    # get Magnitude and Phase Spectrum
    mX, pX = dftAnal(x, w, N)

    # generate output signal without filtering
    y = dftSynth(mX, pX, M) * outputScaleFactor

    # get bin number that is nearest to 70 Hz
    bin_number = int(70.0 / (fs / float(N))) + 1

    # do the 'filtering'
    for i in range(bin_number + 1):
        mX[i] = -120

    # generate the time signal after filtering
    yfilt = dftSynth(mX, pX, M) * outputScaleFactor

    return (y, yfilt)

def computeSNR(inputFile, window, M, N, H):
    """
    Input:
            inputFile (string): wav file name including the path 
            window (string): analysis window type (choice of rectangular, triangular, hanning, hamming, 
                    blackman, blackmanharris)
            M (integer): analysis window length (odd positive integer)
            N (integer): fft size (power of two, > M)
            H (integer): hop size for the stft computation
    Output:
            The function should return a python tuple of both the SNR values (SNR1, SNR2)
            SNR1 and SNR2 are floats.
    """

    #read from the file
    FS, x = UF.wavread(inputFile)

    w = get_window(window, M)
    #do a stft computation
    y = stft.stft(x, FS, w, N, H)

    #compute SNR over complete signal
    diff = y - x
    energy_signal = (y**2).sum()
    energy_noise = (diff**2).sum()
    SNR1 = 10 * np.log10(energy_signal/energy_noise)

    #compute SNR over sliced signal
    energy_signal_sliced = (y[M:-M]**2).sum()
    energy_noise_sliced = (diff[M:-M]**2).sum()
    SNR2 = 10 * np.log10(energy_signal_sliced/energy_noise_sliced)


    return (SNR1, SNR2)

def suppressFreqDFTmodel(x, fs, N):
    """
    Inputs:
        x (numpy array) = input signal of length M (odd)
        fs (float) = sampling frequency (Hz)
        N (positive integer) = FFT size
    Outputs:
        The function should return a tuple (y, yfilt)
        y (numpy array) = Output of the dftSynth() without filtering (M samples long)
        yfilt (numpy array) = Output of the dftSynth() with filtering (M samples long)
    The first few lines of the code have been written for you, do not modify it. 
    """
    M = len(x)
    w = get_window('hamming', M)
    outputScaleFactor = sum(w)
    
    mx,px = dftAnal(x,w,N)

    mx2 = mx.copy()
    mx2[:np.floor(70*N/fs)+1] = -120

    yone = dftSynth(mx,px,M)*sum(w)
    ytwo = dftSynth(mx2,px,M)*sum(w)

    return yone,ytwo

def computeSNR(inputFile, window, M, N, H):
    """
    Input:
            inputFile (string): wav file name including the path 
            window (string): analysis window type (choice of rectangular, triangular, hanning, hamming, 
                    blackman, blackmanharris)
            M (integer): analysis window length (odd positive integer)
            N (integer): fft size (power of two, > M)
            H (integer): hop size for the stft computation
    Output:
            The function should return a python tuple of both the SNR values (SNR1, SNR2)
            SNR1 and SNR2 are floats.
    """
    ## your code here

    def energy(X, k1, k2):
        X2 = np.power(X, 2)
        return np.sum(X2[k1:k2])

    fs, x = UF.wavread(inputFile)
    w = get_window(window, M)
    xsyn = stft.stft(x, fs, w, N, H)
    noise = np.subtract(xsyn, x)
    
    Esignal1 = energy(x, 0, len(x))
    Enoise1 = energy(noise, 0, len(noise))
    SNR1 = 10*np.log10(Esignal1/Enoise1)
    
    Esignal2 = energy(x, M+1, len(x)-M-1)
    Enoise2 = energy(noise, M+1, len(noise)-M-1)
    SNR2 = 10*np.log10(Esignal2/Enoise2)

    return SNR1, SNR2

def stft(data,Fs=100):
    """
    Input: 3 Axis Data
    Output: Power Representation of
    """
    # Remove the DC Component of the function
    data = noDC(data)
    wlen = Fs * 10
    segs = len(data) // wlen # integer divsion to return number of segments of data
    windsegs = []
    numz = nextpow2(wlen) - wlen
    j = 0
    win = []
    for i in data:
        if j < wlen:  # append the value from j=0 to j = (number of samples - 1)
            win.append(i)
        else:
            j = 0
            for i in range(0, numz):  # Zero padding the window
                win.append(0)
            windsegs.append(win)  # Add that window to the segmented dataset
            win = [i]  # Reset the window variable
        j += 1
    b, a = signal.butter(4, [.01, .5], 'bandpass')
    dft=[]
    winlen = int(len(windsegs[0]))
    for seg in windsegs:
        seg = signal.lfilter(b, a, seg)
        wind = signal.get_window(('kaiser', 4.0), winlen)  # Beta Parameter of Kaiser: 4. Num Samples in window: 9.
        snip = seg * wind
        nfft = nextpow2(wlen)
        A = np.fft.fft(snip, nfft)
        dft.append(A)
    return dft

def extractMainLobe(window, M):
    """
    Input:
      window (string): Window type to be used (Either rectangular (‘boxcar’),
          ‘hamming’ or ‘blackmanharris’)
      M (integer): length of the window to be used
    Output:
      The function should return a numpy array containing the main lobe of
      the magnitude spectrum of the window in decibels (dB).
    """
    M = np.asarray(M, int)
    w = get_window(window, M)         # get the window
    N = 8 * M
    W_N = fftshift(fft(w, N))
    W_N = np.abs(W_N)
    W_N[W_N < eps] = eps
    idx_l = int(N / 2)
    lim_l = lim_u = False
    while lim_l is False:
        if W_N[idx_l-1] > W_N[idx_l] and W_N[idx_l+1] > W_N[idx_l]:
            lim_l = True
            print(idx_l)
        else:
            idx_l += -1
    idx_u = int(N / 2)+1
    while lim_u is False:
        if W_N[idx_u-1] > W_N[idx_u] and W_N[idx_u+1] > W_N[idx_u]:
            lim_u = True
            print(idx_u)
        else:
            idx_u += 1
    W_N_dB = 10 * np.log10(W_N)
    # return(W_N_dB, idx_l, idx_u)
    return(W_N_dB[idx_l:idx_u + 1])

def FIR(lowcutoff,highcutoff,fs):

	transwidthratio = 0.25

	fnyq = fs/2
	if highcutoff < lowcutoff :
		temp = highcutoff 
		highcutoff = lowcutoff
		lowcutoff = temp
		print "cutoff frequecy is swapped. lowcutoff cannot be higher than highcutoff"
                
	if fnyq-highcutoff < lowcutoff:
		print "ERROR : (highcutoff - half of sampling frequency (Nyquist)) cannot be higher than lowcutoff"
		return 0

	maxTWBArray = [lowcutoff, fnyq-highcutoff]

	maxDf = min(maxTWBArray)

	df = min([max([maxTWBArray[0] * transwidthratio, 2.0]),maxDf])

	filtorder = 3.3 / (df / fs) # Hamming window
	filtorder = np.ceil(filtorder / 2) * 2  + 1 # Filter order must be even. +1 for centerpoint window
	#print filtorder
	#print np.array([range(int(filtorder))])

	win = get_window('hamming',filtorder)

	cutoffarray = np.array([lowcutoff,highcutoff]) + [-df/2,df/2] # bandpass

	winsinc = firws(filtorder, cutoffarray / fnyq, win)

	return winsinc

def computeSNR(inputFile, window, M, N, H):
    """
    Input:
            inputFile (string): wav file name including the path 
            window (string): analysis window type (choice of rectangular, triangular, hanning, hamming, 
                    blackman, blackmanharris)
            M (integer): analysis window length (odd positive integer)
            N (integer): fft size (power of two, > M)
            H (integer): hop size for the stft computation
    Output:
            The function should return a python tuple of both the SNR values (SNR1, SNR2)
            SNR1 and SNR2 are floats.
    """
    ## your code here
    fs, x = UF.wavread(inputFile)

    w = get_window(window, M)

    y = stft.stft(x, fs, w, N, H)
    noise = np.array(x - y)

    E_x = np.sum( abs(x)**2 )
    E_noise = np.sum( abs(noise)**2 )

    E_xAfterM = np.sum( abs( x[M : x.size-M] )**2 )
    E_nAfterM = np.sum( abs( noise[M : x.size-M] )**2 )

    SNR1 = 10 * np.log10(E_x / E_noise)
    SNR2 = 10 * np.log10(E_xAfterM/E_nAfterM)

    return (SNR1, SNR2)