def iff_filter(sig, scale, plot_show = 0):
    
    order = max(sig.size*scale,90)
    #order = 80
    # Extend signal on both sides for removing boundary effect in convolution
    sig_extend = np.ones(sig.size+int(order/2)*2)
    sig_extend[int(order/2):(sig.size+int(order/2))] = sig
    sig_extend[0:int(order/2)] = sig[(sig.size-int(order/2)):sig.size]
    sig_extend[(sig.size+int(order/2)):sig_extend.size] = sig[0:int(order/2)]
    
    # convolve with hamming window and normalize
    smooth_sig = np.convolve(sig_extend,np.hamming(order),'same')
    smooth_sig = smooth_sig[int(order/2):(sig.size+int(order/2))]
    smooth_sig = np.amax(sig)/np.amax(smooth_sig)*smooth_sig

    # Plot signal for debug
    if(plot_show == 1):
        fig, ax = plt.subplots(ncols=2)
        ax[0].plot(sig)
        ax[0].plot(smooth_sig,'-r')
        ax[0].plot(med_sig,'black')
        ax[1].loglog(rfft(sig))
        ax[1].loglog(rfft(smooth_sig),'-r')
        ax[1].loglog(rfft(med_sig),'black')
        plt.show()
        
    return smooth_sig

 def subtract_original_signal_from_picked_signal(self, original_signal, picked_signal):
     # Note this function assumes that the signals are aligned for the starting point!
     fft_length = max(len(original_signal), len(picked_signal))
     original_f_domain = rfft(original_signal, n= fft_length)
     picked_f_domain = rfft(picked_signal, n= fft_length)
     assert len(original_f_domain) == len(picked_f_domain)
     difference_signal = picked_f_domain - original_f_domain
     return irfft(difference_signal)

 def test_definition(self):
     x = [1,2,3,4,1,2,3,4]
     y = rfft(x)
     y1 = direct_rdft(x)
     assert_array_almost_equal(y,y1)
     x = [1,2,3,4,1,2,3,4,5]
     y = rfft(x)
     y1 = direct_rdft(x)
     assert_array_almost_equal(y,y1)

 def test_random_real(self):
     for size in [1, 51, 111, 100, 200, 64, 128, 256, 1024]:
         x = random([size]).astype(self.rdt)
         y1 = irfft(rfft(x))
         y2 = rfft(irfft(x))
         assert_equal(y1.dtype, self.rdt)
         assert_equal(y2.dtype, self.rdt)
         assert_array_almost_equal(y1, x, decimal=self.ndec, err_msg="size=%d" % size)
         assert_array_almost_equal(y2, x, decimal=self.ndec, err_msg="size=%d" % size)

 def test_random_real(self):
     for size in [1,51,111,100,200,64,128,256,1024]:
         x = random([size]).astype(self.rdt)
         y1 = irfft(rfft(x))
         y2 = rfft(irfft(x))
         self.failUnless(y1.dtype == self.rdt,
                 "Output dtype is %s, expected %s" % (y1.dtype, self.rdt))
         self.failUnless(y2.dtype == self.rdt,
                 "Output dtype is %s, expected %s" % (y2.dtype, self.rdt))
         assert_array_almost_equal (y1, x, decimal=self.ndec)
         assert_array_almost_equal (y2, x, decimal=self.ndec)

    def test_size_accuracy(self):
        # Sanity check for the accuracy for prime and non-prime sized inputs
        if self.rdt == np.float32:
            rtol = 1e-5
        elif self.rdt == np.float64:
            rtol = 1e-10

        for size in LARGE_COMPOSITE_SIZES + LARGE_PRIME_SIZES:
            np.random.seed(1234)
            x = np.random.rand(size).astype(self.rdt)
            y = irfft(rfft(x))
            _assert_close_in_norm(x, y, rtol, size, self.rdt)
            y = rfft(irfft(x))
            _assert_close_in_norm(x, y, rtol, size, self.rdt)

	def sineFit(self,xReal,yReal):
		N=len(xReal)
		OFFSET = (yReal.max()+yReal.min())/2.
		yhat = fftpack.rfft(yReal-OFFSET)
		idx = (yhat**2).argmax()
		freqs = fftpack.rfftfreq(N, d = (xReal[1]-xReal[0])/(2*np.pi))
		frequency = freqs[idx]/(2*np.pi)  #Convert angular velocity to freq

		amplitude = (yReal.max()-yReal.min())/2.0
		phase=0#.5*np.pi*((yReal[0]-offset)/amplitude)
		guess = [amplitude, frequency, phase,0]
		try:
			(amplitude, frequency, phase,offset), pcov = optimize.curve_fit(self.sineFunc, xReal, yReal-OFFSET, guess)
			offset+=OFFSET
			ph = ((phase)*180/(np.pi))
			if(frequency<0):
				#print 'negative frq'
				return False

			if(amplitude<0):
				ph-=180

			if(ph<0):ph = (ph+720)%360
			freq=1e6*abs(frequency)
			amp=abs(amplitude)
			pcov[0]*=1e6
			#print pcov
			if(abs(pcov[-1][0])>1e-6):
				False
			return [amp, freq, offset,ph]
		except:
			return False

def do_gen_random(peakAmpl, durationInMSec, samplingRate, fHigh, stereo=True):
    samples = durationInMSec * samplingRate / 1000
    result = np.zeros(samples * 2 if stereo else samples, dtype=np.int16)
    randomSignal = np.random.normal(scale = peakAmpl * 2 / 3, size=samples)
    fftData = fft.rfft(randomSignal)
    freqSamples = samples/2
    iHigh = freqSamples * fHigh * 2 / samplingRate + 1
    #print len(randomSignal), len(fftData), fLow, fHigh, iHigh
    if iHigh > freqSamples - 1:
        iHigh = freqSamples - 1
    fftData[0] = 0 # DC
    for i in range(iHigh, freqSamples - 1):
        fftData[ 2 * i + 1 ] = 0
        fftData[ 2 * i + 2 ] = 0
    if (samples - 2 *freqSamples) != 0:
        fftData[samples - 1] = 0

    filteredData = fft.irfft(fftData)
    #freq = np.linspace(0.0, samplingRate, num=len(fftData), endpoint=False)
    #plt.plot(freq, abs(fft.fft(filteredData)))
    #plt.plot(filteredData)
    #plt.show()
    if stereo:
        for i in range(len(filteredData)):
            result[2 * i] = filteredData[i]
            result[2 * i + 1] = filteredData[i]
    else:
        for i in range(len(filteredData)):
            result[i] = filteredData[i]
    return result

def show_magnitued(file_name):
    # read audio samples
    input_data = read(file_name)
    audio = input_data[1]
    print(audio)
    # apply a Hanning window
    window = hann(1024)
    audio = audio[0:1024] * window
    # fft
    mags = abs(rfft(audio))
    # convert to dB
    mags = 20 * scipy.log10(mags)
    # normalise to 0 dB max
    # mags -= max(mags)
    file = open('tmp.txt', 'w')
    for i in mags:
        file.write(str(i) + '\n')
    file.close()
    # plot
    plt.plot(mags)
    # label the axes
    plt.ylabel("Magnitude (dB)")
    plt.xlabel("Frequency Bin")
    # set the title
    plt.title(file_name + " Spectrum")
    plt.show()

def bandpass(x, sampling_rate, f_min, f_max, verbose=0):
	"""
	xf = bandpass(x, sampling_rate, f_min, f_max)

	Description
	--------------

	Phasen-treue mit der rueckwaerts-vorwaerts methode!
	Function bandpass-filters a signal without roleoff.  The cutoff frequencies,
	f_min and f_max, are sharp.

	Arguements
	--------------
		x: 			input timeseries
		sampling_rate: 			equidistant sampling with sampling frequency sampling_rate
		f_min, f_max:			filter constants for lower and higher frequency
	
	Returns
	--------------
		xf:		the filtered signal
	"""
	x, N = np.asarray(x, dtype=float), len(x)
	t = np.arange(N)/np.float(sampling_rate)
	xn = detrend_linear(x)
	del t

	yn = np.concatenate((xn[::-1], xn))		# backwards forwards array
	f = np.float(sampling_rate)*np.asarray(np.arange(2*N)/2, dtype=int)/float(2*N)
	s = rfft(yn)*(f>f_min)*(f<f_max)			# filtering

	yf = irfft(s)					# backtransformation
	xf = (yf[:N][::-1]+yf[N:])/2.			# phase average

	return xf

def rfft_freq(data, window_func=signal.hanning):
    w = window_func(data.size)
    sig_fft = fftpack.rfft(data * w)
    freq = fftpack.rfftfreq(sig_fft.size, d=SAMPLING_INTERVAL)
    freq = freq[range(data.size / 2)]
    sig_fft = sig_fft[range(data.size / 2)]
    return sig_fft, freq

 def fft(self):
     t_plot = np.arange(0,self.T,self.dt)
     print len(t_plot)
     amplitude_record = []
     for i in range(int(self.T/self.dt)):
         amplitude_record.append(self.string[i][20])
     print len(amplitude_record)
     '''
     plt.subplot(121)
     plt.title('String signal versus time')
     plt.ylabel('Signal (arbitrary units)')
     plt.xlabel('Time (s)')
     plt.plot(t_plot,amplitude_record)
     '''
     freq = fftfreq(len(amplitude_record), d=self.dt)
     freq = np.array(abs(freq))
     f_signal = rfft(amplitude_record)
     f_signal = np.array(f_signal**2)
     #plt.subplot(122)
     plt.title('Power spectra')
     plt.ylabel('Power (arbitrary units)')
     plt.xlabel('Frequency (Hz)')
     plt.xlim(2000,8000)
     plt.plot(freq,f_signal,label = 'epsilon = '+str(self.e))
     return 0

def select_events(nevents,nfeatures):
    global groups
    fftbins = 8192
    featurewidth = 16
    print "Selecting %d random spectral features.." % nfeatures
    feature_bins = np.random.randint(featurewidth/2,(fftbins/8),nfeatures)
    print "Selecting %d random audio events.." % nevents
    events = np.random.randint(0,len(faudio)-grain_mid,nevents)
    # Initialise features array with the first variable as index
    features = np.zeros((nfeatures+1,nevents))
    features[0] = np.arange(0,nevents)
    print "Computing audio event spectrograms.."
    # For each event..
    for i in range(0,nevents):
        # Calculate spectrogram for the event
        _fftevent = faudio[events[i]:min(events[i]+grain_mid,len(faudio))]*sig.hann(grain_mid)
        mags = abs(rfft(_fftevent,fftbins))
        mags = 20*log10(mags) # dB
        mags -= max(mags) # normalise to 0dB max
        # Calculate each feature for this event
        for j in range(0,nfeatures):
            features[j+1][i] = abs(np.mean(abs(mags[(feature_bins[j]-featurewidth/2):(feature_bins[j]+featurewidth/2)])))
    print "Clustering events with K-Means algorithm.."
    groups = kmeans(np.transpose(features[1:,:]),tracks,minit='points',iter=30)[1]
    return [events,groups]

def compute_fft(fs, ir):

    import scipy.signal, numpy
    from scipy import fftpack

    # creating asymmetric bartlett window for spectral analysis
    window_bart = scipy.signal.bartlett(len(ir), sym=False)

    # windowing the impulse response
    ir_wind = ir * window_bart

    # computing fft
    sig_fft = fftpack.rfft(ir_wind)

    # setting length of fft
    n = sig_fft.size
    timestep = 1 / float(fs)

    # generating frequencies according to fft points
    freq = fftpack.rfftfreq(n, d=timestep)

    # normalizing fft
    sys_fft = abs(sig_fft) / n

    # TODO FFT computing
    return sys_fft, freq

	def calculate_attributes(self):
		source = self.source
		freq = self.frequency
		sampling_rate = float(source.sampling_rate)
		fft_sampling_rate = sampling_rate/float(source.fft_step_size)
		window_length = float(source.fft_window_size)/sampling_rate
		# FIXME real time should be passed in as an extra field
		self.start = window_length
		if self.first_frame > 0:
			self.start += (self.first_frame - 1)/fft_sampling_rate
		self.length = window_length
		if len(freq) > 1:
#			if 'length' not in self.__dict__:
#				print self.__dict__
#			print self.__dict__['length']
			self.length += (len(freq) - 1)/fft_sampling_rate
#		print self.length
		self.end = self.start + self.length
		for k in ('frequency','amplitude'):
			a = getattr(self,k)
			setattr(self, k+'_min', min(a))
			setattr(self, k+'_max', max(a))
			setattr(self, k+'_mean', sum(a)/len(a))
		freq_window = 2 # seconds
		freq_fft_size = 128
		resampled_freq = resample(freq, freq_window*freq_fft_size/fft_sampling_rate, 'sinc_fastest') # FIXME truncate array?
		self.freq_fft = abs(rfft(resampled_freq,n=freq_fft_size,overwrite_x=True))[1:]

def fft_filter(x, fs, band=(9, 14)):
    w = fftfreq(x.shape[0], d=1. / fs * 2)
    f_signal = rfft(x)
    cut_f_signal = f_signal.copy()
    cut_f_signal[(w < band[0]) | (w > band[1])] = 0
    cut_signal = irfft(cut_f_signal)
    return cut_signal

    def test_size_accuracy(self):
        # Sanity check for the accuracy for prime and non-prime sized inputs
        if self.rdt == np.float32:
            rtol = 1e-5
        elif self.rdt == np.float64:
            rtol = 1e-10

        for size in LARGE_COMPOSITE_SIZES + LARGE_PRIME_SIZES:
            np.random.seed(1234)
            x = np.random.rand(size).astype(self.rdt)
            y = irfft(rfft(x))
            self.failUnless(np.linalg.norm(x - y) < rtol*np.linalg.norm(x),
                            (size, self.rdt))
            y = rfft(irfft(x))
            self.failUnless(np.linalg.norm(x - y) < rtol*np.linalg.norm(x),
                            (size, self.rdt))

def fftresample(S, npoints, reflect=False, axis=0):
    """
    Resample a signal using discrete fourier transform. The signal
    is transformed in the fourier domain and then padded or truncated
    to the correct sampling frequency.  This should be equivalent to
    a sinc resampling.
    """
    from scipy.fftpack import rfft, irfft
    from dlab.datautils import flipaxis

    # this may be considerably faster if we do the memory operations in C
    # reflect at the boundaries
    if reflect:
        S = nx.concatenate([flipaxis(S,axis), S, flipaxis(S,axis)],
                           axis=axis)
        npoints *= 3

    newshape = list(S.shape)
    newshape[axis] = int(npoints)

    Sf = rfft(S, axis=axis)
    Sr = (1. * npoints / S.shape[axis]) * irfft(Sf, npoints, axis=axis, overwrite_x=1)
    if reflect:
        return nx.split(Sr,3)[1]
    else:
        return Sr

def fft_filter(x, fs, band=(9, 14)):
    w = fftpack.rfftfreq(x.shape[0], d=1. / fs)
    f_signal = fftpack.rfft(x, axis=0)
    cut_f_signal = f_signal.copy()
    cut_f_signal[(w < band[0]) | (w > band[1])] = 0
    cut_signal = fftpack.irfft(cut_f_signal, axis=0)
    return cut_signal

def get_power2(x, fs, band, n_sec=5):
    n_steps = int(n_sec * fs)
    w = fftpack.fftfreq(n_steps, d=1. / fs * 2)
    print(len(range(0, x.shape[0] - n_steps, n_steps)))
    pows = [2*np.sum(fftpack.rfft(x[k:k+n_steps])[(w > band[0]) & (w < band[1])]**2)/n_steps
            for k in range(0, x.shape[0] - n_steps, n_steps)]
    return np.array(pows)