fit waveforms in any domain
167683
def fit(self):
        
        # solve for estimator of B
        n,p          = np.shape(self._X)
        self._df     = float(n - p)
        self._Cx     = np.linalg.pinv(np.dot(self._X.T,self._X))
        self._Bhat   = np.dot(np.dot(self._Cx, self._X.T), self._A)
        self._Y_rec  = self._compute_prediction(self._X)
62578
def bandpass(ts, low_hz, high_hz, order=3):
    """"""
    orig_ndim = ts.ndim
    if ts.ndim is 1:
        ts = ts[:, np.newaxis]
    channels = ts.shape[1]
    fs = (len(ts) - 1.0) / (ts.tspan[-1] - ts.tspan[0])
    nyq = 0.5 * fs
    low = low_hz/nyq
    high = high_hz/nyq
    b, a = signal.butter(order, [low, high], btype='band')
    if not np.all(np.abs(np.roots(a)) < 1.0):
        raise ValueError('Filter will not be stable with these values.')
    dtype = ts.dtype
    output = np.zeros((len(ts), channels), dtype)
    for i in range(channels):
        output[:, i] = signal.filtfilt(b, a, ts[:, i])
    if orig_ndim is 1:
        output = output[:, 0]
    return Timeseries(output, ts.tspan, labels=ts.labels)
835995
def _check_fit_data(self, X):
        """"""
        X = check_array(X, accept_sparse="csr", dtype=[np.float64, np.float32])
        n_samples, n_features = X.shape
        if X.shape[0] < self.n_clusters:
            raise ValueError(
                "n_samples=%d should be >= n_clusters=%d"
                % (X.shape[0], self.n_clusters)
            )

        for ee in range(n_samples):
            if sp.issparse(X):
                n = sp.linalg.norm(X[ee, :])
            else:
                n = np.linalg.norm(X[ee, :])

            if np.abs(n - 1.) > 1e-4:
                raise ValueError("Data l2-norm must be 1, found {}".format(n))

        return X
49240
def get_waveform_filter_norm(approximant, psd, length, delta_f, f_lower):
    
    if approximant in _filter_norms:
        return _filter_norms[approximant](psd, length, delta_f, f_lower)
    else:
        return None
684536
def Morlet (sampr, freq, width=7.0):
  
  dt = 1. / sampr
  sf = freq / width # sf in Hz
  st = 1. / (2. * np.pi * sf) # st in s
  t = np.arange(-3.5*st, 3.5*st, dt)
  A = 1. / (st * np.sqrt(2.*np.pi)) # A in 1 / s
  # units: 1/s * (exp (s**2 / s**2)) * exp( 1/ s * s)
  return A * np.exp(-t**2. / (2. * st**2.)) * np.exp(1.j * 2. * np.pi * freq * t)
998647
def hilbert(ts):
    """"""
    output = signal.hilbert(signal.detrend(ts, axis=0), axis=0)
    return Timeseries(output, ts.tspan, labels=ts.labels)
53613
def fdfilter(data, *filt, **kwargs):
    
    # parse keyword args
    inplace = kwargs.pop('inplace', False)
    analog = kwargs.pop('analog', False)
    fs = kwargs.pop('sample_rate', None)
    if kwargs:
        raise TypeError("filter() got an unexpected keyword argument '%s'"
                        % list(kwargs.keys())[0])

    # parse filter
    if fs is None:
        fs = 2 * (data.shape[-1] * data.df).to('Hz').value
    form, filt = parse_filter(filt, analog=analog, sample_rate=fs)
    lti = signal.lti(*filt)

    # generate frequency response
    freqs = data.frequencies.value.copy()
    fresp = numpy.nan_to_num(abs(lti.freqresp(w=freqs)[1]))

    # apply to array
    if inplace:
        data *= fresp
        return data
    new = data * fresp
    return new
854186
def calibrate_signal(signal, resp, fs, frange):
    """"""
    # remove dc offset from recorded response (synthesized orignal shouldn't have one)
    dc = np.mean(resp)
    resp = resp - dc

    npts = len(signal)
    f0 = np.ceil(frange[0] / (float(fs) / npts))
    f1 = np.floor(frange[1] / (float(fs) / npts))

    y = resp
    # y = y/np.amax(y) # normalize
    Y = np.fft.rfft(y)

    x = signal
    # x = x/np.amax(x) # normalize
    X = np.fft.rfft(x)

    H = Y / X

    # still issues warning because all of Y/X is executed to selected answers from
    # H = np.where(X.real!=0, Y/X, 1)
    # H[:f0].real = 1
    # H[f1:].real = 1
    # H = smooth(H)

    A = X / H

    return np.fft.irfft(A)
1243246
def _srads2bt(self, data, channel_name):
        """"""
        a__, b__, c__ = BTFIT[channel_name]
        wavenumber = CALIB[self.platform_id][channel_name]["VC"]
        temp = self._tl15(data, wavenumber)

        return a__ * temp * temp + b__ * temp + c__
383186
def fit_transform(self, X, y=None, **fit_params):
        
        y = column_or_1d(X, warn=True)
        _check_numpy_unicode_bug(X)
        self.classes_, X = np.unique(X, return_inverse=True)
        return X
951950
def rfft2d_freqs(h, w):
    """"""

    fy = np.fft.fftfreq(h)[:, None]
    # when we have an odd input dimension we need to keep one additional
    # frequency and later cut off 1 pixel
    if w % 2 == 1:
        fx = np.fft.fftfreq(w)[: w // 2 + 2]
    else:
        fx = np.fft.fftfreq(w)[: w // 2 + 1]
    return np.sqrt(fx * fx + fy * fy)
869194
def fit(self, X):
        
        self.categorical_encoder = self.encoder_class(cols=list(X))
        self.categorical_encoder.fit(X)
        return self
46975
def psd(ts, nperseg=1500, noverlap=1200, plot=True):
    """"""
    ts = ts.squeeze()
    if ts.ndim is 1:
        ts = ts.reshape((-1, 1))
    fs = (len(ts) - 1.0) / (ts.tspan[-1] - ts.tspan[0])
    window = signal.hamming(nperseg, sym=False)
    nfft = max(256, 2**np.int(np.log2(nperseg) + 1))
    freqs, pxx = signal.welch(ts, fs, window, nperseg, noverlap, nfft,
                              detrend='linear', axis=0)
    # Discard estimates for freq bins that are too low for the window size.
    # (require two full cycles to fit within the window)
    index = np.nonzero(freqs >= 2.0*fs/nperseg)[0][0]
    if index > 0:
        freqs = freqs[index:]
        pxx = pxx[index:]
    # Discard estimate for last freq bin as too high for Nyquist frequency:
    freqs = freqs[:-1]
    pxx = pxx[:-1]
    if plot is True:
        _plot_psd(ts, freqs, pxx)
    return freqs, pxx
596719
def sheardec2D(X, shearletsystem):
    """"""
    coeffs = np.zeros(shearletsystem.shearlets.shape, dtype=complex)
    Xfreq = fftshift(fft2(ifftshift(X)))
    for i in range(shearletsystem.nShearlets):
        coeffs[:, :, i] = fftshift(ifft2(ifftshift(Xfreq * np.conj(
                                   shearletsystem.shearlets[:, :, i]))))
    return coeffs.real
864797
def stft_from_frames(frames, window='hann', dtype=np.complex64):
    

    win_length = frames.shape[0]
    n_fft = win_length

    fft_window = filters.get_window(window, win_length, fftbins=True)

    # Reshape so that the window can be broadcast
    fft_window = fft_window.reshape((-1, 1))

    # Pre-allocate the STFT matrix
    stft_matrix = np.empty((int(1 + n_fft // 2), frames.shape[1]),
                           dtype=dtype,
                           order='F')

    # how many columns can we fit within MAX_MEM_BLOCK?
    n_columns = int(util.MAX_MEM_BLOCK / (stft_matrix.shape[0] *
                                          stft_matrix.itemsize))

    for bl_s in range(0, stft_matrix.shape[1], n_columns):
        bl_t = min(bl_s + n_columns, stft_matrix.shape[1])

        # RFFT and Conjugate here to match phase from DPWE code
        stft_matrix[:, bl_s:bl_t] = fft.fft(fft_window *
                                            frames[:, bl_s:bl_t],
                                            axis=0)[:stft_matrix.shape[0]].conj()

    return stft_matrix
516785
def get_interpolated_fd_waveform(dtype=numpy.complex64, return_hc=True,
                                 **params):
    

    def rulog2(val):
        return 2.0 ** numpy.ceil(numpy.log2(float(val)))

    orig_approx = params['approximant']
    params['approximant'] = params['approximant'].replace('_INTERP', '')
    df = params['delta_f']

    if 'duration' not in params:
        duration = get_waveform_filter_length_in_time(**params)
    elif params['duration'] > 0:
        duration = params['duration']
    else:
        err_msg = "Waveform duration must be greater than 0."
        raise ValueError(err_msg)

    #FIXME We should try to get this length directly somehow
    # I think this number should be conservative
    ringdown_padding = 0.5

    df_min = 1.0 / rulog2(duration + ringdown_padding)
    # FIXME: I don't understand this, but waveforms with df_min < 0.5 will chop
    #        off the inspiral when using ringdown_padding - 0.5.
    #        Also, if ringdown_padding is set to a very small
    #        value we can see cases where the ringdown is chopped.
    if df_min > 0.5:
        df_min = 0.5
    params['delta_f'] = df_min
    hp, hc = get_fd_waveform(**params)
    hp = hp.astype(dtype)
    if return_hc:
        hc = hc.astype(dtype)
    else:
        hc = None

    f_end = get_waveform_end_frequency(**params)
    if f_end is None:
        f_end = (len(hp) - 1) * hp.delta_f
    if 'f_final' in params and params['f_final'] > 0:
        f_end_params = params['f_final']
        if f_end is not None:
            f_end = min(f_end_params, f_end)

    n_min = int(rulog2(f_end / df_min)) + 1
    if n_min < len(hp):
        hp = hp[:n_min]
        if hc is not None:
            hc = hc[:n_min]

    offset = int(ringdown_padding * (len(hp)-1)*2 * hp.delta_f)

    hp = interpolate_complex_frequency(hp, df, zeros_offset=offset, side='left')
    if hc is not None:
        hc = interpolate_complex_frequency(hc, df, zeros_offset=offset,
                                           side='left')
    params['approximant'] = orig_approx
    return hp, hc
994634
def freq_filt(bma):
    
    """
    Want to fit linear function to artifact line in freq space,
    Then mask everything near that line at distances of ~5-200 pixels, 
    Or whatever the maximum CCD artifact dimension happens to be, 
    This will depend on scaling - consult CCD map for interval
    """

    #Fill ndv with random data
    bf = malib.randomfill(bma)

    import scipy.fftpack
    f = scipy.fftpack.fft2(bf)
    ff = scipy.fftpack.fftshift(f)

    #Ben suggested a Hahn filter here, remove the low frequency, high amplitude information
    #Then do a second fft?

    #np.log(np.abs(ff))

    #perc = malib.calcperc(np.real(ff), perc=(80, 95))
    #malib.iv(numpy.real(ff), clim=perc)

    #See http://scipy-lectures.github.io/advanced/image_processing/

    #Starting at a,b, compute argmax along vertical axis for restricted range
    #Fit line to the x and y argmax values
    #Mask [argmax[y]-1:argmax[y]+1]

    #Create radial mask
    ff_dim = np.array(ff.shape)
    a,b = ff_dim/2
    n = ff_dim.max()
    y,x = np.ogrid[-a:n-a, -b:n-b]
    r1 = 40 
    r2 = 60 
    ff_mask = np.ma.make_mask(ff)
    radial_mask = (r1**2 <= x**2 + y**2) & (x**2 + y**2 < r2**2)
    #Note issues with rounding indices here
    #Hacked in +1 for testing
    ff_mask[:] = radial_mask[a-ff_dim[0]/2:a+ff_dim[0], b-ff_dim[1]/2:b+1+ff_dim[1]/2]

    #Combine radial and line mask

    #Convert mask to 0-1, then feather

    fm = ff * ff_mask

    #Inverse fft
    bf_filt = scipy.fftpack.ifft2(scipy.fftpack.ifftshift(fm))

    #Apply original mask
    bf_filt = np.ma.masked_array(bf_filt, bma.mask)