Python numpy.log10() Examples

The following are 30 code examples of numpy.log10(). You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may also want to check out all available functions/classes of the module numpy , or try the search function .
Example #1
Source File: freqresp_test.py    From python-control with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def test_mimo(self):
      # MIMO
      B = np.matrix('1,0;0,1')
      D = np.matrix('0,0')
      sysMIMO = ss(self.A,B,self.C,D)

      frqMIMO = sysMIMO.freqresp(self.omega)
      tfMIMO = tf(sysMIMO)

      #bode(sysMIMO) # - should throw not implemented exception
      #bode(tfMIMO) # - should throw not implemented exception

      #plt.figure(3)
      #plt.semilogx(self.omega,20*np.log10(np.squeeze(frq[0])))

      #plt.figure(4)
      #bode(sysMIMO,self.omega) 
Example #2
Source File: GarrettCompleteness.py    From EXOSIMS with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def maxdmag(self, s):
        """Calculates the maximum value of dMag for projected separation
        
        Args:
            s (float):
                Projected separation (AU)
        
        Returns:
            maxdmag (float):
                Maximum planet delta magnitude
        
        """
        
        if s == 0.0:
            maxdmag = self.cdmax - 2.5*np.log10(self.Phi(np.pi))
        elif s < self.rmax:
            maxdmag = self.cdmax - 2.5*np.log10(np.abs(self.Phi(np.pi-np.arcsin(s/self.rmax))))
        else:
            maxdmag = self.cdmax - 2.5*np.log10(self.Phi(np.pi/2.0))

        return maxdmag 
Example #3
Source File: Stark.py    From EXOSIMS with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def calclogf(self):
        """
        # wavelength dependence, from Table 19 in Leinert et al 1998
        # interpolated w/ a quadratic in log-log space
        Returns:
            interpolant (object):
                a 1D quadratic interpolant of intensity vs wavelength

        """
        self.zodi_lam = np.array([0.2, 0.3, 0.4, 0.5, 0.7, 0.9, 1.0, 1.2, 2.2, 3.5,
                4.8, 12, 25, 60, 100, 140]) # um
        self.zodi_Blam = np.array([2.5e-8, 5.3e-7, 2.2e-6, 2.6e-6, 2.0e-6, 1.3e-6,
                1.2e-6, 8.1e-7, 1.7e-7, 5.2e-8, 1.2e-7, 7.5e-7, 3.2e-7, 1.8e-8,
                3.2e-9, 6.9e-10]) # W/m2/sr/um
        x = np.log10(self.zodi_lam)
        y = np.log10(self.zodi_Blam)
        return interp1d(x, y, kind='quadratic') 
Example #4
Source File: test_PlanetPopulation.py    From EXOSIMS with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def test_dist_sma_radius(self):
        """
        Test that sma and radius values outside of the range have zero probability
        """
        
        for mod in self.allmods:
            if 'dist_sma_radius' in mod.__dict__:
                with RedirectStreams(stdout=self.dev_null):
                    pp = mod(**self.spec)
                
                a = np.logspace(np.log10(pp.arange[0].value/10.),np.log10(pp.arange[1].value*100),100)
                Rp = np.logspace(np.log10(pp.Rprange[0].value/10.),np.log10(pp.Rprange[1].value*100),100)
                
                aa, RR = np.meshgrid(a,Rp)
                
                fr = pp.dist_sma_radius(aa,RR)
                self.assertTrue(np.all(fr[aa < pp.arange[0].value] == 0),'dist_sma_radius low bound failed on sma for %s'%mod.__name__)
                self.assertTrue(np.all(fr[aa > pp.arange[1].value] == 0),'dist_sma_radius high bound failed on sma for %s'%mod.__name__)
                self.assertTrue(np.all(fr[RR < pp.Rprange[0].value] == 0),'dist_sma_radius low bound failed on radius for %s'%mod.__name__)
                self.assertTrue(np.all(fr[RR > pp.Rprange[1].value] == 0),'dist_sma_radius high bound failed on radius for %s'%mod.__name__)
                self.assertTrue(np.all(fr[(aa > pp.arange[0].value) & (aa < pp.arange[1].value) & (RR > pp.Rprange[0].value) & (RR < pp.Rprange[1].value)] > 0),'dist_sma_radius is improper pdf for %s'%mod.__name__) 
Example #5
Source File: test_PlanetPopulation.py    From EXOSIMS with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def test_dist_mass(self):
        """
        Test that masses outside of the range have zero probability

        """

        for mod in self.allmods:
            if 'dist_mass' in mod.__dict__:
                with RedirectStreams(stdout=self.dev_null):
                    pp = mod(**self.spec)

                Mp = np.logspace(np.log10(pp.Mprange[0].value/10.),np.log10(pp.Mprange[1].value*100.),100) 

                fr = pp.dist_mass(Mp)
                self.assertTrue(np.all(fr[Mp < pp.Mprange[0].value] == 0),'dist_mass high bound failed for %s'%mod.__name__)
                self.assertTrue(np.all(fr[Mp > pp.Mprange[1].value] == 0),'dist_mass high bound failed for %s'%mod.__name__)
                self.assertTrue(np.all(fr[(Mp >= pp.Mprange[0].value) & (Mp <= pp.Mprange[1].value)] > 0)) 
Example #6
Source File: test_PlanetPopulation.py    From EXOSIMS with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def test_dist_radius(self):
        """
        Test that radii outside of the range have zero probability

        """
        for mod in self.allmods:
            if 'dist_radius' in mod.__dict__:
                with RedirectStreams(stdout=self.dev_null):
                    pp = mod(**self.spec)

                Rp = np.logspace(np.log10(pp.Rprange[0].to('earthRad').value/10.),np.log10(pp.Rprange[1].to('earthRad').value*100.),100) 

                fr = pp.dist_radius(Rp)
                self.assertTrue(np.all(fr[Rp < pp.Rprange[0].to('earthRad').value] == 0),'dist_radius high bound failed for %s'%mod.__name__)
                self.assertTrue(np.all(fr[Rp > pp.Rprange[1].to('earthRad').value] == 0),'dist_radius high bound failed for %s'%mod.__name__)
                self.assertTrue(np.all(fr[(Rp >= pp.Rprange[0].to('earthRad').value) & (Rp <= pp.Rprange[1].to('earthRad').value)] > 0),'dist_radius generates zero probabilities within range for %s'%mod.__name__) 
Example #7
Source File: test_PlanetPopulation.py    From EXOSIMS with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def test_dist_sma(self):
        """
        Test that smas outside of the range have zero probability

        """

        for mod in self.allmods:
            if 'dist_sma' in mod.__dict__:
                with RedirectStreams(stdout=self.dev_null):
                    pp = mod(**self.spec)

                a = np.logspace(np.log10(pp.arange[0].to('AU').value/10.),np.log10(pp.arange[1].to('AU').value*10.),100)

                fa = pp.dist_sma(a)
                self.assertTrue(np.all(fa[a < pp.arange[0].to('AU').value] == 0),'dist_sma high bound failed for %s'%mod.__name__)
                self.assertTrue(np.all(fa[a > pp.arange[1].to('AU').value] == 0),'dist_sma low bound failed for %s'%mod.__name__)
                self.assertTrue(np.all(fa[(a >= pp.arange[0].to('AU').value) & (a <= pp.arange[1].to('AU').value)] >= 0.),'dist_sma generates negative densities within range for %s'%mod.__name__) 
Example #8
Source File: variables.py    From PyRadarMet with GNU General Public License v2.0 6 votes vote down vote up
def zdr(z_h, z_v):
    """
    Differential reflectivity [dB].

    From Rinehart (1997), Eqn 10.3 and Seliga and Bringi (1976)

    Parameters
    ----------
    z_h : float or array
        Horizontal reflectivity [mm^6/m^3]
    z_v : float or array
        Vertical reflectivity [mm^6/m^3]

    Notes
    -----
    Ensure that both powers have the same units!

    Alternating horizontally and linearly polarized pulses are averaged.

    Notes
    -----
    If arrays are used, either for one okay, or both must be the same length.
    """
    return 10. * np.log10(np.asarray(z_h) / np.asarray(z_v)) 
Example #9
Source File: system.py    From PyRadarMet with GNU General Public License v2.0 6 votes vote down vote up
def gain_Pratio(p1, p2):
    """
    Antenna gain [dB] via power ratio.

    From Rinehart (1997), Eqn 2.1

    Parameters
    ----------
    p1 : float or array
        Power on the beam axis [W]
    p2 : float or array
        Power from an isotropic antenna [W]

    Notes
    -----
    Ensure that both powers have the same units!
    If arrays are used, either for one okay, or both must be the same length.
    """
    return 10. * np.log10(np.asarray(p1) / np.asarray(p2)) 
Example #10
Source File: deltaMag.py    From EXOSIMS with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def deltaMag(p,Rp,d,Phi):
    """ Calculates delta magnitudes for a set of planets, based on their albedo, 
    radius, and position with respect to host star.
    
    Args:
        p (ndarray):
            Planet albedo
        Rp (astropy Quantity array):
            Planet radius in units of km
        d (astropy Quantity array):
            Planet-star distance in units of AU
        Phi (ndarray):
            Planet phase function
    
    Returns:
        dMag (ndarray):
            Planet delta magnitudes
    
    """
    dMag = -2.5*np.log10(p*(Rp/d).decompose()**2*Phi).value
    
    return dMag 
Example #11
Source File: testCore.py    From pyGSTi with Apache License 2.0 6 votes vote down vote up
def test_LGST_1overSqrtN_dependence(self):
        my_datagen_gateset = self.model.depolarize(op_noise=0.05, spam_noise=0)
        # !!don't depolarize spam or 1/sqrt(N) dependence saturates!!

        nSamplesList = np.array([ 16, 128, 1024, 8192 ])
        diffs = []
        for nSamples in nSamplesList:
            ds = pygsti.construction.generate_fake_data(my_datagen_gateset, self.lgstStrings, nSamples,
                                                        sampleError='binomial', seed=100)
            mdl_lgst = pygsti.do_lgst(ds, self.fiducials, self.fiducials, self.model, svdTruncateTo=4, verbosity=0)
            mdl_lgst_go = pygsti.gaugeopt_to_target(mdl_lgst, my_datagen_gateset, {'spam':1.0, 'gate': 1.0}, checkJac=True)
            diffs.append( my_datagen_gateset.frobeniusdist(mdl_lgst_go) )

        diffs = np.array(diffs, 'd')
        a, b = polyfit(np.log10(nSamplesList), np.log10(diffs), deg=1)
        #print "\n",nSamplesList; print diffs; print a #DEBUG
        self.assertLess( a+0.5, 0.05 ) 
Example #12
Source File: optim.py    From End-to-end-ASR-Pytorch with MIT License 6 votes vote down vote up
def speech_aug_scheduler(step, s_r, s_i, s_f, peak_lr):
    # Starting from 0, ramp-up to set LR and  converge to 0.01*LR, w/ exp. decay
    final_lr_ratio = 0.01
    exp_decay_lambda = -np.log10(final_lr_ratio)/(s_f-s_i) # Approx. w/ 10-based
    cur_step = step+1

    if cur_step<s_r:
        # Ramp-up
        return peak_lr*float(cur_step)/s_r
    elif cur_step<s_i:
        # Hold
        return peak_lr
    elif cur_step<=s_f:
        # Decay
        return peak_lr*np.power(10,-exp_decay_lambda*(cur_step-s_i))
    else:
        # Converge
        return peak_lr*final_lr_ratio 
Example #13
Source File: error_evaluation.py    From padasip with MIT License 6 votes vote down vote up
def logSE(x1, x2=-1):
    """
    10 * log10(e**2)    
    This function accepts two series of data or directly
    one series with error.

    **Args:**

    * `x1` - first data series or error (1d array)

    **Kwargs:**

    * `x2` - second series (1d array) if first series was not error directly,\\
        then this should be the second series

    **Returns:**

    * `e` - logSE of error (1d array) obtained directly from `x1`, \\
        or as a difference of `x1` and `x2`. The values are in dB!

    """
    e = get_valid_error(x1, x2)
    return 10*np.log10(e**2) 
Example #14
Source File: ctrlutil.py    From python-control with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def db2mag(db):
    """Convert a gain in decibels (dB) to a magnitude

    If A is magnitude,

        db = 20 * log10(A)

    Parameters
    ----------
    db : float or ndarray
        input value or array of values, given in decibels

    Returns
    -------
    mag : float or ndarray
        corresponding magnitudes

    """
    return 10. ** (db / 20.) 
Example #15
Source File: ctrlutil.py    From python-control with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def mag2db(mag):
    """Convert a magnitude to decibels (dB)

    If A is magnitude,

        db = 20 * log10(A)

    Parameters
    ----------
    mag : float or ndarray
        input magnitude or array of magnitudes

    Returns
    -------
    db : float or ndarray
        corresponding values in decibels

    """
    return 20. * np.log10(mag) 
Example #16
Source File: create_spectrograms.py    From Spoken-language-identification with MIT License 6 votes vote down vote up
def plotstft(audiopath, binsize=2**10, plotpath=None, colormap="gray", channel=0, name='tmp.png', alpha=1, offset=0):
    samplerate, samples = wav.read(audiopath)
    samples = samples[:, channel]
    s = stft(samples, binsize)

    sshow, freq = logscale_spec(s, factor=1, sr=samplerate, alpha=alpha)
    sshow = sshow[2:, :]
    ims = 20.*np.log10(np.abs(sshow)/10e-6) # amplitude to decibel
    timebins, freqbins = np.shape(ims)
    
    ims = np.transpose(ims)
    # ims = ims[0:256, offset:offset+768] # 0-11khz, ~9s interval
    ims = ims[0:256, :] # 0-11khz, ~10s interval
    #print "ims.shape", ims.shape
    
    image = Image.fromarray(ims) 
    image = image.convert('L')
    image.save(name) 
Example #17
Source File: test_bayestar.py    From dustmaps with GNU General Public License v2.0 6 votes vote down vote up
def _interp_ebv(self, datum, dist):
        """
        Calculate samples of E(B-V) at an arbitrary distance (in kpc) for one
        test coordinate.
        """
        dm = 5. * (np.log10(dist) + 2.)
        idx_ceil = np.searchsorted(datum['DM_bin_edges'], dm)
        if idx_ceil == 0:
            dist_0 = 10.**(datum['DM_bin_edges'][0]/5. - 2.)
            return dist/dist_0 * datum['samples'][:,0]
        elif idx_ceil == len(datum['DM_bin_edges']):
            return datum['samples'][:,-1]
        else:
            dm_ceil = datum['DM_bin_edges'][idx_ceil]
            dm_floor = datum['DM_bin_edges'][idx_ceil-1]
            a = (dm_ceil - dm) / (dm_ceil - dm_floor)
            return (
                (1.-a) * datum['samples'][:,idx_ceil]
                +    a * datum['samples'][:,idx_ceil-1]
            ) 
Example #18
Source File: system_eq.py    From BiblioPixelAnimations with MIT License 6 votes vote down vote up
def get_audio_data(self):
        frames = self.rec.get_frames()
        result = [0] * self.bins
        if len(frames) > 0:
            # keeps only the last frame
            current_frame = frames[-1]
            # plots the time signal
            # self.line_top.set_data(self.time_vect, current_frame)
            # computes and plots the fft signal
            fft_frame = np.fft.rfft(current_frame)
            if self.auto_gain:
                fft_frame /= np.abs(fft_frame).max()
            else:
                fft_frame *= (1 + self.gain) / 5000000.

            fft_frame = np.abs(fft_frame)
            if self.log_scale:
                fft_frame = np.log10(np.add(1, np.multiply(10, fft_frame)))

            result = [min(int(max(i, 0.) * 1023), 1023) for i in fft_frame][0:self.bins]

        return result 
Example #19
Source File: system_eq.py    From BiblioPixelAnimations with MIT License 6 votes vote down vote up
def get_audio_data(self):
        frames = self.rec.get_frames()
        result = [0] * self.bins
        if len(frames) > 0:
            # keeps only the last frame
            current_frame = frames[-1]
            # plots the time signal
            # self.line_top.set_data(self.time_vect, current_frame)
            # computes and plots the fft signal
            fft_frame = np.fft.rfft(current_frame)
            if self.auto_gain:
                fft_frame /= np.abs(fft_frame).max()
            else:
                fft_frame *= (1 + self.gain) / 5000000.

            fft_frame = np.abs(fft_frame)
            if self.log_scale:
                fft_frame = np.log10(np.add(1, np.multiply(10, fft_frame)))

            result = [min(int(max(i, 0.) * 1023), 1023) for i in fft_frame][0:self.bins]

        return result 
Example #20
Source File: test_backend.py    From kapre with MIT License 6 votes vote down vote up
def test_amplitude_to_decibel():
    """test for backend_keras.amplitude_to_decibel"""
    from kapre.backend_keras import amplitude_to_decibel

    x = np.array([[1e-20, 1e-5, 1e-3, 5e-2], [0.3, 1.0, 20.5, 9999]])  # random positive numbers

    amin = 1e-5
    dynamic_range = 80.0

    x_decibel = 10 * np.log10(np.maximum(x, amin))
    x_decibel = x_decibel - np.max(x_decibel, axis=(1,), keepdims=True)
    x_decibel_ref = np.maximum(x_decibel, -1 * dynamic_range)

    x_var = K.variable(x)
    x_decibel_kapre = amplitude_to_decibel(x_var, amin, dynamic_range)

    assert np.allclose(K.eval(x_decibel_kapre), x_decibel_ref, atol=TOL) 
Example #21
Source File: compute-vad.py    From pykaldi with Apache License 2.0 5 votes vote down vote up
def show_plot(
    key, segment_times, sample_freqs, spec, duration, wav_data, vad_feat
):
    """This function plots the vad against the signal and the spectrogram.

    Args:
        segment_times: the time intervals acting as the x axis
        sample_freqs: the frequency bins acting as the y axis
        spec: the spectrogram
        duration: duration of the wave file
        wav_data: the wave data
        vad_feat: VAD features
    """

    import matplotlib.pyplot as plt
    import matplotlib.mlab as mlb

    plt.subplot(3, 1, 1)
    plt.pcolormesh(segment_times, sample_freqs, 10 * np.log10(spec), cmap="jet")
    plt.ylabel("Frequency [Hz]")
    plt.xlabel("Time [sec]")

    plt.subplot(3, 1, 2)
    axes = plt.gca()
    axes.set_xlim([0, duration])
    tmp_axis = np.linspace(0, duration, wav_data.shape[0])
    plt.plot(tmp_axis, wav_data / np.abs(np.max(wav_data)))
    plt.xlabel("Time [sec]")

    plt.subplot(3, 1, 3)
    axes = plt.gca()
    axes.set_xlim([0, duration])
    tmp_axis = np.linspace(0, duration, vad_feat.shape[0])
    plt.plot(tmp_axis, vad_feat)
    plt.xlabel("Time [sec]")

    plt.savefig("plots/" + key, bbox_inches="tight") 
Example #22
Source File: slowreplib.py    From pyGSTi with Apache License 2.0 5 votes vote down vote up
def composed(cls, terms_to_compose, magnitude):
        logmag = _math.log10(magnitude) if magnitude > 0 else -LARGE
        first = terms_to_compose[0]
        coeffrep = first.coeff
        pre_ops = first.pre_ops[:]
        post_ops = first.post_ops[:]
        for t in terms_to_compose[1:]:
            coeffrep = coeffrep.mult(t.coeff)
            pre_ops += t.pre_ops
            post_ops += t.post_ops
        return SVTermRep(coeffrep, magnitude, logmag, first.pre_state, first.post_state,
                         first.pre_effect, first.post_effect, pre_ops, post_ops) 
Example #23
Source File: slowreplib.py    From pyGSTi with Apache License 2.0 5 votes vote down vote up
def set_magnitude(self, mag):
        self.magnitude = mag
        self.logmagnitude = _math.log10(mag) if mag > 0 else -LARGE 
Example #24
Source File: run_audio_attack.py    From Black-Box-Audio with MIT License 5 votes vote down vote up
def db(audio):
    if len(audio.shape) > 1:
        maxx = np.max(np.abs(audio), axis=1)
        return 20 * np.log10(maxx) if np.any(maxx != 0) else np.array([0])
    maxx = np.max(np.abs(audio))
    return 20 * np.log10(maxx) if maxx != 0 else np.array([0]) 
Example #25
Source File: timeplots.py    From NanoPlot with GNU General Public License v3.0 5 votes vote down vote up
def length_over_time(dfs, path, figformat, title, log_length=False, plot_settings={}):
    if log_length:
        time_length = Plot(path=path + "TimeLogLengthViolinPlot." + figformat,
                           title="Violin plot of log read lengths over time")
    else:
        time_length = Plot(path=path + "TimeLengthViolinPlot." + figformat,
                           title="Violin plot of read lengths over time")
    sns.set(style="white", **plot_settings)
    if log_length:
        length_column = "log_lengths"
    else:
        length_column = "lengths"

    if "length_filter" in dfs:  # produced by NanoPlot filtering of too long reads
        temp_dfs = dfs[dfs["length_filter"]]
    else:
        temp_dfs = dfs

    ax = sns.violinplot(x="timebin",
                        y=length_column,
                        data=temp_dfs,
                        inner=None,
                        cut=0,
                        linewidth=0)
    ax.set(xlabel='Interval (hours)',
           ylabel="Read length",
           title=title or time_length.title)
    if log_length:
        ticks = [10**i for i in range(10) if not 10**i > 10 * np.amax(dfs["lengths"])]
        ax.set(yticks=np.log10(ticks),
               yticklabels=ticks)
    plt.xticks(rotation=45, ha='center', fontsize=8)
    time_length.fig = ax.get_figure()
    time_length.save(format=figformat)
    plt.close("all")
    return time_length 
Example #26
Source File: hrv_nonlinear.py    From NeuroKit with MIT License 5 votes vote down vote up
def _hrv_nonlinear_poincare(rri, out):
    """Compute SD1 and SD2
    - Do existing measures of Poincare plot geometry reflect nonlinear features of heart rate
    variability? - Brennan (2001)
    """

    # HRV and hrvanalysis
    rri_n = rri[:-1]
    rri_plus = rri[1:]
    x1 = (rri_n - rri_plus) / np.sqrt(2)  # Eq.7
    x2 = (rri_n + rri_plus) / np.sqrt(2)
    sd1 = np.std(x1, ddof=1)
    sd2 = np.std(x2, ddof=1)

    out["SD1"] = sd1
    out["SD2"] = sd2

    # SD1 / SD2
    out["SD1SD2"] = sd1 / sd2

    # Area of ellipse described by SD1 and SD2
    out["S"] = np.pi * out["SD1"] * out["SD2"]

    # CSI / CVI
    T = 4 * out["SD1"]
    L = 4 * out["SD2"]
    out["CSI"] = L / T
    out["CVI"] = np.log10(L * T)
    out["CSI_Modified"] = L ** 2 / T

    return out 
Example #27
Source File: rsp_rrv.py    From NeuroKit with MIT License 5 votes vote down vote up
def _rsp_rrv_nonlinear(bbi):
    diff_bbi = np.diff(bbi)
    out = {}

    # Poincaré plot
    out["SD1"] = np.sqrt(np.std(diff_bbi, ddof=1) ** 2 * 0.5)
    out["SD2"] = np.sqrt(2 * np.std(bbi, ddof=1) ** 2 - 0.5 * np.std(diff_bbi, ddof=1) ** 2)
    out["SD2SD1"] = out["SD2"] / out["SD1"]

    # CSI / CVI
    #    T = 4 * out["SD1"]
    #    L = 4 * out["SD2"]
    #    out["CSI"] = L / T
    #    out["CVI"] = np.log10(L * T)
    #    out["CSI_Modified"] = L ** 2 / T

    # Entropy
    out["ApEn"] = entropy_approximate(bbi, dimension=2)
    out["SampEn"] = entropy_sample(bbi, dimension=2, r=0.2 * np.std(bbi, ddof=1))

    # DFA
    if len(bbi) / 10 > 16:
        out["DFA_1"] = fractal_dfa(bbi, windows=np.arange(4, 17))
    if len(bbi) > 65:
        out["DFA_2"] = fractal_dfa(bbi, windows=np.arange(16, 65))

    return out


# =============================================================================
# Internals
# ============================================================================= 
Example #28
Source File: augment_data.py    From Spoken-language-identification with MIT License 5 votes vote down vote up
def plotstft(audiopath, binsize=2**10, plotpath=None, colormap="gray", channel=0, name='tmp.png', alpha=1, offset=0):
    samplerate, samples = wav.read(audiopath)
    samples = samples[:, channel]
    s = stft(samples, binsize)

    sshow, freq = logscale_spec(s, factor=1, sr=samplerate, alpha=alpha)
    sshow = sshow[2:, :]
    ims = 20.*np.log10(np.abs(sshow)/10e-6) # amplitude to decibel
    timebins, freqbins = np.shape(ims)
    
    ims = np.transpose(ims)
    ims = ims[0:256, offset:offset+768] # 0-11khz, ~9s interval
    #print "ims.shape", ims.shape
    
    image = Image.fromarray(ims) 
    image = image.convert('L')
    image.save(name) 
Example #29
Source File: geom.py    From pyscf with Apache License 2.0 5 votes vote down vote up
def search_c2x(self, zaxis, n):
        '''C2 axis which is perpendicular to z-axis'''
        decimals = int(-numpy.log10(TOLERANCE)) - 1
        for lst in self.group_atoms_by_distance:
            if len(lst) > 1:
                r0 = self.atoms[lst,1:]
                zcos = numpy.around(numpy.einsum('ij,j->i', r0, zaxis),
                                    decimals=decimals)
                uniq_zcos = numpy.unique(zcos)
                maybe_c2x = []
                for d in uniq_zcos:
                    if d > TOLERANCE:
                        mirrord = abs(zcos+d)<TOLERANCE
                        if mirrord.sum() == (zcos==d).sum():
                            above = r0[zcos==d]
                            below = r0[mirrord]
                            nelem = len(below)
                            maybe_c2x.extend([above[0] + below[i]
                                              for i in range(nelem)])
                    elif abs(d) < TOLERANCE: # plane which crosses the orig
                        r1 = r0[zcos==d][0]
                        maybe_c2x.append(r1)
                        r2 = numpy.dot(rotation_mat(zaxis, numpy.pi*2/n), r1)
                        if abs(r1+r2).sum() > TOLERANCE:
                            maybe_c2x.append(r1+r2)
                        else:
                            maybe_c2x.append(r2-r1)

                if len(maybe_c2x) > 0:
                    idx = norm(maybe_c2x, axis=1) > TOLERANCE
                    maybe_c2x = _normalize(maybe_c2x)[idx]
                    maybe_c2x = _remove_dupvec(maybe_c2x)
                    for c2x in maybe_c2x:
                        if (not parallel_vectors(c2x, zaxis) and
                            self.has_rotation(c2x, 2)):
                            return c2x 
Example #30
Source File: Forecaster.py    From EXOSIMS with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def calc_radius_from_mass(self, Mp):
        """Forecast the Radius distribution given the mass distribution.
        
        Args:
            Mp (astropy Quantity array):
                Planet mass in units of Earth mass
        
        Returns:
            Rp (astropy Quantity array):
                Planet radius in units of Earth radius
        
        """
        
        mass = Mp.to('earthMass').value
        assert np.min(mass) > 3e-4 and np.max(mass) < 3e5, \
                "Mass range out of model expectation. Returning None."
        
        sample_size = len(mass)
        logm = np.log10(mass)
        prob = np.random.random(sample_size)
        logr = np.ones_like(logm)
        hyper_ind = np.random.randint(low=0, high=np.shape(self.all_hyper)[0], 
                size=sample_size)
        hyper = self.all_hyper[hyper_ind,:]
        
        for i in range(sample_size):
            logr[i] = self.piece_linear(hyper[i], logm[i], prob[i])
        
        Rp = 10.**logr*u.earthRad
        
        return Rp

#####################################################################################
# The following functions where adapted from the func.py file from the FORECASTER   #
# GitHub repository at https://github.com/chenjj2/forecaster (Chen & Kippling 2016) #
#####################################################################################