Perform Short-time Fourier transform to get the spectrogram for the given samples
        :return: short-time Fourier transform of the given signal
977512
def stft(self, samples: np.ndarray):
        
        window = self.window_function(self.window_size)
        hop_size = self.hop_size

        if len(samples) < self.window_size:
            samples = np.append(samples, np.zeros(self.window_size - len(samples)))

        num_frames = max(1, (len(samples) - self.window_size) // hop_size + 1)

        # Get frames as numpy view with stride_tricks to save RAM
        # Same as: frames = [padded_samples[i*hop_size:i*hop_size+self.window_size] for i in range(num_frames)]
        shape = (num_frames, self.window_size)
        strides = (hop_size * samples.strides[-1], samples.strides[-1])
        frames = np.lib.stride_tricks.as_strided(samples, shape=shape, strides=strides)

        result = np.fft.fft(frames * window, self.window_size) / np.atleast_1d(self.window_size)
        return result
1329675
def find_frequencies(data, freq=44100, bits=16):
    
    # Fast fourier transform
    n = len(data)
    p = _fft(data)
    uniquePts = numpy.ceil((n + 1) / 2.0)

    # Scale by the length (n) and square the value to get the amplitude
    p = [(abs(x) / float(n)) ** 2 * 2 for x in p[0:uniquePts]]
    p[0] = p[0] / 2
    if n % 2 == 0:
        p[-1] = p[-1] / 2

    # Generate the frequencies and zip with the amplitudes
    s = freq / float(n)
    freqArray = numpy.arange(0, uniquePts * s, s)
    return zip(freqArray, p)
1163746
def rms(signal, fs):
    
    # if a signal contains a some silence, taking the RMS of the whole
    # signal will be calculated as less loud as a signal without a silent
    # period. I don't like this, so I am going to chunk the signals, and
    # take the value of the most intense chunk
    chunk_time = 0.001  # 1 ms chunk
    chunk_samps = int(chunk_time * fs)
    amps = []
    if chunk_samps > 10:
        for i in range(0, len(signal) - chunk_samps, chunk_samps):
            amps.append(np.sqrt(np.mean(pow(signal[i:i + chunk_samps], 2))))
        amps.append(np.sqrt(np.mean(pow(signal[len(signal) - chunk_samps:], 2))))
        return np.amax(amps)
    else:
        # samplerate low, just rms the whole thing
        return np.sqrt(np.mean(pow(signal, 2)))
1120698
def bandWave(self):
        
        if self.binset.startswith('('):
            return self._computeBandwave(self.binset)
        else:
            return self._getBandwaveFomFile(self.binset)
429608
def subsample(self):
        
        # Create the threads for the analysis
        logging.info('Subsampling FASTQ reads')
        for _ in range(self.cpus):
            threads = Thread(target=self.subsamplethreads, args=())
            threads.setDaemon(True)
            threads.start()
        with progressbar(self.runmetadata.samples) as bar:
            for sample in bar:
                if sample.general.bestassemblyfile != 'NA':
                    # Set the name of the subsampled FASTQ file
                    sample[self.analysistype].subsampledfastq = \
                        os.path.splitext(sample[self.analysistype].baitedfastq)[0] + '_subsampled.fastq'
                    # Set the system call
                    sample[self.analysistype].seqtkcall = 'reformat.sh in={} out={} samplereadstarget=1000'\
                        .format(sample[self.analysistype].baitedfastq,
                                sample[self.analysistype].subsampledfastq)
                    # Add the sample to the queue
                    self.samplequeue.put(sample)
        self.samplequeue.join()
338540
function DFT(bufferSize, sampleRate) {
  FourierTransform.call(this, bufferSize, sampleRate);

  var N = bufferSize/2 * bufferSize;
  var TWO_PI = 2 * Math.PI;

  this.sinTable = new Float32Array(N);
  this.cosTable = new Float32Array(N);

  for (var i = 0; i < N; i++) {
    this.sinTable[i] = Math.sin(i * TWO_PI / bufferSize);
    this.cosTable[i] = Math.cos(i * TWO_PI / bufferSize);
  }
}
113099
def eigenvalues_samples(self):
         """
        res = np.empty((self.nsamples, self.nstates), dtype=config.dtype)
        for i in range(self.nsamples):
            res[i, :] = self._sampled_hmms[i].eigenvalues
        return res
1138010
def welch(timeseries, segmentlength, noverlap=None, **kwargs):
    
    # calculate PSD
    freqs, psd_ = scipy.signal.welch(
        timeseries.value,
        noverlap=noverlap,
        fs=timeseries.sample_rate.decompose().value,
        nperseg=segmentlength,
        **kwargs
    )
    # generate FrequencySeries and return
    unit = scale_timeseries_unit(
        timeseries.unit,
        kwargs.get('scaling', 'density'),
    )
    return FrequencySeries(
        psd_,
        unit=unit,
        frequencies=freqs,
        name=timeseries.name,
        epoch=timeseries.epoch,
        channel=timeseries.channel,
    )
840578
def A(self):
        
        return fft(np.dstack([np.eye(self.m), -self.b]),
                   self.nfft * 2 - 1)[:, :, :self.nfft]
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']))
1310686
def process(self, frames, eod, spec_range=120.0):
        

        samples = frames[:, 0]
        nsamples = len(frames[:, 0])
        if nsamples != self.blocksize:
            self.window = self.window_function(nsamples)
        samples *= self.window

        while nsamples > self.fft_size:
            self.fft_size = 2 * self.fft_size

        zeros_p = numpy.zeros(self.fft_size / 2 - int(nsamples / 2))
        if nsamples % 2:
            zeros_n = numpy.zeros(self.fft_size / 2 - int(nsamples / 2) - 1)
        else:
            zeros_n = numpy.zeros(self.fft_size / 2 - int(nsamples / 2))
        samples = numpy.concatenate((zeros_p, samples, zeros_n), axis=0)

        fft = numpy.fft.fft(samples)
        # normalized abs(FFT) between 0 and 1
        spectrum = numpy.abs(fft[:fft.shape[0] / 2 + 1]) / float(nsamples)
        length = numpy.float64(spectrum.shape[0])

        # scale the db spectrum from [- spec_range db ... 0 db] > [0..1]
        db_spectrum = ((20 * (numpy.log10(spectrum + 1e-30)))
                       .clip(-spec_range, 0.0) + spec_range) / spec_range
        energy = spectrum.sum()
        spectral_centroid = 0

        if energy > 1e-20:
            # calculate the spectral centroid
            if self.spectrum_range is None:
                self.spectrum_range = numpy.arange(length)
            spectral_centroid = (spectrum * self.spectrum_range).sum() / \
                (energy * (length - 1)) * \
                self.samplerate * 0.5
            # clip > log10 > scale between 0 and 1
            spectral_centroid = (math.log10(self.clip(spectral_centroid,
                                                      self.lower,
                                                      self.higher)) -
                                 self.lower_log) / (self.higher_log -
                                                    self.lower_log)

        return (spectral_centroid, db_spectrum)
1140207
def _sum_by_samples(seqs_freq, samples_order):
    
    n = len(seqs_freq[seqs_freq.keys()[0]].freq.keys())
    y = np.array([0] * n)
    for s in seqs_freq:
        x = seqs_freq[s].freq
        exp = [seqs_freq[s].freq[sam] for sam in samples_order]
        y = list(np.array(exp) + y)
    return y
1218555
def plot_spectrogram(self, node_idx=None):
        """
        from pygsp.plotting import _plot_spectrogram
        _plot_spectrogram(self, node_idx=node_idx)
430001
def get_spectrogram(self):
        """"""
        args = self.args

        asd = self.timeseries[0].asd().value
        if (asd.min() == 0):
            self.log(0, 'Input data has a zero in ASD. '
                     'Q-transform not possible.')
            self.got_error = True
            qtrans = None
        else:
            gps = self.qxfrm_args['gps']
            outseg = Segment(gps, gps).protract(args.plot[self.plot_num])

            # This section tries to optimize the amount of data that is
            # processed and the time resolution needed to create a good
            # image. NB:For each time span specified
            # NB: the timeseries h enough data for the longest plot
            inseg = outseg.protract(4) & self.timeseries[0].span
            proc_ts = self.timeseries[0].crop(*inseg)

            #  time resolution is calculated to provide about 4 times
            # the number of output pixels for interpolation
            tres = float(outseg.end - outseg.start) / 4 / self.args.nx
            self.qxfrm_args['tres'] = tres
            self.qxfrm_args['search'] = int(len(proc_ts) * proc_ts.dt.value)

            self.log(3, 'Q-transform arguments:')
            self.log(3, '{0:>15s} = {1}'.format('outseg', outseg))
            for key in sorted(self.qxfrm_args):
                self.log(3, '{0:>15s} = {1}'.format(key, self.qxfrm_args[key]))

            qtrans = proc_ts.q_transform(outseg=outseg, **self.qxfrm_args)

            if args.ymin is None:  # set before Spectrogram.make_plot
                args.ymin = qtrans.yspan[0]

        return qtrans
394133
def synthesize(self, duration):
        
        sr = self.samplerate.samples_per_second
        seconds = duration / Seconds(1)
        samples = np.random.uniform(low=-1., high=1., size=int(sr * seconds))
        return AudioSamples(samples, self.samplerate)
763641
def get_samples(self, sample_count):
        
        if self.amplitude.value <= 0:
            return None
        # Build samples by rolling the period cache through the buffer
        rolled_array = numpy.roll(self.wave_cache,
                                  -1 * self.last_played_sample)
        # Append remaining partial period
        full_count, remainder = divmod(sample_count, self.cache_length)
        final_subarray = rolled_array[:int(remainder)]
        return_array = numpy.concatenate((numpy.tile(rolled_array, full_count),
                                          final_subarray))
        # Keep track of where we left off to prevent popping between chunks
        self.last_played_sample = int(((self.last_played_sample + remainder) %
                                       self.cache_length))
        # Multiply output by amplitude
        return return_array * (self.amplitude.value *
                               self.amplitude_multiplier)
1039051
def plot_welch_perdiogram(x, fs, nperseg):
    '''
    '''
    import scipy.signal
    import numpy

    # Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by
    # 0.001V**2/Hz of white noise sampled at 10 kHz.
    N = len(x)
    time = numpy.arange(N) / fs

    # Compute and plot the power spectral density.
    f, Pxx_den = scipy.signal.welch(x, fs, nperseg=nperseg)

    plt.semilogy(f, Pxx_den)
    plt.ylim([0.5e-3, 1])
    plt.xlabel('frequency [Hz]')
    plt.ylabel('PSD [V**2/Hz]')
    plt.show()

    # If we average the last half of the spectral density, to exclude the peak,
    # we can recover the noise power on the signal.
    numpy.mean(Pxx_den[256:])  # 0.0009924865443739191

    # compute power spectrum
    f, Pxx_spec = scipy.signal.welch(x, fs, 'flattop', 1024,
                                     scaling='spectrum')

    plt.figure()
    plt.semilogy(f, numpy.sqrt(Pxx_spec))
    plt.xlabel('frequency [Hz]')
    plt.ylabel('Linear spectrum [V RMS]')
    plt.show()

    return None