def demod_QAM16(QAM, t, f0=1800, fs=48000):
    r = QAM*np.cos(2*np.pi*f0*t)
    i = -QAM*np.sin(2*np.pi*f0*t)

    #plot(r+i)
    num_taps = 100
    lp = signal.firwin(num_taps, np.pi*f0/4,nyq=fs/2.0)
    r_lp = signal.fftconvolve(lp,r)
    i_lp = signal.fftconvolve(lp, i)

    #fig = figure(figsize = (16,4))
    #frange = np.linspace(-fs/2,fs/2,len(r))
    #frange_filt = np.linspace(-fs/2,fs/2,len(r_lp))
    #plt.plot(frange_filt, abs(fft.fftshift(fft.fft(lp))))
    '''
    ylim(-3,3)

    fig = figure(figsize = (16,4))
    plt.plot(frange, abs(fft.fftshift(fft.fft(i))))
    plt.plot(frange_filt, abs(fft.fftshift(fft.fft(i_lp))))

    #r_env = abs(r_lp)
    #i_env = abs(i_lp)
    '''
    r_lp = r_lp[num_taps/2:-num_taps/2+1]
    i_lp = i_lp[num_taps/2:-num_taps/2+1]
    return r_lp, i_lp

def detrend(EEG):
    window_size = 207*10
    filt = np.ones((window_size,))/float(window_size)
    trend0 = signal.fftconvolve(EEG[:,0], filt, 'same')
    trend1 = signal.fftconvolve(EEG[:,1], filt, 'same')
    trend = np.vstack([trend0,trend1]).T
    return EEG-trend

def golayIR(respa, respb, a, b):
    # Comptute impulse response h for Signle Channel answers a and b
    L = len(a)
    h = np.array(np.zeros(respa.shape))
    h = fftconvolve(a[-1::-1], respa, mode="same") + fftconvolve(b[-1::-1], respb, mode="same")
    h = h / (2 * L)
    return h

def conv_mul(lin_op, rh_val, transpose=False, is_abs=False):
    """Multiply by a convolution operator.

    arameters
    ----------
    lin_op : LinOp
        The root linear operator.
    rh_val : NDArray
        The vector being convolved.
    transpose : bool
        Is the transpose of convolution being applied?
    is_abs : bool
        Is the absolute value of convolution being applied?

    Returns
    -------
    NumPy NDArray
        The convolution.
    """
    constant = mul(lin_op.data, {}, is_abs)
    # Convert to 2D
    constant, rh_val = map(intf.from_1D_to_2D, [constant, rh_val])
    if transpose:
        constant = np.flipud(constant)
        # rh_val always larger than constant.
        return fftconvolve(rh_val, constant, mode='valid')
    else:
        # First argument must be larger.
        if constant.size >= rh_val.size:
            return fftconvolve(constant, rh_val, mode='full')
        else:
            return fftconvolve(rh_val, constant, mode='full')

def energy(traces, duration, dt):
    """
    Compute an mean-squared energy measurement for each point of a
    seismic section.
        
    :param traces: The data array to use for calculating MS energy.
                   Must be 1D or 2D numpy array.
    :param duration: the time duration of the window (in seconds)
    :param dt: the sample interval of the data (in seconds) 
    :returns An array the same dimensions as the input array.
    """

    energy_data = numpy.zeros(traces.shape)
    signal = traces * traces
    n_samples = int(duration / dt)

    window = numpy.ones(n_samples)

    if (len(signal.shape)) == 1:
        ## Compute the sliding average using a convolution
        energy_data = fftconvolve(signal, window, mode="same") / n_samples

    elif len(signal.shape) == 2:
        for trace in range(signal.shape[1]):
            energy_data[:, trace] = fftconvolve(signal[:, trace], window, mode="same")

    else:
        raise ValueError("Array must be 1D or 2D")

    return energy_data

    def demodulate2(self,samples):
        # DEMODULATION CODE

        # LIMITER goes here

        # low pass & down sampling
        h = signal.firwin(128,80000,nyq=1.2e5)
        lp_samples = signal.fftconvolve(samples, h)

    # polar discriminator

        A = lp_samples[1:lp_samples.size]
        B = lp_samples[0:lp_samples.size-1]

        dphase = ( A * np.conj(B) ) / np.pi

        dphase.resize(dphase.size+1)
        dphase[dphase.size-1] = dphase[dphase.size-2]

        h = signal.firwin(128,16000,nyq=1.2e5)
        rebuilt = signal.fftconvolve(dphase,h)

        output = rebuilt[::self.decim_r2]

        output = self.lowpass(output, self.audioFilterSize)

        return np.real(output)

def getData():    #function called to get data
    
    startTime = time.time()
    raw650 = np.array(session[s1Vals])
    raw950 = np.array(session[s2Vals])
  
 #   while True:

   #     if time.time() - startTime >= 5:
    startTime = time.time()
    print('got data')
    working950 = reject_outliers(raw950)
    working650 = reject_outliers(raw650)
    sig950 = np.std(working950)
    sig650 = np.std(working650)
    print(sig650)
    window950 = signal.general_gaussian(51, p=1.5, sig= sig950)
    filtered950 = signal.fftconvolve(window950, working950)
    filtered950 = (np.average(working950) / np.average(filtered950)) * filtered950
    window650 = signal.general_gaussian(51, p=1.5, sig= sig650)
    filtered650 = signal.fftconvolve(window650, working650)
    filtered650 = (np.average(working650) / np.average(filtered650)) * filtered650

  #  filtered = np.roll(filtered, -25)
 #    plt.plot(working950)
 # #   plt.plot(window950)
 #    plt.plot(filtered950)
 #    plt.plot(raw650)
 #    plt.show()
    
    print(filtered950)
    
    print(working950)
    concentration(filtered650,filtered950)

def acovf_fft(x, demean=True):
    '''autocovariance function with call to fftconvolve, biased

    Parameters
    ----------
    x : array_like
        timeseries, signal
    demean : boolean
        If true, then demean time series

    Returns
    -------
    acovf : array
        autocovariance for data, same length as x

    might work for nd in parallel with time along axis 0

    '''
    from scipy import signal
    x = np.asarray(x)

    if demean:
        x = x - x.mean()

    signal.fftconvolve(x,x[::-1])[len(x)-1:len(x)+10]/x.shape[0]

def SIC_method(X,Y,order=100):   
    #low-passing to take LFP only
    
    h_for=AR_fit(X,Y,order)
    y_new=signal.fftconvolve(h_for,X)
    h_back=AR_fit(Y,X,order)
    x_new=signal.fftconvolve(h_back,Y)

    #Sx=welch(x_new,nperseg=1000)[1]
    #Sy=welch(y_new,nperseg=1000)[1]

#    Sy=welch(Y,nperseg=1000)[1]
#    Sx=welch(X,nperseg=1000)[1]
#
#    X_Y=delta_estimator_3(Sy/Sx,Sx)
#    Y_X=delta_estimator_3(Sx/Sy,Sy)
            
    #mask1=Sy!=0
    #mask2=Sx[mask1]!=0
    #X_Y=eval('delta_estimator_'+str(method_no))(Sy[mask1][mask2][1:-1]/Sx[mask1][mask2][1:-1],Sx[mask1][mask2][1:-1])
    #Y_X=eval('delta_estimator_'+str(method_no))(Sx[mask1][mask2][1:-1]/Sy[mask1][mask2][1:-1],Sy[mask1][mask2][1:-1])
    #X_Y=eval('delta_estimator_'+str(method_no))(Sy[mask1][mask2][1:-1]/Sx[mask1][mask2][1:-1],Sx[mask1][mask2][1:-1])
    #Y_X=eval('delta_estimator_'+str(method_no))(Sx[mask1][mask2][1:-1]/Sy[mask1][mask2][1:-1],Sy[mask1][mask2][1:-1])

    X_Y=np.var(y_new)/float(np.sum(h_for**2)*np.var(X))
    Y_X=np.var(x_new)/float(np.sum(h_back**2)*np.var(Y))
    
    return X_Y,Y_X

def nc_afskDemod(sig, TBW=2.0, N=74, fs=48000):
    bw = float(TBW*fs)/N
    t = np.r_[:N]/float(fs)
    h = signal.firwin(N,bw/2.0,nyq=float(fs)/2)
    b0=h*np.exp(2*np.pi*1.0j*1200.0*t)
    b1=h*np.exp(2*np.pi*1.0j*2200.0*t)
    return np.abs(signal.fftconvolve(sig,b1)) - np.abs(signal.fftconvolve(sig,b0))

def dwt(X, L, H, n):
    '''
    Compute the discrete wavelet transform of f with respect to 
    the wavelet filters lo and hi.
    Inputs:
        X -- numpy array corresponding to the signal
        L -- numpy array giving the lo-pass filter
        H -- numpy array giving the hi-pass filter
        n -- integer, giving what level of decomposition
    Returns:
        list of the form [A, D1, D2, ..., Dn] where each entry
        is a numpy array. These are the approximation frame (A)
        and the detail coefficients.
    '''
    coeffs = []
    A = X
    i=0
    while i < n:
        D = fftconvolve(A,H)[1::2]
        A = fftconvolve(A,L)[1::2]
        coeffs.append(D)
        i += 1
    coeffs.append(A)
    coeffs.reverse()
    return coeffs

def energy(traces, duration, dt=1):
    """
    Compute an mean-squared energy measurement for each point of a
    seismic section.

    :param traces: The data array to use for calculating MS energy.
                   Must be 1D or 2D numpy array.
    :param duration: the time duration of the window (in seconds), or
                     samples if dt=1.
    :param dt: the sample interval of the data (in seconds). Defaults
               to 1 so duration can be in samples.
    :returns: An array the same dimensions as the input array.
    """

    energy_data = np.zeros(traces.shape)
    signal = traces * traces
    n_samples = int(duration / dt)

    window = np.ones(n_samples)

    if np.ndim(signal) == 1:
        # Compute the sliding average using a convolution
        energy_data = fftconvolve(signal, window, mode='same') \
                     / n_samples

    elif np.ndim(signal) == 2:
        for trace in range(signal.shape[1]):
            energy_data[:, trace] = (fftconvolve(signal[:, trace],
                                                 window,
                                                 mode='same'))

    else:
        raise ValueError('Array must be 1D or 2D')

    return energy_data

  def __init__ (self, var, saxis, kernel, fft):
  # {{{
    ''' __init__()'''
    import numpy as np

    assert len(kernel) <= var.shape[saxis], 'Kernel must not be longer than dimension being smoothed.'

    # Construct new variable
    self.saxis = saxis
    self.var = var
    self.kernel = kernel
    self.fft = fft
    self.klen = len(kernel)

    # Normalize and reshape kernel
    self.kernel /= np.sum(self.kernel)
    self.kernel.shape = [self.klen if i == saxis else 1 for i in range(var.naxes)]

    # Determine which convolution function to use
    from scipy import signal as sg
    tdata = np.ones(len(kernel), 'd')
    if self.fft:
      try:
        sg.fftconvolve(tdata, kernel, 'same', old_behaviour=False)
        self._convolve = lambda x, y, z: sg.fftconvolve(x, y, z, old_behaviour=False)
      except TypeError:
        self._convolve = sg.fftconvolve
    else:
      try:
        sg.convolve(tdata, kernel, 'same', old_behaviour=False)
        self._convolve = lambda x, y, z: sg.convolve(x, y, z, old_behaviour=False)
      except TypeError:
        self._convolve = sg.convolve

    Var.__init__(self, var.axes, var.dtype, name=var.name, atts=var.atts, plotatts=var.plotatts)

    def contexttrack(self):
        print("[%s] tracking" % QThread.currentThread().objectName())
        if self._isRunning:
            self.wait.emit()
            t_trackerstart = time.time()
            print('\t start context tracking')
            dx, dy = np.random.randint(0, 20, size=2)
            self.data = np.roll(self.data, dx, axis=0)
            self.data = np.roll(self.data, dy, axis=1)
            self.data[0:dy, :] = 0.
            self.data[:, 0:dx] = 0.
            self.correlation = sig.fftconvolve(self.data, self.referencedata[::-1, ::-1], 'same')
            maxposy, maxposx = np.unravel_index(self.correlation.argmax(), self.correlation.shape)
            xcorrect = int(np.round(self.resolution1/2-maxposx))
            ycorrect = int(np.round(self.resolution2/2-maxposy))
            self.data = np.roll(self.data, ycorrect, axis=0)
            self.data = np.roll(self.data, xcorrect, axis=1)

            # check result
            self.correlation1 = sig.fftconvolve(self.data, self.referencedata[::-1, ::-1], 'same')
            maxposx, mayposy = np.unravel_index(self.correlation1.argmax(), self.correlation1.shape)

            self.update.emit(self.data, self.correlation, xcorrect, ycorrect, self.correlation.max(), time.strftime('%H:%M:%S')+' - ok')

            print('tracking took ', time.time()-t_trackerstart, 's')
            print('\t context tracking done')

            self.goon.emit()

    def fir_filter(self, fir_ac=None, fir_dc=None, f_ac=None, f_dc=None,
                   a_ac=10, a_dc=10, alpha=None, filter_name=None, **kwargs):
        """Apply filters to generate the lock-in and dc components of phi"""

        if filter_name == 'bessel_matched':
            N_pts = kwargs.get('N_pts', int(self.ks / self.k0_dc * 6))
            dec = kwargs.get('dec', 32)
            n_pts_eval_fir = kwargs.get('n_pts_eval_fir', 2**16)
            window = kwargs.get('window', 'hann')

            fir_ac, fir_dc = _matched_filters(self.ks, self.x_m, N_pts, dec, window,
                                              n_pts_eval_fir)

            self.fir_ac = fir_ac
            self.fir_dc = fir_dc
        else:
            if fir_ac is None:
                if f_ac is None and alpha is None:
                    f_ac = self.fx * 0.5
                elif alpha is not None:
                    f_ac = self.v_tip/self.x_m * alpha
                self.fir_ac = signal.firwin(self.fs / (f_ac) * a_ac,
                                            f_ac, nyq=0.5 * self.fs,
                                            window='blackman')
            else:
                self.fir_ac = fir_ac

            if fir_dc is None:
                if f_dc is None and alpha is None:
                    f_dc = self.fx * 0.5
                elif alpha is not None:
                    f_dc = self.v_tip/self.x_m * alpha
                self.fir_dc = signal.firwin(self.fs/(f_dc) * a_dc,
                                            f_dc, nyq=0.5*self.fs,
                                            window='blackman')
            else:
                self.fir_dc = fir_dc

        indices = np.arange(self.phi.size)
        fir_ac_size = self.fir_ac.size
        fir_dc_size = self.fir_dc.size

        fir_max_size = max(fir_ac_size, fir_dc_size)

        self.m = indices[fir_max_size//2: -fir_max_size//2]
        self.tm = self.t[self.m]

        self._lock = np.exp(np.pi * 2j * self.fx * self.t)

        self.phi_lock = signal.fftconvolve(self.phi * self._lock * 2,
                                           self.fir_ac,
                                           mode='same')

        self.V_lock = self.phi_lock

        self.phi_lock_a = np.abs(self.phi_lock)
        self.phi_lock_phase = np.angle(self.phi_lock)

        self.phi_dc = signal.fftconvolve(self.phi, self.fir_dc, mode='same')
        self.V_dc = self.phi_dc

def convolutionFilter(image,vORh):
    if vORh == 'v':
        # return sg.convolve(image, [[1,-1]], "valid")         # Vertical Data
        return sg.fftconvolve(image, [[1,-1]], "valid")         # Vertical Data
    else:
        # return sg.convolve(image, [[1],[-1]], "valid")       # Horizontal Data
        return sg.fftconvolve(image, [[1],[-1]], "valid")       # Horizontal Data

    def _transform(self, data):
        """
        Do the transform itself.

        The transform is made by using `scipy.signal.fftconvolve`.

        TODO: document.

        Parameters
        ----------
        data : `~gammapy.detect.CWTData`
            Images for transform.
        """
        from scipy.signal import fftconvolve

        total_background = data._model + data._background + data._approx
        excess = data._counts - total_background
        log.debug('Excess sum: {0:.4f}'.format(excess.sum()))
        log.debug('Excess max: {0:.4f}'.format(excess.max()))

        log.debug('Computing transform and error')
        for idx_scale, kern in self.kernels.kern_base.items():
            data._transform_3d[idx_scale] = fftconvolve(excess, kern, mode='same')
            data._error[idx_scale] = np.sqrt(fftconvolve(total_background, kern ** 2, mode='same'))
        log.debug('Error sum: {0:.4f}'.format(data._error.sum()))
        log.debug('Error max: {0:.4f}'.format(data._error.max()))

        log.debug('Computing approx and approx_bkg')
        data._approx = fftconvolve(data._counts - data._model - data._background,
                                   self.kernels.kern_approx, mode='same')
        data._approx_bkg = fftconvolve(data._background, self.kernels.kern_approx, mode='same')
        log.debug('Approximate sum: {0:.4f}'.format(data._approx.sum()))
        log.debug('Approximate background sum: {0:.4f}'.format(data._approx_bkg.sum()))

def ssim(img1, img2, cs_map=False):
    """Return the Structural Similarity Map corresponding to input images img1 
    and img2 (images are assumed to be uint8)
    
    This function attempts to mimic precisely the functionality of ssim.m a 
    MATLAB provided by the author's of SSIM
    https://ece.uwaterloo.ca/~z70wang/research/ssim/ssim_index.m
    """
    img1 = img1.astype(numpy.float64)
    img2 = img2.astype(numpy.float64)
    size = 11
    sigma = 1.5
    window = gauss.fspecial_gauss(size, sigma)
    K1 = 0.01
    K2 = 0.03
    L = 255 #bitdepth of image
    C1 = (K1*L)**2
    C2 = (K2*L)**2
    mu1 = signal.fftconvolve(window, img1, mode='valid')
    mu2 = signal.fftconvolve(window, img2, mode='valid')
    mu1_sq = mu1*mu1
    mu2_sq = mu2*mu2
    mu1_mu2 = mu1*mu2
    sigma1_sq = signal.fftconvolve(window, img1*img1, mode='valid') - mu1_sq
    sigma2_sq = signal.fftconvolve(window, img2*img2, mode='valid') - mu2_sq
    sigma12 = signal.fftconvolve(window, img1*img2, mode='valid') - mu1_mu2
    if cs_map:
        return (((2*mu1_mu2 + C1)*(2*sigma12 + C2))/((mu1_sq + mu2_sq + C1)*
                    (sigma1_sq + sigma2_sq + C2)), 
                (2.0*sigma12 + C2)/(sigma1_sq + sigma2_sq + C2))
    else:
        return ((2*mu1_mu2 + C1)*(2*sigma12 + C2))/((mu1_sq + mu2_sq + C1)*
                    (sigma1_sq + sigma2_sq + C2))

def nc_afsk1200Demod(sig, fs=48000.0, TBW=2.0):
    #  non-coherent demodulation of afsk1200
    # function returns the NRZ (without rectifying it)
    # 
    # sig  - signal
    # baud - The bitrate. Default 1200
    # fs   - sampling rate in Hz
    # TBW  - TBW product of the filters
    #
    # Returns:
    #     NRZ  
    # your code here
    taps = fs/1200-1
    bandpass = signal.firwin(taps, 1200, nyq=fs/2)
    spacepass = bandpass * np.exp(1j*2*np.pi*1200*np.r_[0.0:taps]/fs)
    markpass = bandpass * np.exp(1j*2*np.pi*3600*np.r_[0.0:taps]/fs)
    spaces = signal.fftconvolve(sig, spacepass, mode='same')
    marks = signal.fftconvolve(sig, markpass, mode='same')

    analog = np.abs(spaces)-np.abs(marks)
    lowpass = signal.firwin(taps, 2400*1.2, nyq=fs/2)
    filtered = signal.fftconvolve(analog, lowpass, mode='same')
    NRZ = filtered
    
    return NRZ

def run():
    print_complex = get_print_complex()
    convolutionCPU = get_convolution_cpu()
    check_results = get_check_results()

    #data = np.ones((3,3)).astype('complex64')
    data = np.asfortranarray(np.random.randn(3,3).astype('complex64'))
    #kernel = np.ones((3,3)).astype('complex64')
    kernel = np.asfortranarray(np.random.randn(3,3).astype('complex64'))
    result = np.asfortranarray(np.zeros_like(data).astype('complex64'))

    convolutionCPU(_get_float2_ptr(result), _get_float2_ptr(data), _get_float2_ptr(kernel), data.shape[1], data.shape[0], kernel.shape[1], kernel.shape[0], 1, 6)

    print
    print kernel
    print
    print data
    print

    s1 = np.array(data.shape)
    s2 = np.array(kernel.shape)

    print result
    print 
    print fftconvolve(data.real, kernel.real, mode='full').astype('complex64')