Calculate the central wavelength or the central wavenumber,
    depending on what is input
777175
def get_central_wave(wavl, resp):
    
    return np.trapz(resp * wavl, wavl) / np.trapz(resp, wavl)
777449
def calc_spectrum(signal, rate):
    """"""
    npts = len(signal)
    padto = 1 << (npts - 1).bit_length()
    # print 'length of signal {}, pad to {}'.format(npts, padto)
    npts = padto

    sp = np.fft.rfft(signal, n=padto) / npts
    # print('sp len ', len(sp))
    freq = np.arange((npts / 2) + 1) / (npts / rate)
    # print('freq len ', len(freq))
    return freq, abs(sp)
852266
function powerSpectrum(amplitudes) {
  var N = amplitudes.length;

  return amplitudes.map(function (a) {
    return (a * a) / N;
  });
}
1037287
def _read_wave(self, path):
        
        with contextlib.closing(wave.open(path, 'rb')) as wf:
            num_channels = wf.getnchannels()
            assert num_channels == 1
            sample_width = wf.getsampwidth()
            assert sample_width == 2
            sample_rate = wf.getframerate()
            assert sample_rate in (8000, 16000, 32000)
            pcm_data = wf.readframes(wf.getnframes())
            return pcm_data, sample_rate
974407
def _check(self):
        """"""
        # Run the super method
        super(Spectrum, self)._check()

        err_str = None
        has_data = self._KEYS.DATA in self
        has_wave = self._KEYS.WAVELENGTHS in self
        has_flux = self._KEYS.FLUXES in self
        has_filename = self._KEYS.FILENAME in self

        if not has_data:
            if (not has_wave or not has_flux) and not has_filename:
                err_str = (
                    "If `{}` not given".format(self._KEYS.DATA) +
                    "; `{}` or `{}` needed".format(
                        self._KEYS.WAVELENGTHS, self._KEYS.FLUXES))

        if err_str is not None:
            raise ValueError(err_str)

        return
806685
def wave_units(self, units):
        
        # Make sure it's length units
        if not units.is_equivalent(q.m):
            raise ValueError(units, ": New wavelength units must be a length.")

        # Update the units
        self._wave_units = units

        # Update all the wavelength values
        self._wave = self.wave.to(self.wave_units).round(5)
        self.wave_min = self.wave_min.to(self.wave_units).round(5)
        self.wave_max = self.wave_max.to(self.wave_units).round(5)
        self.wave_eff = self.wave_eff.to(self.wave_units).round(5)
        self.wave_center = self.wave_center.to(self.wave_units).round(5)
        self.wave_mean = self.wave_mean.to(self.wave_units).round(5)
        self.wave_peak = self.wave_peak.to(self.wave_units).round(5)
        self.wave_phot = self.wave_phot.to(self.wave_units).round(5)
        self.wave_pivot = self.wave_pivot.to(self.wave_units).round(5)
        self.width_eff = self.width_eff.to(self.wave_units).round(5)
        self.fwhm = self.fwhm.to(self.wave_units).round(5)
399352
def _process_wave_param(self, pval):
        """"""
        return self._process_generic_param(
            pval, self._internal_wave_unit, equivalencies=u.spectral())
935828
def get_reflectance(self, sun_zenith, sat_zenith, azidiff, bandname, redband=None):
        """"""
        # Get wavelength in nm for band:
        if isinstance(bandname, float):
            LOG.warning('A wavelength is provided instead of band name - ' +
                        'disregard the relative spectral responses and assume ' +
                        'it is the effective wavelength: %f (micro meter)', bandname)
            wvl = bandname * 1000.0
        else:
            wvl = self.get_effective_wavelength(bandname)
            wvl = wvl * 1000.0

        rayl, wvl_coord, azid_coord, satz_sec_coord, sunz_sec_coord = self.get_reflectance_lut()

        # force dask arrays
        compute = False
        if HAVE_DASK and not isinstance(sun_zenith, Array):
            compute = True
            sun_zenith = from_array(sun_zenith, chunks=sun_zenith.shape)
            sat_zenith = from_array(sat_zenith, chunks=sat_zenith.shape)
            azidiff = from_array(azidiff, chunks=azidiff.shape)
            if redband is not None:
                redband = from_array(redband, chunks=redband.shape)

        clip_angle = rad2deg(arccos(1. / sunz_sec_coord.max()))
        sun_zenith = clip(sun_zenith, 0, clip_angle)
        sunzsec = 1. / cos(deg2rad(sun_zenith))
        clip_angle = rad2deg(arccos(1. / satz_sec_coord.max()))
        sat_zenith = clip(sat_zenith, 0, clip_angle)
        satzsec = 1. / cos(deg2rad(sat_zenith))
        shape = sun_zenith.shape

        if not(wvl_coord.min() < wvl < wvl_coord.max()):
            LOG.warning(
                "Effective wavelength for band %s outside 400-800 nm range!",
                str(bandname))
            LOG.info(
                "Set the rayleigh/aerosol reflectance contribution to zero!")
            if HAVE_DASK:
                chunks = sun_zenith.chunks if redband is None else redband.chunks
                res = zeros(shape, chunks=chunks)
                return res.compute() if compute else res
            else:
                return zeros(shape)

        idx = np.searchsorted(wvl_coord, wvl)
        wvl1 = wvl_coord[idx - 1]
        wvl2 = wvl_coord[idx]

        fac = (wvl2 - wvl) / (wvl2 - wvl1)
        raylwvl = fac * rayl[idx - 1, :, :, :] + (1 - fac) * rayl[idx, :, :, :]
        tic = time.time()

        smin = [sunz_sec_coord[0], azid_coord[0], satz_sec_coord[0]]
        smax = [sunz_sec_coord[-1], azid_coord[-1], satz_sec_coord[-1]]
        orders = [
            len(sunz_sec_coord), len(azid_coord), len(satz_sec_coord)]
        f_3d_grid = atleast_2d(raylwvl.ravel())

        if HAVE_DASK and isinstance(smin[0], Array):
            # compute all of these at the same time before passing to the interpolator
            # otherwise they are computed separately
            smin, smax, orders, f_3d_grid = da.compute(smin, smax, orders, f_3d_grid)
        minterp = MultilinearInterpolator(smin, smax, orders)
        minterp.set_values(f_3d_grid)

        if HAVE_DASK:
            ipn = map_blocks(self._do_interp, minterp, sunzsec, azidiff,
                             satzsec, dtype=raylwvl.dtype, chunks=azidiff.chunks)
        else:
            ipn = self._do_interp(minterp, sunzsec, azidiff, satzsec)

        LOG.debug("Time - Interpolation: {0:f}".format(time.time() - tic))

        ipn *= 100
        res = ipn
        if redband is not None:
            res = where(redband < 20., res,
                        (1 - (redband - 20) / 80) * res)

        res = clip(res, 0, 100)
        if compute:
            res = res.compute()
        return res
188612
def wavelength_match(a, b):
        
        if type(a) == (type(b) or
                       isinstance(a, numbers.Number) and
                       isinstance(b, numbers.Number)):
            return a == b
        elif a is None or b is None:
            return False
        elif isinstance(a, (list, tuple)) and len(a) == 3:
            return a[0] <= b <= a[2]
        elif isinstance(b, (list, tuple)) and len(b) == 3:
            return b[0] <= a <= b[2]
        else:
            raise ValueError("Can only compare wavelengths of length 1 or 3")
205660
def convert(self):
        """"""

        from pyspectral.utils import (convert2wavenumber, get_central_wave)
        if self._wavespace == WAVE_LENGTH:
            rsr, info = convert2wavenumber(self.rsr)
            for band in rsr.keys():
                for det in rsr[band].keys():
                    self.rsr[band][det][WAVE_NUMBER] = rsr[
                        band][det][WAVE_NUMBER]
                    self.rsr[band][det]['response'] = rsr[
                        band][det]['response']
                    self.unit = info['unit']
                    self.si_scale = info['si_scale']
            self._wavespace = WAVE_NUMBER
            for band in rsr.keys():
                for det in rsr[band].keys():
                    self.rsr[band][det]['central_wavenumber'] = \
                        get_central_wave(self.rsr[band][det][WAVE_NUMBER], self.rsr[band][det]['response'])
                    del self.rsr[band][det][WAVE_LENGTH]
        else:
            errmsg = "Conversion from {wn} to {wl} not supported yet".format(wn=WAVE_NUMBER, wl=WAVE_LENGTH)
            raise NotImplementedError(errmsg)
592288
def unitResponse(self, band):
        """"""
        #sumfilt(wave,0,band)
        # SUMFILT = Sum [ FILT(I) * WAVE(I) ** NPOW * DWAVE(I) ]
        total = band.trapezoidIntegration(band.wave, band.throughput)
        return 1.0/total
361026
def model(self):
        """"""
        if self.z == 0:
            m = self._model
        else:
            # wavelength
            if self._internal_wave_unit.physical_type == 'length':
                rs = self._redshift_model.inverse
            # frequency or wavenumber
            # NOTE: This will never execute as long as internal wavelength
            #       unit remains Angstrom.
            else:  # pragma: no cover
                rs = self._redshift_model

            if self.z_type == 'wavelength_only':
                m = rs | self._model
            else:  # conserve_flux
                m = rs | self._model | self._redshift_flux_model

        return m
1139667
def waveform_length(X):
    '''  '''
    return np.sum(np.abs(np.diff(X, axis=1)), axis=1)
586688
def LaplaceCentreWeight(self):
        """"""

        sz = [1,] * self.S.ndim
        for ax in self.axes:
            sz[ax] = self.S.shape[ax]
        lcw = 2*len(self.axes)*np.ones(sz, dtype=self.dtype)
        for ax in self.axes:
            lcw[(slice(None),)*ax + ([0, -1],)] -= 1.0
        return lcw
576023
def _getWaveProp(self):
        """"""
        wave = units.Angstrom().Convert(self._wavetable, self.waveunits.name)
        return wave
378130
def pivot_wavelength_ee(bpass):
    
    from scipy.integrate import simps
    return np.sqrt(simps(bpass.resp, bpass.wlen) /
                   simps(bpass.resp / bpass.wlen**2, bpass.wlen))
819243
def Throughput(self):
        
        try:
            throughput = spectrum.TabularSpectralElement()

            product = self._multiplyThroughputs(0)

            throughput._wavetable = product.GetWaveSet()
            throughput._throughputtable = product(throughput._wavetable)
            throughput.waveunits = product.waveunits
            throughput.name='*'.join([str(x) for x in self.components])

##            throughput = throughput.resample(spectrum._default_waveset)

            return throughput

        except IndexError:   # graph table is broken.
            return None