def THDN(signal, sample_rate, filename):

    signal = signal - mean(signal)  # this is so that the DC offset is not the peak in the case of PDM
    windowed = signal * blackmanharris(len(signal)) 

    

    # Find the peak of the frequency spectrum (fundamental frequency), and filter 
    # the signal by throwing away values between the nearest local minima
    f = rfft(windowed)

    #limit the bandwidth
    if len(f) > 24000:
       f[24000:len(f)] = 0

    bandwidth_limited = irfft(f)
    total_rms = rms_flat(bandwidth_limited)

    #for i in range(6000):
       #print abs(f[i])
    i = argmax(abs(f))

    #This will be exact if the frequency under test falls into the middle of an fft bin 
    print 'Frequency: %f Hz' % (sample_rate * (i / len(windowed)))
    lowermin, uppermin = find_range(abs(f), i)


    #notch out the signal
    f[lowermin: uppermin] = 0



    noise = irfft(f)
    THDN = rms_flat(noise) / total_rms
    print "THD+N:     %.12f%% or %.12f dB" % (THDN * 100, 20 * log10(THDN))

def freq_from_HPS(sig, fs):
    """
    Estimate frequency using harmonic product spectrum (HPS)
    
    """
    windowed = sig * blackmanharris(len(sig))

    from pylab import subplot, plot, log, copy, show

    #harmonic product spectrum:
    c = abs(rfft(windowed))
    maxharms = 8
    subplot(maxharms,1,1)
    plot(log(c))
    for x in range(2,maxharms):
        a = copy(c[::x]) #Should average or maximum instead of decimating
        # max(c[::x],c[1::x],c[2::x],...)
        c = c[:len(a)]
        i = argmax(abs(c))
        true_i = parabolic(abs(c), i)[0]
        print 'Pass %d: %f Hz' % (x, fs * true_i / len(windowed))
        c *= a
        subplot(maxharms,1,x)
        plot(log(c))
    show()

def window(f,start,stop,type='blackman'):
    """
    runs the data through a hamming window.
    @param f: The data matrix
    @param start: The start index of the hamming window.
    @param stop: The end index of the hamming window.
    """
    h=numpy.zeros(f.shape,dtype=float)

    if len(h.shape)==1:
        if type=='hamming':
            h[start:stop]=signal.hamming(stop-start)
        elif type=='blackman':
            h[start:stop]=signal.blackman(stop-start)
        elif type=='hann':
            h[start:stop]=signal.hann(stop-start)
        elif type=='blackmanharris':
            h[start:stop]=signal.blackmanharris(stop-start)
        elif type=='rectangular' or type=='rect' or type=='boxcar':
            h[start:stop]=signal.boxcar(stop-start)
    else:
        if type=='hamming':
            h[:,start:stop]=signal.hamming(stop-start)
        elif type=='blackman':
            h[:,start:stop]=signal.blackman(stop-start)
        elif type=='hann':
            h[:,start:stop]=signal.hann(stop-start)
        elif type=='blackmanharris':
            h[:,start:stop]=signal.blackmanharris(stop-start)
        elif type=='rectangular' or type=='rect' or type=='boxcar':
            h[:,start:stop]=signal.boxcar(stop-start)
    return numpy.multiply(f,h)

def freq_from_fft(signal, fs):
    # Compute Fourier transform of windowed signal
    windowed = signal * blackmanharris(len(signal))
    f = rfft(windowed)
    i = argmax(abs(f))
    
    # Convert to equivalent frequency
    return fs*i/len(windowed)

def delayTransformVisibilities(visList,df,kernel=None,wind='Blackman-Harris'):
    '''
     ARGS:
     visList: nfreqxnvis complex matrix of visibilities
     df: float indicating channel width
     RETURNS:
     ndelayxnvis grid of delay-transformed visibilities
    '''
    nf=visList.shape[0]
    if kernel is None:
        kernel=np.ones(nf)

    nv=visList.shape[1]
    windowGrid=np.zeros_like(visList)
    if wind=='Blackman-Harris':
        window=signal.blackmanharris(nf)
    elif wind=='Nithya':
        window=signal.blackmanharris(nf/2)
        window=signal.fftconvolve(window,window,mode='full')
        window=np.append(window,window[-1])
    window/=window.max()
    window/=np.sqrt(np.mean(np.abs(kernel*window)**2.)) 
    for vNum in range(nv):
        windowGrid[:,vNum]=window*kernel
    return df*fft.fftshift(fft.fft(fft.fftshift(visList*windowGrid,axes=[0]),axis=0),axes=[0])

 def test_basic(self):
     assert_allclose(signal.blackmanharris(6, False),
                     [6.0e-05, 0.055645, 0.520575, 1.0, 0.520575, 0.055645])
     assert_allclose(signal.blackmanharris(6),
                     [6.0e-05, 0.1030114893456638, 0.7938335106543362,
                      0.7938335106543364, 0.1030114893456638, 6.0e-05])
     assert_allclose(signal.blackmanharris(7, sym=True),
                     [6.0e-05, 0.055645, 0.520575, 1.0, 0.520575, 0.055645,
                      6.0e-05])

 def test_basic(self):
     assert_allclose(signal.blackmanharris(6, False),
                     [6.0e-05, 0.055645, 0.520575, 1.0, 0.520575, 0.055645])
     assert_allclose(signal.blackmanharris(7, sym=False),
                     [6.0e-05, 0.03339172347815117, 0.332833504298565,
                      0.8893697722232837, 0.8893697722232838,
                      0.3328335042985652, 0.03339172347815122])
     assert_allclose(signal.blackmanharris(6),
                     [6.0e-05, 0.1030114893456638, 0.7938335106543362,
                      0.7938335106543364, 0.1030114893456638, 6.0e-05])
     assert_allclose(signal.blackmanharris(7, sym=True),
                     [6.0e-05, 0.055645, 0.520575, 1.0, 0.520575, 0.055645,
                      6.0e-05])

def sineModelSynth(tfreq, tmag, tphase, N, H, fs):
    """
    Synthesis of a sound using the sinusoidal model
    tfreq, tmag, tphase: frequencies, magnitudes and phases of sinusoids
    N: synthesis FFT size, H: hop size, fs: sampling rate
    returns y: output array sound
    """

    hN = N / 2 # half of FFT size for synthesis
    L = tfreq.shape[0] # number of frames
    pout = 0
    ysize = H * (L+3)
    y = np.zeros(ysize)
    # synthesis window
    sw = np.zeros(N)
    ow = triang(2 * H)
    sw[hN-H:hN+H] = ow
    bh = blackmanharris(N)
    bh = bh / sum(bh)
    sw[hN-H:hN+H] = sw[hN-H:hN+H] / bh[hN-H:hN+H]
    lastytfreq = tfreq[0,:]
    ytphase = 2*np.pi*np.random.rand(tfreq[0,:].size)
    for l in range(L): # iterate over all frames
        if (tphase.size > 0):
            ytphase = tphase[l, :]

def freq_power_from_fft(sig, fs):
	# Compute Fourier transform of windowed signal
	windowed = sig * blackmanharris(CHUNK)
	fftData = abs(rfft(windowed))
	# powerData = 20*log10(fftData)
	# Find the index of the peak and interpolate to get a more accurate peak
	i = argmax(fftData) # Just use this for less-accurate, naive version
	# make sure i is not an endpoint of the bin
	if i != len(fftData)-1 and i != 0:
		#print("data: " + str(parabolic(log(abs(fftData)), i)))
		true_i,logmag = parabolic(log(abs(fftData)), i)
		# Convert to equivalent frequency
		freq=  fs * true_i / len(windowed)
		#print("frequency="+ str(freq) + " log of magnitude=" + str(logmag))
		if logmag < 0:
			logmag = 0	
		return freq,logmag
	
	else:
		freq = fs * i / len(windowed)
		logmag = log(abs(fftData))[i]
		if logmag < 0:
			logmag = 0
		#print("frequency="+ str(freq) + "log of magnitude not interp=" + str(logmag))
		return freq,logmag

    def run(self, tf, dt, K, c1, c2):
        """ Run a simulation """
        n, m, k = self.n, self.m, self.k

        # Total simulation time
        simTime = int(tf / dt)

        # Returns the three synaptic connections kernels
        W12, W21, W22, delays = self.build_kernels(K)

        # Compute delays by dividing distances by axonal velocity
        delays12 = np.floor(delays[0] / c2)
        delays21 = np.floor(delays[1] / c1)
        delays22 = np.floor(delays[2] / c2)
        maxDelay = int(max(delays12[0].max(), delays21[0].max(), delays22[0].max()))

        # Set the initial conditions and the history
        self.initial_conditions(simTime)

        # Initialize the cortical and striatal inputs
        Cx = 0.5
        Str = 0.4

        # Presynaptic activities
        pre12, pre21, pre22 = np.empty((m,)), np.empty((m,)), np.empty((m,))

        # Simulation
        for i in range(maxDelay, simTime):
            # Take into account the history of rate for each neuron according
            # to its axonal delay
            for idxi, ii in enumerate(range(m)):
                mysum = 0.0
                for jj in range(k, n):
                    mysum += (W12[ii, jj] * self.X2[i - delays12[ii, jj], jj]) * self.dx
                pre12[idxi] = mysum

            for idxi, ii in enumerate(range(k, n)):
                mysum = 0.0
                for jj in range(0, m):
                    mysum += (W21[ii, jj] * self.X1[i - delays21[ii, jj], jj]) * self.dx
                pre21[idxi] = mysum

            for idxi, ii in enumerate(range(k, n)):
                mysum = 0.0
                for jj in range(k, n):
                    mysum += (W22[ii, jj] * self.X2[i - delays22[ii, jj], jj]) * self.dx
                pre22[idxi] = mysum

            # Forward Euler step
            self.X1[i, :m] = self.X1[i - 1, :m] + (-self.X1[i - 1, :m] + self.S1(-pre12 + Cx)) * dt / self.tau1
            self.X2[i, k:] = self.X2[i - 1, k:] + (-self.X2[i - 1, k:] + self.S2(pre21 - pre22 - Str)) * dt / self.tau2
        dx = 1.0 / float(m)
        fr = self.X1.sum(axis=1) * dx / 1.0

        signal = detrend(fr)
        windowed = signal * blackmanharris(len(signal))
        f = rfft(windowed)
        i = np.argmax(np.abs(f))
        # true_i = parabolic(np.log(np.abs(f)), i)[0]
        return i

def _BLACKMANHARRIS_Window(signal):
    """ blackman-harris window process
        Use Windows to smooth transient data in the begin/end
    """
    windowed = signal * blackmanharris(len(signal))  # TODO Kaiser?

    return windowed

def stochasticResidualAnal(x, N, H, sfreq, smag, sphase, fs, stocf):
	"""
	Subtract sinusoids from a sound and approximate the residual with an envelope
	x: input sound, N: fft size, H: hop-size
	sfreq, smag, sphase: sinusoidal frequencies, magnitudes and phases
	fs: sampling rate; stocf: stochastic factor, used in the approximation
	returns stocEnv: stochastic approximation of residual
	"""

	hN = N/2                                              # half of fft size
	x = np.append(np.zeros(hN),x)                         # add zeros at beginning to center first window at sample 0
	x = np.append(x,np.zeros(hN))                         # add zeros at the end to analyze last sample
	bh = blackmanharris(N)                                # synthesis window
	w = bh/ sum(bh)                                       # normalize synthesis window
	L = sfreq.shape[0]                                    # number of frames, this works if no sines
	pin = 0
	for l in range(L):
		xw = x[pin:pin+N] * w                               # window the input sound
		X = fft(fftshift(xw))                               # compute FFT
		Yh = UF_C.genSpecSines(N*sfreq[l,:]/fs, smag[l,:], sphase[l,:], N)   # generate spec sines
		Xr = X-Yh                                           # subtract sines from original spectrum
		mXr = 20*np.log10(abs(Xr[:hN]))                     # magnitude spectrum of residual
		mXrenv = resample(np.maximum(-200, mXr), mXr.size*stocf)  # decimate the mag spectrum
		if l == 0:                                          # if first frame
			stocEnv = np.array([mXrenv])
		else:                                               # rest of frames
			stocEnv = np.vstack((stocEnv, np.array([mXrenv])))
		pin += H                                            # advance sound pointer
	return stocEnv

def sineSubtraction(x, N, H, sfreq, smag, sphase, fs):
	"""
	Subtract sinusoids from a sound
	x: input sound, N: fft-size, H: hop-size
	sfreq, smag, sphase: sinusoidal frequencies, magnitudes and phases
	returns xr: residual sound
	"""

	hN = N/2                                           # half of fft size
	x = np.append(np.zeros(hN),x)                      # add zeros at beginning to center first window at sample 0
	x = np.append(x,np.zeros(hN))                      # add zeros at the end to analyze last sample
	bh = blackmanharris(N)                             # blackman harris window
	w = bh/ sum(bh)                                    # normalize window
	sw = np.zeros(N)                                   # initialize synthesis window
	sw[hN-H:hN+H] = triang(2*H) / w[hN-H:hN+H]         # synthesis window
	L = sfreq.shape[0]                                 # number of frames, this works if no sines
	xr = np.zeros(x.size)                              # initialize output array
	pin = 0
	for l in range(L):
		xw = x[pin:pin+N]*w                              # window the input sound
		X = fft(fftshift(xw))                            # compute FFT
		Yh = UF_C.genSpecSines(N*sfreq[l,:]/fs, smag[l,:], sphase[l,:], N)   # generate spec sines
		Xr = X-Yh                                        # subtract sines from original spectrum
		xrw = np.real(fftshift(ifft(Xr)))                # inverse FFT
		xr[pin:pin+N] += xrw*sw                          # overlap-add
		pin += H                                         # advance sound pointer
	xr = np.delete(xr, range(hN))                      # delete half of first window which was added in stftAnal
	xr = np.delete(xr, range(xr.size-hN, xr.size))     # delete half of last window which was added in stftAnal
	return xr

def sinewaveSynth(freq, mag, N, H, fs):
  # Synthesis of a time-varying sinusoid
  # freq,mag, phase: frequency, magnitude and phase of sinusoid,
  # N: synthesis FFT size, H: hop size, fs: sampling rate
  # returns y: output array sound
  hN = N/2                                                # half of FFT size for synthesis
  L = freq.size                                           # number of frames
  pout = 0                                                # initialize output sound pointer         
  ysize = H*(L+3)                                         # output sound size
  y = np.zeros(ysize)                                     # initialize output array
  sw = np.zeros(N)                                        # initialize synthesis window
  ow = triang(2*H);                                       # triangular window
  sw[hN-H:hN+H] = ow                                      # add triangular window
  bh = blackmanharris(N)                                  # blackmanharris window
  bh = bh / sum(bh)                                       # normalized blackmanharris window
  sw[hN-H:hN+H] = sw[hN-H:hN+H]/bh[hN-H:hN+H]             # normalized synthesis window
  lastfreq = freq[0]                                      # initialize synthesis frequencies
  phase = 0                                               # initialize synthesis phases 
  for l in range(L):                                      # iterate over all frames
    phase += (np.pi*(lastfreq+freq[l])/fs)*H            # propagate phases
    Y = UF.genSpecSines(freq[l], mag[l], phase, N, fs)    # generate sines in the spectrum         
    lastfreq = freq[l]                                    # save frequency for phase propagation
    yw = np.real(fftshift(ifft(Y)))                       # compute inverse FFT
    y[pout:pout+N] += sw*yw                               # overlap-add and apply a synthesis window
    pout += H                                             # advance sound pointer
  y = np.delete(y, range(hN))                             # delete half of first window
  y = np.delete(y, range(y.size-hN, y.size))              # delete half of the last window 
  return y

def plot_specgram(ax, data, fs, nfft=256, noverlap=128, window='hann',
                  cmap='jet', interpolation='bilinear', rasterized=True):

    if window not in SPECGRAM_WINDOWS:
        raise ValueError("Window not supported")

    elif window == "boxcar":
        mwindow = signal.boxcar(nfft)
    elif window == "hamming":
        mwindow = signal.hamming(nfft)
    elif window == "hann":
        mwindow = signal.hann(nfft)
    elif window == "bartlett":
        mwindow = signal.bartlett(nfft)
    elif window == "blackman":
        mwindow = signal.blackman(nfft)
    elif window == "blackmanharris":
        mwindow = signal.blackmanharris(nfft)

    specgram, freqs, time = mlab.specgram(data, NFFT=nfft, Fs=fs,
                                          window=mwindow,
                                          noverlap=noverlap)
    specgram = 10 * np.log10(specgram[1:, :])
    specgram = np.flipud(specgram)

    freqs = freqs[1:]
    halfbin_time = (time[1] - time[0]) / 2.0
    halfbin_freq = (freqs[1] - freqs[0]) / 2.0
    extent = (time[0] - halfbin_time, time[-1] + halfbin_time,
              freqs[0] - halfbin_freq, freqs[-1] + halfbin_freq)

    ax.imshow(specgram, cmap=cmap, interpolation=interpolation,
                            extent=extent, rasterized=rasterized)
    ax.axis('tight')

def stochasticResidual(x, N, H, sfreq, smag, sphase, fs, stocf):
  # subtract sinusoids from a sound
  # x: input sound, N: fft-size, H: hop-size
  # sfreq, smag, sphase: sinusoidal frequencies, magnitudes and phases
  # returns mYst: stochastic approximation of residual 
  hN = N/2  
  x = np.append(np.zeros(hN),x)                    # add zeros at beginning to center first window at sample 0
  x = np.append(x,np.zeros(hN))                    # add zeros at the end to analyze last sample
  bh = blackmanharris(N)                           # synthesis window
  w = bh/ sum(bh)                                  # normalize synthesis window
  L = sfreq[:,0].size                              # number of frames   
  pin = 0
  for l in range(L):
    xw = x[pin:pin+N]*w                            # window the input sound                               
    X = fft(fftshift(xw))                          # compute FFT 
    Yh = UF_C.genSpecSines(N*sfreq[l,:]/fs, smag[l,:], sphase[l,:], N)   # generate spec sines          
    Xr = X-Yh                                      # subtract sines from original spectrum
    mXr = 20*np.log10(abs(Xr[:hN]))                # magnitude spectrum of residual
    mXrenv = resample(np.maximum(-200, mXr), mXr.size*stocf)  # decimate the mag spectrum                        
    if l == 0: 
      mYst = np.array([mXrenv])
    else:
      mYst = np.vstack((mYst, np.array([mXrenv])))
    pin += H                                       # advance sound pointer
  return mYst

def sineModelSynth(tfreq, tmag, tphase, N, H, fs):
	"""
	Synthesis of a sound using the sinusoidal model
	tfreq,tmag,tphase: frequencies, magnitudes and phases of sinusoids
	N: synthesis FFT size, H: hop size, fs: sampling rate
	returns y: output array sound
	"""
	
	hN = N/2                                                # half of FFT size for synthesis
	L = tfreq.shape[0]                                      # number of frames
	pout = 0                                                # initialize output sound pointer         
	ysize = H*(L+3)                                         # output sound size
	y = np.zeros(ysize)                                     # initialize output array
	sw = np.zeros(N)                                        # initialize synthesis window
	ow = triang(2*H)                                        # triangular window
	sw[hN-H:hN+H] = ow                                      # add triangular window
	bh = blackmanharris(N)                                  # blackmanharris window
	bh = bh / sum(bh)                                       # normalized blackmanharris window
	sw[hN-H:hN+H] = sw[hN-H:hN+H]/bh[hN-H:hN+H]             # normalized synthesis window
	lastytfreq = tfreq[0,:]                                 # initialize synthesis frequencies
	ytphase = 2*np.pi*np.random.rand(tfreq[0,:].size)       # initialize synthesis phases 
	for l in range(L):                                      # iterate over all frames
		if (tphase.size > 0):                                 # if no phases generate them
			ytphase = tphase[l,:] 
		else:
			ytphase += (np.pi*(lastytfreq+tfreq[l,:])/fs)*H     # propagate phases
		Y = UF.genSpecSines(tfreq[l,:], tmag[l,:], ytphase, N, fs)  # generate sines in the spectrum         
		lastytfreq = tfreq[l,:]                               # save frequency for phase propagation
		ytphase = ytphase % (2*np.pi)                         # make phase inside 2*pi
		yw = np.real(fftshift(ifft(Y)))                       # compute inverse FFT
		y[pout:pout+N] += sw*yw                               # overlap-add and apply a synthesis window
		pout += H                                             # advance sound pointer
	y = np.delete(y, range(hN))                             # delete half of first window
	y = np.delete(y, range(y.size-hN, y.size))              # delete half of the last window 
	return y

def waveform_to_file(station,clen=10000,oversample=10, filt=False, outpath=None,verbose=False):
    """
    lets use 0.1 s code cycle and coherence assumption
    our transmit bandwidth is 100 kHz, and with a clen=10e3 baud code,
    that is 0.1 seconds per cycle as a coherence assumption.
    furthermore, we use a 1 MHz bandwidth, so we oversample by a factor of 10.

    NOTE: this writing method doesn't store any metadata - have to know the sample rate
    """

    a = rep_seq(create_pseudo_random_code(clen,station,verbose), rep=oversample)

    if filt == True:
        w = np.zeros([oversample*clen], dtype='complex64') # yes, zeros for zero-padded
        fl = int(oversample+(0.1*oversample))

        w[:fl]= signal.blackmanharris(fl) # W[fl:] \equiv 0

        aa = np.fft.ifft(np.fft.fft(w) * np.fft.fft(a))

        a = (aa/abs(aa).max()).astype('complex64') #normalized, single prec complex

    if outpath:
        p = Path(outpath).expanduser()
        p.mkdir(parents=True, exist_ok=True)

        ofn = p / f"code-l{clen}-b{oversample}-{station:06d}.bin"
        print('writing',ofn)

        # https://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.tofile.html
        a.tofile(ofn) # this binary format is OK for GNU Radio to read

    return a

def harmonicModel(x, fs, w, N, t, nH, minf0, maxf0, f0et):
	"""
	Analysis/synthesis of a sound using the sinusoidal harmonic model
	x: input sound, fs: sampling rate, w: analysis window, 
	N: FFT size (minimum 512), t: threshold in negative dB, 
	nH: maximum number of harmonics, minf0: minimum f0 frequency in Hz, 
	maxf0: maximim f0 frequency in Hz, 
	f0et: error threshold in the f0 detection (ex: 5),
	returns y: output array sound
	"""

	hN = N/2                                                # size of positive spectrum
	hM1 = int(math.floor((w.size+1)/2))                     # half analysis window size by rounding
	hM2 = int(math.floor(w.size/2))                         # half analysis window size by floor
	x = np.append(np.zeros(hM2),x)                          # add zeros at beginning to center first window at sample 0
	x = np.append(x,np.zeros(hM1))                          # add zeros at the end to analyze last sample
	Ns = 512                                                # FFT size for synthesis (even)
	H = Ns/4                                                # Hop size used for analysis and synthesis
	hNs = Ns/2      
	pin = max(hNs, hM1)                                     # init sound pointer in middle of anal window          
	pend = x.size - max(hNs, hM1)                           # last sample to start a frame
	fftbuffer = np.zeros(N)                                 # initialize buffer for FFT
	yh = np.zeros(Ns)                                       # initialize output sound frame
	y = np.zeros(x.size)                                    # initialize output array
	w = w / sum(w)                                          # normalize analysis window
	sw = np.zeros(Ns)                                       # initialize synthesis window
	ow = triang(2*H)                                        # overlapping window
	sw[hNs-H:hNs+H] = ow      
	bh = blackmanharris(Ns)                                 # synthesis window
	bh = bh / sum(bh)                                       # normalize synthesis window
	sw[hNs-H:hNs+H] = sw[hNs-H:hNs+H] / bh[hNs-H:hNs+H]     # window for overlap-add
	hfreqp = []
	f0t = 0
	f0stable = 0
	while pin<pend:             
	#-----analysis-----             
		x1 = x[pin-hM1:pin+hM2]                               # select frame
		mX, pX = DFT.dftAnal(x1, w, N)                        # compute dft
		ploc = UF.peakDetection(mX, t)                        # detect peak locations     
		iploc, ipmag, ipphase = UF.peakInterp(mX, pX, ploc)   # refine peak values
		ipfreq = fs * iploc/N
		f0t = UF.f0Twm(ipfreq, ipmag, f0et, minf0, maxf0, f0stable)  # find f0
		if ((f0stable==0)&(f0t>0)) \
				or ((f0stable>0)&(np.abs(f0stable-f0t)<f0stable/5.0)):
			f0stable = f0t                                     # consider a stable f0 if it is close to the previous one
		else:
			f0stable = 0
		hfreq, hmag, hphase = harmonicDetection(ipfreq, ipmag, ipphase, f0t, nH, hfreqp, fs) # find harmonics
		hfreqp = hfreq
	#-----synthesis-----
		Yh = UF.genSpecSines(hfreq, hmag, hphase, Ns, fs)     # generate spec sines          
		fftbuffer = np.real(ifft(Yh))                         # inverse FFT
		yh[:hNs-1] = fftbuffer[hNs+1:]                        # undo zero-phase window
		yh[hNs-1:] = fftbuffer[:hNs+1] 
		y[pin-hNs:pin+hNs] += sw*yh                           # overlap-add
		pin += H                                              # advance sound pointer
	y = np.delete(y, range(hM2))                            # delete half of first window which was added in stftAnal
	y = np.delete(y, range(y.size-hM1, y.size))             # add zeros at the end to analyze last sample
	return y

def sprModel(x, fs, w, N, t):
	"""
	Analysis/synthesis of a sound using the sinusoidal plus residual model, one frame at a time
	x: input sound, fs: sampling rate, w: analysis window, 
	N: FFT size (minimum 512), t: threshold in negative dB, 
	returns y: output sound, ys: sinusoidal component, xr: residual component
	"""

	hN = N/2                                                      # size of positive spectrum
	hM1 = int(math.floor((w.size+1)/2))                           # half analysis window size by rounding
	hM2 = int(math.floor(w.size/2))                               # half analysis window size by floor
	Ns = 512                                                      # FFT size for synthesis (even)
	H = Ns/4                                                      # Hop size used for analysis and synthesis
	hNs = Ns/2      
	pin = max(hNs, hM1)                                           # initialize sound pointer in middle of analysis window          
	pend = x.size - max(hNs, hM1)                                 # last sample to start a frame
	fftbuffer = np.zeros(N)                                       # initialize buffer for FFT
	ysw = np.zeros(Ns)                                            # initialize output sound frame
	xrw = np.zeros(Ns)                                            # initialize output sound frame
	ys = np.zeros(x.size)                                         # initialize output array
	xr = np.zeros(x.size)                                         # initialize output array
	w = w / sum(w)                                                # normalize analysis window
	sw = np.zeros(Ns)     
	ow = triang(2*H)                                              # overlapping window
	sw[hNs-H:hNs+H] = ow      
	bh = blackmanharris(Ns)                                       # synthesis window
	bh = bh / sum(bh)                                             # normalize synthesis window
	wr = bh                                                       # window for residual
	sw[hNs-H:hNs+H] = sw[hNs-H:hNs+H] / bh[hNs-H:hNs+H]
	while pin<pend:  
  #-----analysis-----             
		x1 = x[pin-hM1:pin+hM2]                                     # select frame
		mX, pX = DFT.dftAnal(x1, w, N)                              # compute dft
		ploc = UF.peakDetection(mX, t)                              # find peaks 
		iploc, ipmag, ipphase = UF.peakInterp(mX, pX, ploc)         # refine peak values		iploc, ipmag, ipphase = UF.peakInterp(mX, pX, ploc)          # refine peak values
		ipfreq = fs*iploc/float(N)                                  # convert peak locations to Hertz
		ri = pin-hNs-1                                              # input sound pointer for residual analysis
		xw2 = x[ri:ri+Ns]*wr                                        # window the input sound                                       
		fftbuffer = np.zeros(Ns)                                    # reset buffer
		fftbuffer[:hNs] = xw2[hNs:]                                 # zero-phase window in fftbuffer
		fftbuffer[hNs:] = xw2[:hNs]                           
		X2 = fft(fftbuffer)                                         # compute FFT for residual analysis
  #-----synthesis-----
		Ys = UF.genSpecSines(ipfreq, ipmag, ipphase, Ns, fs)        # generate spec of sinusoidal component          
		Xr = X2-Ys;                                                 # get the residual complex spectrum
		fftbuffer = np.zeros(Ns)
		fftbuffer = np.real(ifft(Ys))                               # inverse FFT of sinusoidal spectrum
		ysw[:hNs-1] = fftbuffer[hNs+1:]                             # undo zero-phase window
		ysw[hNs-1:] = fftbuffer[:hNs+1] 
		fftbuffer = np.zeros(Ns)
		fftbuffer = np.real(ifft(Xr))                               # inverse FFT of residual spectrum
		xrw[:hNs-1] = fftbuffer[hNs+1:]                             # undo zero-phase window
		xrw[hNs-1:] = fftbuffer[:hNs+1]
		ys[ri:ri+Ns] += sw*ysw                                      # overlap-add for sines
		xr[ri:ri+Ns] += sw*xrw                                      # overlap-add for residual
		pin += H                                                    # advance sound pointer
	y = ys+xr                                                     # sum of sinusoidal and residual components
	return y, ys, xr