Return a fourier domain waveform approximant, using interpolation
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
969310
def _call_numpy(self, x):
        
        assert isinstance(x, np.ndarray)

        if self.halfcomplex:
            return np.fft.rfftn(x, axes=self.axes)
        else:
            if self.sign == '-':
                return np.fft.fftn(x, axes=self.axes)
            else:
                # Need to undo Numpy IFFT scaling
                return (np.prod(np.take(self.domain.shape, self.axes)) *
                        np.fft.ifftn(x, axes=self.axes))
286407
def apply(self, im):
        
        from scipy.ndimage.interpolation import shift
        return shift(im, map(lambda x: -x, self.delta), mode='nearest')
257987
def fft_freqs(n_fft, fs):
    
    return np.arange(0, (n_fft // 2 + 1)) / float(n_fft) * float(fs)
924908
def _ndtri(y):
    
    # approximation for 0 <= abs(z - 0.5) <= 3/8
    P0 = [
        -5.99633501014107895267E1,
        9.80010754185999661536E1,
        -5.66762857469070293439E1,
        1.39312609387279679503E1,
        -1.23916583867381258016E0,
    ]

    Q0 = [
        1.95448858338141759834E0,
        4.67627912898881538453E0,
        8.63602421390890590575E1,
        -2.25462687854119370527E2,
        2.00260212380060660359E2,
        -8.20372256168333339912E1,
        1.59056225126211695515E1,
        -1.18331621121330003142E0,
    ]

    # Approximation for interval z = sqrt(-2 log y ) between 2 and 8
    # i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
    P1 = [
        4.05544892305962419923E0,
        3.15251094599893866154E1,
        5.71628192246421288162E1,
        4.40805073893200834700E1,
        1.46849561928858024014E1,
        2.18663306850790267539E0,
        -1.40256079171354495875E-1,
        -3.50424626827848203418E-2,
        -8.57456785154685413611E-4,
    ]

    Q1 = [
        1.57799883256466749731E1,
        4.53907635128879210584E1,
        4.13172038254672030440E1,
        1.50425385692907503408E1,
        2.50464946208309415979E0,
        -1.42182922854787788574E-1,
        -3.80806407691578277194E-2,
        -9.33259480895457427372E-4,
    ]

    # Approximation for interval z = sqrt(-2 log y ) between 8 and 64
    # i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
    P2 = [
        3.23774891776946035970E0,
        6.91522889068984211695E0,
        3.93881025292474443415E0,
        1.33303460815807542389E0,
        2.01485389549179081538E-1,
        1.23716634817820021358E-2,
        3.01581553508235416007E-4,
        2.65806974686737550832E-6,
        6.23974539184983293730E-9,
    ]

    Q2 = [
        6.02427039364742014255E0,
        3.67983563856160859403E0,
        1.37702099489081330271E0,
        2.16236993594496635890E-1,
        1.34204006088543189037E-2,
        3.28014464682127739104E-4,
        2.89247864745380683936E-6,
        6.79019408009981274425E-9,
    ]

    sign_flag = 1

    if y > (1 - EXP_NEG2):
        y = 1 - y
        sign_flag = 0

    # Shortcut case where we don't need high precision
    # between -0.135 and 0.135
    if y > EXP_NEG2:
        y -= 0.5
        y2 = y ** 2
        x = y + y * (y2 * _polevl(y2, P0, 4) / _p1evl(y2, Q0, 8))
        x = x * ROOT_2PI
        return x

    x = math.sqrt(-2.0 * math.log(y))
    x0 = x - math.log(x) / x

    z = 1.0 / x
    if x < 8.0:                 # y > exp(-32) = 1.2664165549e-14
        x1 = z * _polevl(z, P1, 8) / _p1evl(z, Q1, 8)
    else:
        x1 = z * _polevl(z, P2, 8) / _p1evl(z, Q2, 8)

    x = x0 - x1
    if sign_flag != 0:
        x = -x

    return x
1180377
def function(self,p):
        """"""
        return 0.5 + 0.5*np.sin(p.frequency*2*pi*self.pattern_y + p.phase)
1371562
function fft2d(re, im) {
    var tre = [],
        tim = [],
        i = 0;
    // x-axis
    for (var y = 0; y < _n; y++) {
        i = y * _n;
        for (var x1 = 0; x1 < _n; x1++) {
            tre[x1] = re[x1 + i];
            tim[x1] = im[x1 + i];
        }
        fft1d(tre, tim);
        for (var x2 = 0; x2 < _n; x2++) {
            re[x2 + i] = tre[x2];
            im[x2 + i] = tim[x2];
        }
    }
    // y-axis
    for (var x = 0; x < _n; x++) {
        for (var y1 = 0; y1 < _n; y1++) {
            i = x + y1 * _n;
            tre[y1] = re[i];
            tim[y1] = im[i];
        }
        fft1d(tre, tim);
        for (var y2 = 0; y2 < _n; y2++) {
            i = x + y2 * _n;
            re[i] = tre[y2];
            im[i] = tim[y2];
        }
    }
}
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)
508213
def _fourier(self):
        

        freq_bin_upper = 2000
        freq_bin_lower = 40
        fs = self._metadata['fs']
        Y_transformed = {}
        for key in self.Y_dict.keys():
            # normalize by fs, bins have units strain/Hz
            fs = self._metadata["fs"]
            Y_transformed[key] = (1./fs)*np.fft.fft(self.Y_dict[key])[freq_bin_lower:freq_bin_upper]
        self.Y_transformed = Y_transformed
        # add in transform metadata
        self._metadata["dF"] = 1./self._metadata["T"]

        # because we are in the fourier domain, we will need the psd
        self.psd = load_psd()[freq_bin_lower:freq_bin_upper]
        dF = self._metadata['dF']
        self.sigma = convert_psd_to_sigma(self.psd, dF)
801789
def apply_magnitude_interpolation(self, mag, iml_table):
        
        # do not allow "mag" to exceed maximum table magnitude
        if mag > self.m_w[-1]:
            mag = self.m_w[-1]

        # Get magnitude values
        if mag < self.m_w[0] or mag > self.m_w[-1]:
            raise ValueError("Magnitude %.2f outside of supported range "
                             "(%.2f to %.2f)" % (mag,
                                                 self.m_w[0],
                                                 self.m_w[-1]))
        # It is assumed that log10 of the spectral acceleration scales
        # linearly (or approximately linearly) with magnitude
        m_interpolator = interp1d(self.m_w, numpy.log10(iml_table), axis=1)
        return 10.0 ** m_interpolator(mag)
1047471
def inverse(self):
        
        return WaveletTransformInverse(
            range=self.domain, wavelet=self.pywt_wavelet, nlevels=self.nlevels,
            pad_mode=self.pad_mode, pad_const=self.pad_const, impl=self.impl,
            axes=self.axes)
905814
def _fftconv(a, b, axes=(0, 1)):
    """"""

    if cp.isrealobj(a) and cp.isrealobj(b):
        fft = cp.fft.rfftn
        ifft = cp.fft.irfftn
    else:
        fft = cp.fft.fftn
        ifft = cp.fft.ifftn
    dims = cp.maximum(cp.asarray([a.shape[i] for i in axes]),
                      cp.asarray([b.shape[i] for i in axes]))
    dims = [int(d) for d in dims]
    af = fft(a, dims, axes)
    bf = fft(b, dims, axes)
    return ifft(af * bf, dims, axes)
203711
def _ncc_c_2dim(x, y):
    
    den = np.array(norm(x, axis=1) * norm(y))
    den[den == 0] = np.Inf
    x_len = x.shape[-1]
    fft_size = 1 << (2*x_len-1).bit_length()
    cc = ifft(fft(x, fft_size) * np.conj(fft(y, fft_size)))
    cc = np.concatenate((cc[:,-(x_len-1):], cc[:,:x_len]), axis=1)
    return np.real(cc) / den[:, np.newaxis]
172685
def interpolation_bilinear(x, y, x1, x2, y1, y2, z11, z21, z22, z12):
    '''
    
    '''
    t = (x - x1) / (x2 - x1)
    s = (y - y1) / (y2 - y1)
    v1 = (1.0 - t) * (1.0 - s) * z11
    v2 = t * (1.0 - s) * z21
    v3 = t * s * z22
    v4 = (1.0 - t) * s * z12
    ret = v1 + v2 + v3 + v4
    return ret
222120
def fourier_series(x, *a):
    
    output = 0
    output += a[0]/2
    w = a[1]
    for n in range(2, len(a), 2):
        n_ = n/2
        val1 = a[n]
        val2 = a[n+1]
        output += val1*np.sin(n_*x*w)
        output += val2*np.cos(n_*x*w)
    return output
654557
def get_waveform_filter(out, template=None, **kwargs):
    
    n = len(out)

    input_params = props(template, **kwargs)

    if input_params['approximant'] in filter_approximants(_scheme.mgr.state):
        wav_gen = filter_wav[type(_scheme.mgr.state)]
        htilde = wav_gen[input_params['approximant']](out=out, **input_params)
        htilde.resize(n)
        htilde.chirp_length = get_waveform_filter_length_in_time(**input_params)
        htilde.length_in_time = htilde.chirp_length
        return htilde

    if input_params['approximant'] in fd_approximants(_scheme.mgr.state):
        wav_gen = fd_wav[type(_scheme.mgr.state)]

        duration = get_waveform_filter_length_in_time(**input_params)
        hp, _ = wav_gen[input_params['approximant']](duration=duration,
                                               return_hc=False, **input_params)

        hp.resize(n)
        out[0:len(hp)] = hp[:]
        hp = FrequencySeries(out, delta_f=hp.delta_f, copy=False)

        hp.length_in_time = hp.chirp_length = duration
        return hp

    elif input_params['approximant'] in td_approximants(_scheme.mgr.state):
        wav_gen = td_wav[type(_scheme.mgr.state)]
        hp, _ = wav_gen[input_params['approximant']](**input_params)
        # taper the time series hp if required
        if 'taper' in input_params.keys() and \
                input_params['taper'] is not None:
            hp = wfutils.taper_timeseries(hp, input_params['taper'],
                                          return_lal=False)
        return td_waveform_to_fd_waveform(hp, out=out)

    else:
        raise ValueError("Approximant %s not available" %
                            (input_params['approximant']))
698783
def notch(ts, freq_hz, bandwidth_hz=1.0):
    
    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
    freq = freq_hz/nyq
    bandwidth = bandwidth_hz/nyq
    R = 1.0 - 3.0*(bandwidth/2.0)
    K = ((1.0 - 2.0*R*np.cos(np.pi*freq) + R**2) /
         (2.0 - 2.0*np.cos(np.pi*freq)))
    b, a = np.zeros(3), np.zeros(3)
    a[0] = 1.0
    a[1] = -2.0*R*np.cos(np.pi*freq)
    a[2] = R**2
    b[0] = K
    b[1] = -2*K*np.cos(np.pi*freq)
    b[2] = K
    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)