Phase of the analytic signal, using the Hilbert transform
765764
def hilbert_phase(ts):
    """"""
    output = np.angle(signal.hilbert(signal.detrend(ts, axis=0), axis=0))
    return Timeseries(output, ts.tspan, labels=ts.labels)
14402
public void transform(double[] a) {
        int n = a.length;

        if (!Math.isPower2(n)) {
            throw new IllegalArgumentException("The data vector size is not a power of 2.");
        }

        if (n < ncof) {
            throw new IllegalArgumentException("The data vector size is less than wavelet coefficient size.");
        }

        for (int nn = n; nn >= ncof; nn >>= 1) {
            forward(a, nn);
        }
    }
305581
def get_layer_timing_signal_sinusoid_1d(channels, layer, num_layers):
  

  signal = get_timing_signal_1d(num_layers, channels)
  layer_signal = tf.expand_dims(signal[:, layer, :], axis=1)

  return layer_signal
867693
def _analytical_forbild_phantom(resolution, ear):
    
    sha = 0.2 * np.sqrt(3)
    y016b = -14.294530834372887
    a16b = 0.443194085308632
    b16b = 3.892760834372886

    E = [[-4.7, 4.3, 1.79989, 1.79989, 0, 0.010, 0],  # 1
         [4.7, 4.3, 1.79989, 1.79989, 0, 0.010, 0],  # 2
         [-1.08, -9, 0.4, 0.4, 0, 0.0025, 0],  # 3
         [1.08, -9, 0.4, 0.4, 0, -0.0025, 0],  # 4
         [0, 0, 9.6, 12, 0, 1.800, 0],  # 5
         [0, 8.4, 1.8, 3.0, 0, -1.050, 0],  # 7
         [1.9, 5.4, 0.41633, 1.17425, -31.07698, 0.750, 0],  # 8
         [-1.9, 5.4, 0.41633, 1.17425, 31.07698, 0.750, 0],  # 9
         [-4.3, 6.8, 1.8, 0.24, -30, 0.750, 0],  # 10
         [4.3, 6.8, 1.8, 0.24, 30, 0.750, 0],  # 11
         [0, -3.6, 1.8, 3.6, 0, -0.005, 0],  # 12
         [6.39395, -6.39395, 1.2, 0.42, 58.1, 0.005, 0],  # 13
         [0, 3.6, 2, 2, 0, 0.750, 4],  # 14
         [0, 9.6, 1.8, 3.0, 0, 1.800, 4],  # 15
         [0, 0, 9.0, 11.4, 0, 0.750, 3],  # 16a
         [0, y016b, a16b, b16b, 0, 0.750, 1],  # 16b
         [0, 0, 9.0, 11.4, 0, -0.750, ear],  # 6
         [9.1, 0, 4.2, 1.8, 0, 0.750, 1]]  # R_ear
    E = np.array(E)

    # generate the air cavities in the right ear
    cavity1 = np.arange(8.8, 5.6, -0.4)[:, None]
    cavity2 = np.zeros([9, 1])
    cavity3_7 = np.ones([53, 1]) * [0.15, 0.15, 0, -1.800, 0]

    for j in range(1, 4):
        kj = 8 - 2 * int(np.floor(j / 3))
        dj = 0.2 * int(np.mod(j, 2))

        cavity1 = np.vstack((cavity1,
                             cavity1[0:kj] - dj,
                             cavity1[0:kj] - dj))
        cavity2 = np.vstack((cavity2,
                             j * sha * np.ones([kj, 1]),
                             -j * sha * np.ones([kj, 1])))

    E_cavity = np.hstack((cavity1, cavity2, cavity3_7))

    # generate the left ear (resolution pattern)
    x0 = -7.0
    y0 = -1.0
    d0_xy = 0.04

    d_xy = [0.0357, 0.0312, 0.0278, 0.0250]
    ab = 0.5 * np.ones([5, 1]) * d_xy
    ab = ab.T.ravel()[:, None] * np.ones([1, 4])
    abr = ab.T.ravel()[:, None]

    leftear4_7 = np.hstack([abr, abr, np.ones([80, 1]) * [0, 0.75, 0]])

    x00 = np.zeros([0, 1])
    y00 = np.zeros([0, 1])
    for i in range(1, 5):
        y00 = np.vstack((y00,
                         (y0 + np.arange(0, 5) * 2 * d_xy[i - 1])[:, None]))
        x00 = np.vstack((x00,
                         (x0 + 2 * (i - 1) * d0_xy) * np.ones([5, 1])))

    x00 = x00 * np.ones([1, 4])
    x00 = x00.T.ravel()[:, None]
    y00 = np.vstack([y00, y00 + 12 * d0_xy,
                     y00 + 24 * d0_xy, y00 + 36 * d0_xy])

    leftear = np.hstack([x00, y00, leftear4_7])
    C = [[1.2, 1.2, 0.27884, 0.27884, 0.60687, 0.60687, 0.2,
          0.2, -2.605, -2.605, -10.71177, y016b + 10.71177, 8.88740, -0.21260],
         [0, 180, 90, 270, 90, 270, 0,
          180, 15, 165, 90, 270, 0, 0]]
    C = np.array(C)

    if not resolution and not ear:
        phantomE = E[:17, :]
        phantomC = C[:, :12]
    elif not resolution and ear:
        phantomE = np.vstack([E, E_cavity])
        phantomC = C
    elif resolution and not ear:
        phantomE = np.vstack([leftear, E[:17, :]])
        phantomC = C[:, :12]
    else:
        phantomE = np.vstack([leftear, E, E_cavity])
        phantomC = C

    return phantomE, phantomC
166895
def cwt_haar(x: np.ndarray, scale=10):
    
    next_power_two = 2 ** int(np.log2(len(x)))

    x = x[0:next_power_two]
    num_data = len(x)

    # get FFT of x (eq. (3) in paper)
    x_hat = np.fft.fft(x)

    # Get omega (eq. (5) in paper)
    f = (2.0 * np.pi / num_data)
    omega = f * np.concatenate((np.arange(0, num_data // 2), np.arange(num_data // 2, num_data) * -1))

    # get psi hat (eq. (6) in paper)
    psi_hat = np.sqrt(2.0 * np.pi * scale) * normalized_haar_wavelet(scale * omega, scale)

    # get W (eq. (4) in paper)
    W = np.fft.ifft(x_hat * psi_hat)

    return W[2 * scale:-2 * scale]
786495
def dephasing_operators(p):
    
    k0 = np.eye(2) * np.sqrt(1 - p / 2)
    k1 = np.sqrt(p / 2) * Z
    return k0, k1
878896
def fftconvolve(in1, in2, mode="full", axis=None):
    
    s1 = np.array(in1.shape)
    s2 = np.array(in2.shape)
    complex_result = (np.issubdtype(in1.dtype, np.complexfloating) or
                      np.issubdtype(in2.dtype, np.complexfloating))

    if axis is None:
        size = s1 + s2 - 1
        fslice = tuple([slice(0, int(sz)) for sz in size])
    else:
        equal_shapes = s1 == s2
        # allow equal_shapes[axis] to be False
        equal_shapes[axis] = True
        assert equal_shapes.all(), 'Shape mismatch on non-convolving axes'
        size = s1[axis] + s2[axis] - 1
        fslice = [slice(l) for l in s1]
        fslice[axis] = slice(0, int(size))
        fslice = tuple(fslice)

    # Always use 2**n-sized FFT
    fsize = 2 ** int(np.ceil(np.log2(size)))
    if axis is None:
        IN1 = fftpack.fftn(in1, fsize)
        IN1 *= fftpack.fftn(in2, fsize)
        ret = fftpack.ifftn(IN1)[fslice].copy()
    else:
        IN1 = fftpack.fft(in1, fsize, axis=axis)
        IN1 *= fftpack.fft(in2, fsize, axis=axis)
        ret = fftpack.ifft(IN1, axis=axis)[fslice].copy()
    del IN1
    if not complex_result:
        ret = ret.real
    if mode == "full":
        return ret
    elif mode == "same":
        if np.product(s1, axis=0) > np.product(s2, axis=0):
            osize = s1
        else:
            osize = s2
        return signaltools._centered(ret, osize)
    elif mode == "valid":
        return signaltools._centered(ret, abs(s2 - s1) + 1)
1164261
def mixing_phases(U):
    """"""
    f = {}
    # angles
    f['t13'] = asin(abs(U[0,2]))
    if U[0,0] == 0:
        f['t12'] = pi/2
    else:
        f['t12'] = atan(abs(U[0,1])/abs(U[0,0]))
    if U[2,2] == 0:
        f['t23'] = pi/2
    else:
        f['t23'] = atan(abs(U[1,2])/abs(U[2,2]))
    s12 = sin(f['t12'])
    c12 = cos(f['t12'])
    s13 = sin(f['t13'])
    c13 = cos(f['t13'])
    s23 = sin(f['t23'])
    c23 = cos(f['t23'])
    # standard phase
    if (s12*s23) == 0 or (c12*c13**2*c23*s13) == 0:
        f['delta'] = 0
    else:
        f['delta'] = -phase((U[0,0].conj()*U[0,2]*U[2,0]*U[2,2].conj()/(c12*c13**2*c23*s13) + c12*c23*s13)/(s12*s23))
    # Majorana phases
    f['delta1']  = phase(exp(1j*f['delta']) * U[0, 2])
    f['delta2']  = phase(U[1, 2])
    f['delta3']  = phase(U[2, 2])
    f['phi1'] = 2*phase(exp(1j*f['delta1']) * U[0, 0].conj())
    f['phi2'] = 2*phase(exp(1j*f['delta1']) * U[0, 1].conj())
    return f
352938
def xstep(self):
        
        """

        ngsit = 0
        gsrrs = np.inf
        YU = self.Y - self.U
        SYU = self.S + YU[..., -1]
        YU[..., -1] = 0.0
        ATYU = self.cnst_AT(YU)
        while gsrrs > self.opt['GSTol'] and ngsit < self.opt['MaxGSIter']:
            self.X = self.GaussSeidelStep(
                SYU, self.X, ATYU, 1.0, self.lcw, 1.0)
            gsrrs = sl.rrs(
                self.cnst_AT(self.cnst_A(self.X)),
                self.cnst_AT(self.cnst_c() - self.cnst_B(self.Y) - self.U)
            )
            ngsit += 1

        self.xs = (ngsit, gsrrs)
660752
public static void fft(boolean forward, long m, double[] x , double[] y)
    {
        int n, i, i1, j, k, i2, l, l1, l2;
        double c1, c2, tx, ty, t1, t2, u1, u2, z;

        /* Calculate the number of points */
        n = 1;
        for (i = 0; i < m; i++)
        {
            n *= 2;
        }

        /* Do the bit reversal */
        i2 = n >> 1;
        j = 0;
        for (i = 0; i < n - 1; i++)
        {
            if (i < j)
            {
                tx = x[i];
                ty = y[i];
                x[i] = x[j];
                y[i] = y[j];
                x[j] = tx;
                y[j] = ty;
            }
            k = i2;
            while (k <= j)
            {
                j -= k;
                k >>= 1;
            }
            j += k;
        }

        /* Compute the FFT */
        c1 = -1.0;
        c2 = 0.0;
        l2 = 1;
        for (l = 0; l < m; l++)
        {
            l1 = l2;
            l2 <<= 1;
            u1 = 1.0;
            u2 = 0.0;
            for (j = 0; j < l1; j++)
            {
                for (i = j; i < n; i += l2)
                {
                    i1 = i + l1;
                    t1 = u1 * x[i1] - u2 * y[i1];
                    t2 = u1 * y[i1] + u2 * x[i1];
                    x[i1] = x[i] - t1;
                    y[i1] = y[i] - t2;
                    x[i] += t1;
                    y[i] += t2;
                }
                z = u1 * c1 - u2 * c2;
                u2 = u1 * c2 + u2 * c1;
                u1 = z;
            }
            c2 = Math.sqrt((1.0 - c1) / 2.0);
            if (forward)
            {
                c2 = -c2;
            }
            c1 = Math.sqrt((1.0 + c1) / 2.0);
        }

        /* Scaling for forward transform */
        if (forward)
        {
            for (i = 0; i < n; i++)
            {
                x[i] /= n;
                y[i] /= n;
            }
        }
    }
1373103
def _getAnalysis(self, axis, analysis, ref=None):
        
        cache = self.cache.get(str(ref))
        if cache is None:
            cache = {'x': {}, 'y': {}, 'z': {}, 'sum': {}}
            self.cache[str(ref)] = cache
        if axis in cache:
            data = self.cache['raw'].get(axis, None)
            cachedAxis = cache.get(axis)
            if cachedAxis.get(analysis) is None:
                if axis == 'sum':
                    if self._canSum(analysis):
                        fx, Pxx = self._getAnalysis('x', analysis)
                        fy, Pxy = self._getAnalysis('y', analysis)
                        fz, Pxz = self._getAnalysis('z', analysis)
                        # calculate the sum of the squares with an additional weighting for x and y
                        Psum = (((Pxx * 2.2) ** 2) + ((Pxy * 2.4) ** 2) + (Pxz ** 2)) ** 0.5
                        if ref is not None:
                            Psum = librosa.amplitude_to_db(Psum, ref)
                        cachedAxis[analysis] = (fx, Psum)
                    else:
                        return None
                else:
                    cachedAxis[analysis] = getattr(data.highPass(), analysis)(ref=ref)
            return cachedAxis[analysis]
        else:
            return None
762529
def calculate_wave(i, lst_bin, wpm=WPM, frequency=FREQUENCY, framerate=FRAMERATE, amplitude=AMPLITUDE, seconds_per_dot=SECONDS_PER_DOT):
    
    bit = morse_bin(i=i, lst_bin=lst_bin, wpm=wpm, framerate=framerate,
                    default_value=0.0, seconds_per_dot=seconds_per_dot)
    sine = sine_wave(i=i, frequency=frequency, framerate=framerate, amplitude=amplitude)
    return bit * sine
844847
def hilbertrot(n, x, y, rx, ry):
    
    if ry == 0:
        if rx == 1:
            x = n - 1 - x
            y = n - 1 - y
        return y, x
    return x, y
1076179
def main(self, x, y, phase):
        
        self.initial_step(phase, x, y)

        # pipelined CORDIC
        for i in range(self.ITERATIONS - 1):
            if self.MODE == CordicMode.ROTATION:
                direction = self.phase[i] > 0
            elif self.MODE == CordicMode.VECTORING:
                direction = self.y[i] < 0

            if direction:
                self.x[i + 1] = self.x[i] - (self.y[i] >> i)
                self.y[i + 1] = self.y[i] + (self.x[i] >> i)
                self.phase[i + 1] = self.phase[i] - self.PHASE_LUT[i]
            else:
                self.x[i + 1] = self.x[i] + (self.y[i] >> i)
                self.y[i + 1] = self.y[i] - (self.x[i] >> i)
                self.phase[i + 1] = self.phase[i] + self.PHASE_LUT[i]

        return self.x[-1], self.y[-1], self.phase[-1]
335525
function hilbert(x, y) {
    let a = x ^ y;
    let b = 0xFFFF ^ a;
    let c = 0xFFFF ^ (x | y);
    let d = x & (y ^ 0xFFFF);

    let A = a | (b >> 1);
    let B = (a >> 1) ^ a;
    let C = ((c >> 1) ^ (b & (d >> 1))) ^ c;
    let D = ((a & (c >> 1)) ^ (d >> 1)) ^ d;

    a = A; b = B; c = C; d = D;
    A = ((a & (a >> 2)) ^ (b & (b >> 2)));
    B = ((a & (b >> 2)) ^ (b & ((a ^ b) >> 2)));
    C ^= ((a & (c >> 2)) ^ (b & (d >> 2)));
    D ^= ((b & (c >> 2)) ^ ((a ^ b) & (d >> 2)));

    a = A; b = B; c = C; d = D;
    A = ((a & (a >> 4)) ^ (b & (b >> 4)));
    B = ((a & (b >> 4)) ^ (b & ((a ^ b) >> 4)));
    C ^= ((a & (c >> 4)) ^ (b & (d >> 4)));
    D ^= ((b & (c >> 4)) ^ ((a ^ b) & (d >> 4)));

    a = A; b = B; c = C; d = D;
    C ^= ((a & (c >> 8)) ^ (b & (d >> 8)));
    D ^= ((b & (c >> 8)) ^ ((a ^ b) & (d >> 8)));

    a = C ^ (C >> 1);
    b = D ^ (D >> 1);

    let i0 = x ^ y;
    let i1 = b | (0xFFFF ^ (i0 | a));

    i0 = (i0 | (i0 << 8)) & 0x00FF00FF;
    i0 = (i0 | (i0 << 4)) & 0x0F0F0F0F;
    i0 = (i0 | (i0 << 2)) & 0x33333333;
    i0 = (i0 | (i0 << 1)) & 0x55555555;

    i1 = (i1 | (i1 << 8)) & 0x00FF00FF;
    i1 = (i1 | (i1 << 4)) & 0x0F0F0F0F;
    i1 = (i1 | (i1 << 2)) & 0x33333333;
    i1 = (i1 | (i1 << 1)) & 0x55555555;

    return ((i1 << 1) | i0) >>> 0;
}
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)
927223
function(signal1, signal2, options){
    options = options || {};
    options = _.defaults(options, defaultOptions);
    
    // ensure signal2 is not longer than signal2
    if(signal1.length < signal2.length){
      var tmp = signal1;
      signal1 = signal2;
      signal2 = tmp;
    }
    
    if(options.method == "cauchy"){
      return cauchy(signal1,signal2,options);
    } else {
      return specConvolution(signal1, signal2,options);
    }
  }