def sanity(wavobj):
    transformed = wavobj['transformed']
    transformed_raw = wavobj['raw']
    rate = wavobj['rate']
    output = irfft(transformed)
    san_output = irfft(transformed)* transformed_raw.max(axis=0)

    write('sanity.wav', rate, np.array(san_output, dtype=np.int16))
    print("Sanity.wav written")

 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_definition(self):
     x = [1,2,3,4,1,2,3,4]
     x1 = [1,2+3j,4+1j,2+3j,4,2-3j,4-1j,2-3j]
     y = irfft(x)
     y1 = direct_irdft(x)
     assert_array_almost_equal(y,y1)
     assert_array_almost_equal(y,ifft(x1))
     x = [1,2,3,4,1,2,3,4,5]
     x1 = [1,2+3j,4+1j,2+3j,4+5j,4-5j,2-3j,4-1j,2-3j]
     y = irfft(x)
     y1 = direct_irdft(x)
     assert_array_almost_equal(y,y1)
     assert_array_almost_equal(y,ifft(x1))

    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 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 savefft(wavfile, wavobj, filtered):
    transformed = wavobj['transformed']
    transformed_raw = wavobj['raw']
    rate = wavobj['rate']
    data = filtered * transformed_raw.max(axis=0)
    output = irfft(data)
    write(wavfile, rate, np.array(output, dtype=np.int16))

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 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 = 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(x, xr):
     y = irfft(np.array(x, dtype=self.rdt))
     y1 = direct_irdft(x)
     self.assertTrue(y.dtype == self.rdt,
             "Output dtype is %s, expected %s" % (y.dtype, self.rdt))
     assert_array_almost_equal(y,y1, decimal=self.ndec)
     assert_array_almost_equal(y,ifft(xr), decimal=self.ndec)

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 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 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 testFFTComp(self):
     data = np.array([1, 0, 0, 1, 1, 0, 0, 1], dtype=np.float32)
     self.stft.performStft(data)
     dfts = self.stft.getDFTs()
     transformed = mat.to_numpy_format(dfts)
     ifftout = fftp.irfft(transformed[:, 0] / 2)
     print ifftout
     self.assertListFloatEqual(data, ifftout)

    def applyParameter(self, data):
        
        #low pass filter
        fftdata = rfft(data)
        fftdata[int(float(self._variables["FREQUENCY"])*np.pi*2):] = 0

        ##removes negative values
        filtered_data = irfft(fftdata).clip(0)

        return filtered_data

    def irfft(self, X):
        """
        Apply real FFT to X

        Parameters
        ----------
        X     : list of lists or ndArrays
        """
        X = np.array([irfft(x) for x in X])
        return X

def bandpass(data_list, min_hz, max_hz):

    fft_list = rfft(data_list)

    # Filter
    for i in range(len(fft_list)):
        if not (min_hz < i/2+1 < max_hz): fft_list[i] = 0

    result_vals = irfft(fft_list)

    return result_vals

def remove_freq_range(data_list, min_hz, max_hz):

    fft_list = rfft(data_list)

    # Filter
    for i in range(len(fft_list)):
        if (min_hz < i / 2 + 1 < max_hz): fft_list[i] = 0

    result_vals = irfft(fft_list)

    return result_vals

def ourFFT2(signal,ffTtime):
    W = fftfreq(len(signal), d=ffTtime[1]-ffTtime[0])
    f_signal = rfft(signal)

    # If our original signal time was in seconds, this is now in Hz    
    cut_f_signal = f_signal.copy()
    #cut_f_signal[(W>0.004)] = 0
    cut_f_signal[(W<=0)] = 0
    cut_f_signal[(W>0.001)] = 0

    cut_signal = irfft(cut_f_signal)
    return W, ffTtime, cut_signal, cut_f_signal, f_signal