def getMSD(d1, d2):
    diff = 0
    d11 = []
    d22 = []
    
    index = 0
    for index in xrange(9):
        chunk = d1[index * 22050: index * 22050 + 22050]
        x = cp.fft(chunk)
        maximum1 = max(x)
        d11.append(maximum1)
         
    index = 0
    for index in xrange(9):
        chunk = d2[index * 22050: index * 22050 + 22050]
        x = cp.fft(chunk)
        maximum1 = max(x)
        d22.append(maximum1)
    
    currentMax1 = max(d11)
    currentMax2 = max(d22)
    
    for x in xrange(len(d11)):
#         diff += abs(d11[x]/currentMax1 - d22[x]/currentMax2) * abs(d11[x]/currentMax1 - d22[x]/currentMax2)
        diff += abs(d11[x] - d22[x]) * abs(d11[x] - d22[x])

    diff /= len(d11)
    
    return diff

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 itakura_saito_spectrum_distance(s,shat,winfunc):
    size = min(len(s),len(shat))
    window = winfunc(size)
    s = s[0:size]
    shat = shat[0:size]
    s_amp = sp.absolute(sp.fft(s*window))
    shat_amp = sp.absolute(sp.fft(shat*window))
    return sp.mean(sp.log(s_amp / shat_amp) + (shat_amp/s_amp) - 1.0)

def log_spectrum_distance(s,shat,winfunc):
    size = min(len(s),len(shat))
    window = winfunc(size)
    s = s[0:size]
    shat = shat[0:size]
    s_amp = sp.absolute(sp.fft(s*window))
    shat_amp = sp.absolute(sp.fft(shat*window))
    return sp.sqrt(sp.mean((sp.log10(s_amp / shat_amp)*10.0)**2.0))

def prob5():
	rate, sig = wavfile.read('tada.wav')
	sig = sp.float32(sig)
	noise = sp.float32(sp.random.randint(-32767,32767,sig.shape))
	out = sp.ifft(sp.fft(sig)*sp.fft(noise))
	out = sp.real(out)
	out = sp.int16(out/sp.absolute(out).max() * 32767)
	wavfile.write('white-conv.wav',rate,out)

    def fourierTransform(self, fromPos, toPos, only = []):
        self.checkToPos(toPos)
        if len(only) > 0:
            self.allFid[toPos] = np.array([fftshift(fft(self.allFid[fromPos][fidIndex])) for fidIndex in only])
        else:
            self.allFid[toPos] = np.array([fftshift(fft(fid)) for fid in self.allFid[fromPos]])

        self.frequency = np.linspace(-self.sweepWidthTD2/2,self.sweepWidthTD2/2,len(self.allFid[fromPos][0]))

def fftconv(x, y):
    """ Convolution of x and y using the FFT convolution theorem. """
    n = np.int(np.round(2 ** np.ceil(np.log2(len(x))))) + 1
    X, Y, x_y = fft(x, n), fft(y, n), []
    for i in range(n):
        x_y.append(X[i] * Y[i])

    # Returns the inverse Fourier transform with padding correction
    return fft.ifft(x_y)[4:len(x)+4]

def algoChannelSelection(left, right):

    ''' Algorithm which automatically selects the channel with dominant vocals from a stereo flamenco recording
    based on spectral band energies as described in section 2-A-I of

    Kroher, N. & Gomez, E. (2016). Automatic Transcription of Flamenco Singing from Polyphonic Music Recordings.
    ACM / IEEE Transactions on Audio, Speech and Language Processing, 24(5), pp. 901-913.

    :param left: samples of the left audio channel in 44.1kHz
    :param right: samples of the right audio channel in 44.1kHz
    :return: index of the dominant vocal channel (0 = left, 1 = right)
    '''

    # PARAMETERS
    fs = 44100 # sample rate
    wSize = 2048 # window size in samples
    hSize = 2048 # hop size in samples
    fftSize = 2048 # FFT size
    freqGuitLow = 80.0 # lower bound for guitar band
    freqGuitHigh = 400.0 # upper bound for guitar band
    freqVocLow = 500.0 # lower bound for vocal band
    freqVocHigh = 6000.0 # higher bound for vocal band

    # INIT
    window = hanning(wSize)
    numFrames = int(math.floor(float(len(left))/float(wSize)))
    # bin indices corresponding to freqeuncy band limits
    indGuitLow = int(round((freqGuitLow/fs)*fftSize))
    indGuitHigh = int(round((freqGuitHigh/fs)*fftSize))
    indVocLow = int(round((freqVocLow/fs)*fftSize))
    indVocHigh = int(round((freqVocHigh/fs)*fftSize))

    # frame-wise computation of the spectral band ratio
    sbrL = []
    sbrR = []
    for i in range(0,numFrames-100):
        frameL = left[i*hSize:i*hSize+wSize]
        specL = fft(frameL*window) / fftSize
        specL = abs(specL * conj(specL))
        guitMag = sum(specL[indGuitLow:indGuitHigh],0)
        vocMag = sum(specL[indVocLow:indVocHigh],0)
        sbrL.append(20*math.log10(vocMag/guitMag))
        frameR = right[i*hSize:i*wSize+wSize]
        specR = fft(frameR*window) / fftSize
        specR = abs(specR * conj(specR))
        guitMag = sum(specR[indGuitLow:indGuitHigh],0)
        vocMag = sum(specR[indVocLow:indVocHigh],0)
        sbrR.append(20*math.log10(vocMag/guitMag))

    # select channel based on mean SBR
    if mean(sbrL)>=mean(sbrR):
        ind = 0
    else:
        ind = 1

    return ind

 def calculation(self, data):
     if(isinstance(data[0], types.ListType)):
         freq = [np.abs(fft(l)/len(l))[1:len(l)/2] for l in data]
         freq = [n.tolist() for n in freq]
     else:
         l = data
         freq = np.abs(fft(l)/len(l))[1:len(l)/2]
         #freq = [[v] for v in freq]
         #freq = freq.tolist()
     return freq  # .tolist()

def prob3():
	rate1,sig1 = wavfile.read('chopinw.wav')
	n = sig1.shape[0]
	rate2,sig2 = wavfile.read('balloon.wav')
	m = sig2.shape[0]
	sig1 = sp.append(sig1,sp.zeros((m,2)))
	sig2 = sp.append(sig2,sp.zeros((n,2)))
	f1 = sp.fft(sig1)
	f2 = sp.fft(sig2)
	out = sp.ifft((f1*f2))
	out = sp.real(out)
	scaled = sp.int16(out/sp.absolute(out).max() * 32767)
	wavfile.write('test.wav',rate1,scaled)

def makeIR(wav_in,wav_out,fs,duration,noise=0.025):
    """ measures the response of a speaker (+amp+mic) and build an IR """
    # step 1: full duplex playback and recording. Input: provided sweep wav file
    # output: recorded time response
    ecasound_cmd="ecasound -f:16,1,%i -a:1 -i jack,system,capture " + \
    " -o /tmp/capture.wav -a:2 -i %s -o jack,system -t %i"
    ecasound_cmd=ecasound_cmd%(int(fs),wav_in,int(duration))
    # run capture    
    os.system(ecasound_cmd)
    # load input and capture wave files 
    time.sleep(3)
    f=wave.open(wav_in,'rb')
    len1=f.getnframes()
    #nc1=f.getnchannels()
    #bp1=f.getsampwidth()
    data=f.readframes(len1)
    f.close()
    Y1=scipy.float32(scipy.fromstring(data,dtype='int16'))
    f=wave.open('/tmp/capture.wav','rb')
    len2=f.getnframes()
    #nc1=f.getnchannels()
    #bp1=f.getsampwidth()
    data=f.readframes(len2)
    f.close()    
    Y2=scipy.float32(scipy.fromstring(data,dtype='int16'))
    # truncate and normalize wave file 
    #(or we could pad the shortest to the longest... TODO!)
    minlen = min([len1,len2])
    Y2=Y2[0:minlen]
    Y2=Y2/max(abs(Y2))
    Y1=Y1[0:minlen]
    Y1=Y1/max(abs(Y1))
    # compute frequency response function as ration of both spectra
    FRF=scipy.fft(Y2)/scipy.fft(Y1)
    # compute impulse response as inverse FFT of FRF
    IRraw=scipy.real(scipy.ifft(FRF))
    # get rid of initial lag in IR
    thr=max(abs(IRraw))*noise
    offset=max([0 , min(min(scipy.where(abs(IRraw)>thr)))-5 ])
    IR=IRraw[offset:-1] 
    IRnorm=IR/max(abs(IR))
    # TODO: add post pro options such as low/high pass and decay
    # write output IR
    f = wave.open(wav_out, 'w')
    f.setparams((1, 2, fs, 0, 'NONE', 'not compressed'))
    maxVol=2**15-1.0 #maximum amplitude
    wvData=""
    for i in range(len(IRnorm)):
        wvData+=pack('h', maxVol*IRnorm[i])
    f.writeframes(wvData)
    f.close()

def load_validation_set():
    """
    Output
        a tuple of features: (fft features, mfcc features, mean-std features)
    Description
        extracts three types of features from validation set.
    """
    ffts = dict()
    mfccs = dict()
    mean_stds = dict()

    for i in validation_ids:
        path = './validation/validation.{i}.wav'.format(i=i)

        _, X = read_wav(path)

        # FFT
        fft = np.array(abs(sp.fft(X)[:1000]))
        ffts.update({i: fft})

        # MFCC
        ceps, mspec, spec = mfcc(X)
        num_ceps = len(ceps)
        x = np.mean(ceps[int(num_ceps*1/10):int(num_ceps*9/10)], axis=0)
        mfccs.update({i: x})


        # Mean-Std
        [Fs, x] = audioBasicIO.readAudioFile(path);
        F = audioFeatureExtraction.stFeatureExtraction(x, Fs, 0.050*Fs, 0.025*Fs);
        mean_std = []
        for f in F:
            mean_std.extend([f.mean(), f.std()])
        mean_stds.update({i: np.array(mean_std)})
    return (ffts, mfccs, mean_stds)

 def _fft(self, data):
     data_length = len(data)
     frequencies = scipy.fft(data) / data_length
     frequencies = frequencies[range(data_length / 2)]
     frequencies[0] = 0
     frequencies = np.abs(frequencies)
     return frequencies

def printer(function):
    window[function] = eval("scipy.signal.%s(size)" % function)
    pylab.plot(numpy.abs(window[function]))
    pylab.xlim(0, size - 1)
    pylab.savefig("%s-%03d.png" % (function, speed))
    pylab.close()

    fft[function] = scipy.fft(window[function])
    pylab.loglog(numpy.abs(fft[function])[0 : size / 2 + 1])
    pylab.savefig("%s-fft-%03d.png" % (function, speed))
    pylab.close()

    sinus[function] = scipy.fft(y * window[function])
    pylab.loglog(numpy.abs(sinus[function])[0 : size / 2 + 1])
    pylab.savefig("%s-sinus-%03d.png" % (function, speed))
    pylab.close()

 def get_amplitude(self, signal, l):
     if self.amplitude.has_key(l):
         return self.amplitude[l]
     else:
         amp = sp.absolute(sp.fft(get_frame(signal, self.winsize, l) * self.window))
         self.amplitude[l] = amp
         return amp

def save_plotSpectrum(y,Fs,image_name):
    """
    Plots a Single-Sided Amplitude Spectrum of y(t)
    """
    fig = Figure(linewidth=0.0)
    fig.set_size_inches(fig_width,fig_length, forward=True)
    Figure.subplots_adjust(fig, left = fig_left, right = fig_right, bottom = fig_bottom, top = fig_top, hspace = fig_hspace)
    n = len(y) # length of the signal

    _subplot = fig.add_subplot(2,1,1)        
    print "Fi"
    _subplot.plot(arange(0,n),y)
    xlabel('Time')
    ylabel('Amplitude')
    _subploti_2=fig.add_subplot(2,1,2)
    k = arange(n)
    T = n/Fs
    frq = k/T # two sides frequency range
    frq = frq[range(n/2)] # one side frequency range

    Y = fft(y)/n # fft computing and normalization
    Y = Y[range(n/2)]

    _subplot_2.plot(frq,abs(Y),'r') # plotting the spectrum
    xlabel('Freq (Hz)')
    ylabel('|Y(freq)|')
    print "here"
    canvas = FigureCanvasAgg(fig)
    if '.eps' in outfile_name:
        canvas.print_eps(outfile_name, dpi = 110)
    if '.png' in outfile_name:
        canvas.print_figure(outfile_name, dpi = 110)

def test():
    run_c("fft.c")
    x_re = [float(i) for i in open("x_re")]
    x_im = [float(i) for i in open("x_im")]
    fft_x_re = [float(i) for i in open("fft_x_re")]
    fft_x_im = [float(i) for i in open("fft_x_im")]

    time_complex = [i + (j*1.0) for i, j in zip(x_re, x_im)]
    numpy_complex = s.fft(time_complex)
    numpy_magnitude = n.abs(numpy_complex)

    chips_complex = [i + (j*1.0j) for i, j in zip(fft_x_re, fft_x_im)]
    chips_magnitude = n.abs(chips_complex)

    f, subplot = pyplot.subplots(3, sharex=True)

    pyplot.subplots_adjust(hspace=1.0)
    subplot[0].plot(x_re, 'g')
    subplot[1].plot(numpy_magnitude, 'r')
    subplot[2].plot(chips_magnitude, 'b')
    pyplot.xlim(0, 1023)
    subplot[0].set_title("Time Domain Signal (64 point sine)")
    subplot[1].set_title("Frequency Spectrum - Numpy")
    subplot[2].set_title("Frequency Spectrum - Chips")
    subplot[0].set_xlabel("Sample")
    subplot[1].set_xlabel("Sample")
    subplot[2].set_xlabel("Sample")
    pyplot.savefig("../docs/source/examples/images/example_5.png")
    pyplot.show()

def compute_spectrum(signal, rate):
    n = len(signal)
    k = np.arange(n)
    T = n/rate
    freq = k/T
    Y = sp.fft(signal)/n
    return freq, Y

def returnSpectrum(y,Fs):
    """
    Plots a Single-Sided Amplitude Spectrum of y(t)
    y = wavefile in array form
    Fs = sample frequency, sometimes called sampling rate; likely 44100
    """

    n = len(y) # length of the signal
    k = arange(n) # array from [0, 1, 2, ... n-1]
    T = n/float(Fs) # need to make Fs a float so T can be 0.xx decimal

    if y.ndim > 1: # stereo music has more than one dimension for stereo sound
        y = y[:,0] # taking 1st dimension, slice of matrix

    frq = k/T # two sides frequency range
    # take array, divide by num secs
    frq = frq[range(n/2)] # one side frequency range

    Y = fft(y)/n # fft computing and normalization
    Y = Y[range(n/2)]

    zipped =  zip(frq, abs(Y)) # takes 2 arrays of same size & combines

    frequencies = dict(zipped)
    return frequencies