def problem4():
	# read in tada.wav
	rate, tada = wavfile.read('tada.wav')
	
	# upon inspection, we find that tada.wav is a stereo audio file. 
	# we create stereo white noise that lasts 10 seconds
	L_white = sp.int16(sp.random.randint(-32767,32767,rate*10))
	R_white = sp.int16(sp.random.randint(-32767,32767,rate*10))
	white = sp.zeros((len(L_white),2))
	white[:,0] = L_white
	white[:,1] = R_white
	
	# pad tada signal with zeros
	padded_tada = sp.zeros_like(white)
	padded_tada[:len(tada)] = tada
	ptada = padded_tada
	
	# fourier transforms
	ftada = sp.fft(ptada,axis=0)
	fwhite = sp.fft(white,axis=0)
	
	# inverse transform of convolution
	out = sp.ifft((ftada*fwhite),axis=0)
	
	# prepping output and writing file
	out = sp.real(out)
	scaled = sp.int16(out / sp.absolute(out).max() * 32767)
	wavfile.write('my_tada_conv.wav',rate,scaled)

def get_envelope(R,dim=1):
    """
    Returns the complex version of the input signal R.
    @param R: The input data matrix.
    @param dim: The dimension along which the envelope is to be taken. default: dim=1
    """
    if dim==0:
        R=R.T
    if len(R.shape)==1:
        freqs=scipy.fft(R)
        length=len(R)/2
        freqs[length:]=0
        freqs[1:length]=2*freqs[1:length]
        ## freqs[1:length]=freqs[1:length]
        
        env=scipy.ifft(freqs)
    else:
        freqs=scipy.fft(R)
        length=R.shape[dim]/2
        #Something is fishy here:
        freqs[:,length:]=0
        freqs[:,1:length]=2*freqs[0,1:length]
        ## freqs[:,1:length]=freqs[0,1:length]
        
        env=scipy.ifft(freqs)
        if dim==0:
            return env.T

    return env

def istft(X, fs, T, hop):
    x = scipy.zeros(T * fs)
    framesamp = X.shape[1]
    hopsamp = int(hop * fs)
    for n, i in enumerate(range(0, len(x)-framesamp, hopsamp)):
        x[i:i+framesamp] += scipy.real(scipy.ifft(X[n]))
    return x

def mmse_stsa(infile, outfile, noise_sum):
    signal, params = read_signal(infile, WINSIZE)
    nf = len(signal)/(WINSIZE/2) - 1
    sig_out=sp.zeros(len(signal),sp.float32)

    G = sp.ones(WINSIZE)
    prevGamma = G
    alpha = 0.98
    window = sp.hanning(WINSIZE)
    gamma15=spc.gamma(1.5)
    lambdaD = noise_sum / 5.0
    percentage = 0
    for no in xrange(nf):
        p = int(math.floor(1. * no / nf * 100))
        if (p > percentage):
            percentage = p
            print "{}%".format(p),

        y = get_frame(signal, WINSIZE, no)
        Y = sp.fft(y*window)
        Yr = sp.absolute(Y)
        Yp = sp.angle(Y)
        gamma = Yr**2/lambdaD
        xi = alpha * G**2 * prevGamma + (1-alpha)*sp.maximum(gamma-1, 0)
        prevGamma = gamma
        nu = gamma * xi / (1+xi)
        G = (gamma15 * sp.sqrt(nu) / gamma ) * sp.exp(-nu/2) * ((1+nu)*spc.i0(nu/2)+nu*spc.i1(nu/2))
        idx = sp.isnan(G) + sp.isinf(G)
        G[idx] = xi[idx] / (xi[idx] + 1)
        Yr = G * Yr
        Y = Yr * sp.exp(Yp*1j)
        y_o = sp.real(sp.ifft(Y))
        add_signal(sig_out, y_o, WINSIZE, no)
    
    write_signal(outfile, params, sig_out)

def ffacorr(a):
 """Returns the autocorrelation of a. Expects raw data"""
 z=np.zeros(2*len(a))
 z[:len(a)]=a
 fft=sc.fft(z)
 out=sc.ifft(fft*sc.conj(fft))
 return (out[:len(out)/2])

def cwt_freq(data, wavelet, widths, dt, axis):
    # compute in frequency
    # next highest power of two for padding
    N = data.shape[axis]
    pN = int(2 ** np.ceil(np.log2(N)))
    # N.B. padding in fft adds zeros to the *end* of the array,
    # not equally either end.
    fft_data = scipy.fft(data, n=pN, axis=axis)
    # frequencies
    w_k = np.fft.fftfreq(pN, d=dt) * 2 * np.pi

    # sample wavelet and normalise
    norm = (2 * np.pi * widths / dt) ** .5
    wavelet_data = norm[:, None] * wavelet(w_k, widths[:, None])

    # Convert negative axis. Add one to account for
    # inclusion of widths axis above.
    axis = (axis % data.ndim) + 1

    # perform the convolution in frequency space
    slices = [slice(None)] + [None for _ in data.shape]
    slices[axis] = slice(None)

    out = scipy.ifft(fft_data[None] * wavelet_data.conj()[slices],
                     n=pN, axis=axis)

    # remove zero padding
    slices = [slice(None) for _ in out.shape]
    slices[axis] = slice(None, N)

    if data.ndim == 1:
        return out[slices].squeeze()
    else:
        return out[slices]

    def propagate(self):
        r"""
        Given the wavefunction values :math:`\Psi` at time :math:`t`, calculate new
        values at time :math:`t + \tau`. We perform exactly one timestep :math:`\tau` here.
        """
        # How many states we have
        nst = self.Psi.get_number_components()

        # Read values out of current WaveFunction state
        vals = self.Psi.get_values()

        # Do the propagation
        tmp = [ zeros(vals[0].shape, dtype=complexfloating) for item in vals ]
        for row in xrange(0, nst):
            for col in xrange(0, nst):
                tmp[row] = tmp[row] + self.VE[row*nst+col] * vals[col]

        tmp = tuple([ fft(item) for item in tmp ])

        tmp = tuple([ self.TE * item for item in tmp ])

        tmp = tuple([ ifft(item) for item in tmp ])

        values = [ zeros(tmp[0].shape, dtype=complexfloating) for item in tmp ]
        for row in xrange(0, nst):
            for col in xrange(0, nst):
                values[row] = values[row] + self.VE[row*nst+col] * tmp[col]

        # Write values back to WaveFunction object
        self.Psi.set_values(values)

def istft(spectrogram, w_size, step):
    
    if spectrogram.shape[0] != w_size:
        print ("Mismatch w_size and spectrogram")

    eps = np.finfo(float).eps
    window = np.hanning(w_size)
    # spectrogram.shape = w_size , bins
    spectr_len = spectrogram.shape[1]
    reconst_len = w_size + (spectr_len - 1) * step

    reconst_x = np.zeros(reconst_len, dtype=float)
    windowsum = np.zeros(reconst_len, dtype=float)
    windowsq = window * window

    # Overlap add
    for i in range(0, spectr_len):
        s = i * step
        e = i * step + w_size
        r = ifft(spectrogram[:, i]).real
        # r = abs(ifft(spectrogram[:, i]))

        reconst_x[s:e] += r * window
        windowsum[s:e] += windowsq

    # Normalize by window

    # for i in range(0,reconst_len)
    #     if windowsum[i] > eps
    #         reconst_x[i] /= windowsum[i]
    pos = (windowsum != 0)
    reconst_x[pos] /= windowsum[pos]
    return reconst_x.astype("int16")

def pCodePhaseSearch(sample, LOFreq, PRNSpectrumConjugate):

    """ Search in a given LO freq space all the code phases
        at once

    """

    I, Q, _ = generateIQ(sample.size, LOFreq)

    # mix is down with the LO
    sampledMixedI = sample * I
    sampledMixedQ = sample * Q

    # munge them into a single array of complex numbers for the fft
    combinedMixed = sampledMixedI + 1j*sampledMixedQ

    # do the fft
    signalSpectrum = scipy.fft(combinedMixed)

    # circulator correlation in da frequency space
    correlatedSpectrum = signalSpectrum * PRNSpectrumConjugate

    # and back to time domain
    timeDomainReconstructed = np.abs(scipy.ifft(correlatedSpectrum))**2

    return timeDomainReconstructed

def bandpass_ifft(X, Low_cutoff, High_cutoff, F_sample, M=None):
    """Bandpass filtering on a real signal using inverse FFT
    
    Inputs
    =======
    
    X: 1-D numpy array of floats, the real time domain signal (time series) to be filtered
    Low_cutoff: float, frequency components below this frequency will not pass the filter (physical frequency in unit of Hz)
    High_cutoff: float, frequency components above this frequency will not pass the filter (physical frequency in unit of Hz)
    F_sample: float, the sampling frequency of the signal (physical frequency in unit of Hz)    
    
    Notes
    n=====
    1. The input signal must be real, not imaginary nor complex
    2. The Filtered_signal will have only half of original amplitude. Use abs() to restore. 
    3. In Numpy/Scipy, the frequencies goes from 0 to F_sample/2 and then from negative F_sample to 0. 
    
    """        
    import numpy, scipy
   
    if M == None: # if the number of points for FFT is not specified
        M = X.size # let M be the length of the time series
    
    #Spectrum = scipy.fft(X, n=M) 
    Spectrum = X

    [Low_cutoff, High_cutoff, F_sample] = map(float, [Low_cutoff, High_cutoff, F_sample])
    
    #Convert cutoff frequencies into points on spectrum
    #[Low_point, High_point] = map(lambda F: F/F_sample * M /2, [Low_cutoff, High_cutoff])# the division by 2 is because the spectrum is symmetric 
    [Low_point, High_point] = map(lambda F: F/F_sample * M /1, [Low_cutoff, High_cutoff])# the division by 2 is because the spectrum is symmetric 

    Filtered_spectrum = [Spectrum[i] if i >= Low_point and i <= High_point else 0.0 for i in xrange(M)] # Filtering
    Filtered_signal = scipy.ifft(Filtered_spectrum, n=M)  # Construct filtered signal 
    return Spectrum, Filtered_spectrum, Filtered_signal, Low_point, High_point

def istft(X, chunk_size, hop, w=None):
    """
    Naively inverts the short time fourier transform using an overlap and add
    method. The overlap is defined by hop

    Args:
      X: STFT windows to invert, overlap and add. 
      chunk_size: size of analysis window.
      hop: hop distance between analysis windows
      w: windowing function to apply. Must be of length chunk_size

    Returns:
      ISTFT of X using an overlap and add method. Windowing used to smooth.

    Raises:
      ValueError if window w is not of size chunk_size
    """

    if not w:
        w = sp.hanning(chunk_size)
    else:
        if len(w) != chunk_size:
            raise ValueError("window w is not of the correct length {0}.".format(chunk_size))

    x = sp.zeros(len(X) * (hop))
    i_p = 0
    for n, i in enumerate(range(0, len(x)-chunk_size, hop)):
        x[i:i+chunk_size] += w*sp.real(sp.ifft(X[n]))
    return x

def filter_wav(input_file, output_file = "filtered.wav"):
    """ filters input_file and save resulting data to output_file """
    raw = invert(readwave(input_file))  
    fftr = numpy.fft.fft(raw)  
    fft = fftr[:]  
    fftx= numpy.fft.fftfreq(len(raw), d=(1.0/(rate)))
    
    lowfreq = 58
    highfreq = 62
    lowspectrum = [(lowfreq+(item*60)) for item in range(5)]
    highspectrum = [(highfreq+(item*60)) for item in range(5)]

    fft=bandStop(fft,fftx,0,20)
    fft=bandStop(fft,fftx,-20,0)
    for i in range(5):
         fft = bandStop(fft,fftx,lowspectrum[i],highspectrum[i])
    fft=bandStop(fft,fftx,300,max(fftx))

    fix = scipy.ifft(fft)
    smoothed = Hanning(fix)
    gain = 1.0/max(numpy.absolute(smoothed))
    nsamples = len(smoothed)
    asamples = numpy.zeros(nsamples,dtype=numpy.int16)
    numpy.multiply(smoothed,32767.0 * gain,asamples)
    wavfile.write(output_file,rate,asamples)

def delayedsignalF(x,t0_pts):
#==============================================================
    """
     Delay a signal with a non integer value
     (computation in frequency domain)
     
     Synopsis:
              y=delayedsignalF(x,t0_pts)
     
     Inputs: x vector of length N
             t0_pts is a REAL delay
             expressed wrt the sampling time Ts=1:
               t0_pts = 1 corresponds to one time dot
               t0_pts may be positive, negative, non integer
               t0_pts>0: shift to the right
               t0_pts<0: shift to the left
     Rk: the length of FFT is 2^(nextpow2(N)+1
    """
    #
    # M. Charbit, Jan. 2010
    #==============================================================
    N         = len(x)
    p         = ceil(log2(N))+1;
    Lfft      = int(2.0**p);
    Lffts2    = Lfft/2;
    fftx      = fft(x, Lfft);
    ind       = concatenate((range(Lffts2+1), 
            range(Lffts2+1-Lfft,0)),axis=0)
    fftdelay  = exp(-2j*pi*t0_pts*ind/Lfft);
    fftdelay[Lffts2] = real(fftdelay[Lffts2]);
    ifftdelay        = ifft(fftx*fftdelay);
    y                = ifftdelay[range(N)];
    if isreal(any(x)):
        y=real(y)
    return y

def IMRpeakAmp(m1,m2,spin1z,spin2z,d):

    """
	IMRpeakAmp finds the peak amplitude of the waveform for a given source parameters and the source distance.
	usage: IMRpeakAmp(m1,m2,spin1z,spin2z,distance)
	e.g.
	spawaveApp.IMRpeakAmp(30,40,0.45,0.5,100)

	"""

    chi = spawaveform.computechi(m1, m2, spin1z, spin2z)   
    imrfFinal = spawaveform.imrffinal(m1, m2, chi, 'fcut')
    fLower = 10.0
    order = 7
    dur = 2**numpy.ceil(numpy.log2(spawaveform.chirptime(m1,m2,order,fLower)))
    sr = 2**numpy.ceil(numpy.log2(imrfFinal*2))
    deltaF = 1.0 / dur
    deltaT = 1.0 / sr
    s = numpy.empty(sr * dur, 'complex128')	
    spawaveform.imrwaveform(m1, m2, deltaF, fLower, s, spin1z, spin2z)
    s = scipy.ifft(s)
    #s = numpy.abs(s)
    s = numpy.real(s)
    max = numpy.max(s)/d
    return max

def HT(y):
	"""
	HT(y)

	Computes the Hilbert transform of the array y using the 
	frequency-domain approach, as in
	http://library2.usask.ca/theses/available/etd-09272006-135609/unrestricted/XianglingWang.pdf, pages 36-37

	For a good transform, the following requirements must be satisfied:
	1. The numbers of zeros and local extrema in y must differ by at most one.
	2. y must be symmetric about zero (i.e. the mean value of the envelopes 
		defined the by the local maxima and minima of y must be zero).
	"""

	from scipy import fft, ifft
	from numpy import zeros_like

	#Take the Fourier transform of the function
	ffty = fft(y)

	#Write as the FFT of the Hilbert transform
	Y = zeros_like(ffty)
	N = len(ffty)
	for i in range(-N/2+1,N/2+1):
		Y[i] = -1j*sgn(i)*ffty[i]

	#Take the inverse Fourier transform
	HT = ifft(Y)

	return HT

def start():
    
        global signal
        global heartrate
        global spo2
        
        easypulse = EasyPulse()
        beats = easypulse.readPulse()
        heartrate = easypulse.computeHeartrate(beats)
        spo2 = spo2()
        
        #print ("Number of beats: " + str(len(beats)))
        #print( "Heartrate: " + str(heartrate))
        
        signal=[]
        
        for i in range(1,100):
            
            reading = easypulse.readadc(2, 11, 10, 9, 8)
            signal.append(reading*3.3/1023)
            
        fft=scipy.fft(signal)
        bp=fft[:]
        
        for i in range(len(bp)): 
         if i>=10:bp[i]=0  

        ibp=scipy.ifft(bp)

def down_sample(filename, new_rate, outputfile=None):
    """
    Create a down-sampled copy of the provided .wav file.  Unless overridden, the output
        file will be of the form "down_<orginalname>.wav"
        
    Parameters
    ----------
    filename : string
        input .wav file
    new_rate : int
        sample rate of output file
    outputfile : string
        name of output file
    """

    if outputfile is None:
        outputfile = "down_" + filename

    old_rate, in_sig = wavfile.read(filename)
    in_sig = sp.float32(in_sig)
    fin = sp.fft(in_sig)
    nsiz = sp.floor(in_sig.size * new_rate / old_rate)
    nsizh = sp.floor(nsiz / 2)
    fout = sp.zeros(nsiz)
    fout = fout + 0j
    fout[0:nsizh] = fin[0:nsizh]
    fout[nsiz - nsizh + 1 :] = sp.conj(sp.flipud(fout[1:nsizh]))
    out = sp.ifft(fout)
    out = sp.real(out)  # Take the real component of the signal
    out = sp.int16(out / sp.absolute(out).max() * 32767)
    wavfile.write(outputfile, new_rate, out)

def get_noise_freq_domain_1NSD( P , df , inittime , parityN , seed ) :
    """
    returns the noise time-sereis given the power spectral density
    INPUT:
    P --- power spectral density. numpy array
    df --- frequency resolution
    inittime --- initial time of the noise time-series
    parityN --- is the length of the time-series 'Odd' or 'Even'
    seed --- seed for the noise
    OUPUT:
    t --- time [s]
    n --- noise time-series
    """
    Nf = P.shape[0]
    
    if parityN == 'Odd' :
        N = 2 * Nf + 1
    elif parityN == 'Even' :
        N = 2 * ( Nf + 1 ) 
    else :
        raise InputError , "parityN must be either 'Odd' or 'Even'!"
    stime = 1 / ( N*df )
    t = inittime + stime * np.arange( N )
    
    np.random.seed( seed )
    z = ( np.random.standard_normal( P.shape ) + 1j * np.random.standard_normal( P.shape ) ) / np.sqrt( 2 )

    ntilde_fplus = np.sqrt( N / ( 2*stime ) * P ) * z 

    if N % 2 == 0 :
        ntilde = np.array( [0] + list( ntilde_fplus ) + [0] + list( np.flipud(np.conj(ntilde_fplus)) ) )
    else :
        ntilde = np.array( [0] + list( ntilde_fplus ) + list( np.flipud(np.conj(ntilde_fplus)) ) )
    n = np.real( sp.ifft( ntilde ) )
    return t , n

def violent_multi_band_pass(signal, mask):
    """ A filter that lets you define a bunch of bands:
    everything outside of these frequency ranges gets mercylessly filtered"""
    assert(len(signal) == len(mask))
    f = fft(signal)
    filtered = [mask[i] * f[i] for i in range(len(signal))]
    return scipy.ifft(filtered).real.tolist()

 def ApplyFiltersWall(self, data, lowFreq, highFreq, SampleFreq=1250000.0):
     '''Applies FFT software filters.
     In -> Data, lowFreq, highFreq
     Out-> filteredSignal
     
     Warning: Appling a blank wall filter is the same as
     convolving with sinc function!!
     THIS PRODUCES RINGING!!
     '''
     #Perform FFT
     fftsigs=numpy.fft.fft(data)
     fftfreqs = numpy.fft.fftfreq(len(data),1.0/SampleFreq)
     
     
     #Apply blank wall filter by cutting off the unwanted frequencies
     if lowFreq is not None:
         lowfreq_pos =  [i for i,x in enumerate(fftfreqs) if lowFreq-10<=x<=lowFreq+10]
         for i in range(lowfreq_pos[0]):
                 fftsigs[i]=0
         #for i in range(len(fftsigs)):
         #    if abs(fftfreqs[i])<=lowFreq:
         #        fftsigs[i]=0
                 
     if highFreq is not None:
         for i in range(len(fftsigs)):
             if fftfreqs[i]>=highFreq:
                 fftsigs[i]=0
     #Regenerate Signal
     return scipy.ifft(fftsigs, n=len(fftsigs)/2)