``B. C. J. Moore and B. R. Glasberg, "Suggested formulae for calculating
    auditory filter bandwidths and excitation patterns". J. Acoust. Soc.
    Am., 74, 1983, pp. 750-753.``
341610
def erb(freq, Hz=None):
  
  if Hz is None:
    if freq < 7: # Perhaps user tried something up to 2 * pi
      raise ValueError("Frequency out of range.")
    Hz = 1
  fHz = freq / Hz
  result = 6.23e-6 * fHz ** 2 + 93.39e-3 * fHz + 28.52
  return result * Hz
317964
def calc_uncorrelated_signal_power(params):
    
    ''''''


    for i, data_type in enumerate(['LFP','CSD']):
        if i % SIZE == RANK:
    
            # Determine size of PSD matrix
    
            f = h5py.File(os.path.join(params.savefolder, data_type + 'sum.h5'),'r')
            data = f['data'].value
            srate = f['srate'].value
            if ana_params.mlab:
                Psum, freqs = plt.mlab.psd(data[0], NFFT=ana_params.NFFT, Fs=srate,
                                    noverlap=ana_params.noverlap, window=ana_params.window)
            else:
                [freqs, Psum] = hlp.powerspec([data[0]], tbin= 1./srate*1000.,
                                       Df=ana_params.Df, pointProcess=False)
            f.close()
            P = np.zeros((data.shape[0],Psum.shape[0]))
               
            for y in params.y:
    
                print 'processing ', y
                # Load data
                f = h5py.File(os.path.join(params.populations_path, '%s_%ss' %
                                           (y,data_type) + '.h5'),'r')
                data_y = f['data'].value[:,:, ana_params.transient:]
                # subtract mean
                for j in range(len(data_y)):
                    data_yT = data_y[j].T - data_y[j].mean(axis=1)
                    data_y[j] = data_yT.T
                srate = f['srate'].value
                tvec = np.arange(data_y.shape[2]) * 1000. / srate
                f.close()
    
            
                for j in range(len(data_y)): # loop over cells
                    if ana_params.mlab:
                        for ch in range(len(params.electrodeParams['z'])): # loop over channels
                            P_j_ch, freqs = plt.mlab.psd(data_y[j,ch],
                                                         NFFT=ana_params.NFFT, Fs=srate,
                                                         noverlap=ana_params.noverlap,
                                                         window=ana_params.window)
                                       
                            P[ch] += P_j_ch
                    else:
                        [freqs, P_j] = hlp.powerspec(data_y[j], tbin= 1./srate*1000.,
                                                     Df=ana_params.Df, pointProcess=False)
                        mask = np.where(freqs >= 0.)
                        freqs = freqs[mask]
                        P_j = P_j[:,mask][:,0,:]
                        P_j = P_j/tvec[tvec >= ana_params.transient].size**2
                        
                        P += P_j
    
            #rescale PSD as they may be computed from a fraction of single cell LFPs
            P /= params.recordSingleContribFrac
            
            # Save data
            f = h5py.File(os.path.join(params.savefolder, ana_params.analysis_folder,
                                       data_type +  ana_params.fname_psd_uncorr),'w')
            f['freqs']=freqs
            f['psd']=P
            f['transient']=ana_params.transient
            f['mlab']=ana_params.mlab
            f['NFFT']=ana_params.NFFT
            f['noverlap']=ana_params.noverlap
            f['window']=str(ana_params.window)
            f['Df']=str(ana_params.Df)        
            f.close()

    return
142622
def stSpectralFlux(X, X_prev):
    
    # compute the spectral flux as the sum of square distances:
    sumX = numpy.sum(X + eps)
    sumPrevX = numpy.sum(X_prev + eps)
    F = numpy.sum((X / sumX - X_prev/sumPrevX) ** 2)

    return F
630732
def filter(self, s, method='exact', order=None):
        """"""
        return super(Gabor, self).filter(s, method='exact')
591202
def lowpass_cheby2_coeffs(freq, sr, maxorder=12):
    # type: (float, int, int) -> Tuple[np.ndarray, np.ndarray, float]
    
    
    nyquist = sr * 0.5
    # rp - maximum ripple of passband, rs - attenuation of stopband
    rp, rs, order = 1, 96, 1e99
    ws = freq / nyquist  # stop band frequency
    wp = ws              # pass band frequency
    # raise for some bad scenarios
    if ws > 1:
        ws = 1.0
        msg = "Selected corner frequency is above Nyquist. " + \
              "Setting Nyquist as high corner."
        warnings.warn(msg)
    while True:
        if order <= maxorder:
            break
        wp = wp * 0.99
        order, wn = signal.cheb2ord(wp, ws, rp, rs, analog=0)
    b, a = signal.cheby2(order, rs, wn, btype='low', analog=0, output='ba')
    return b, a, wp*nyquist
698492
function _getDurationsFromProperties(dimensions, coefficients, speedOfSound) {
  let durations = new Float32Array(Utils.NUMBER_REVERB_FREQUENCY_BANDS);

  // Sanitize inputs.
  dimensions = _sanitizeDimensions(dimensions);
  coefficients = _sanitizeCoefficients(coefficients);
  if (speedOfSound == undefined) {
    speedOfSound = Utils.DEFAULT_SPEED_OF_SOUND;
  }

  // Acoustic constant.
  let k = Utils.TWENTY_FOUR_LOG10 / speedOfSound;

  // Compute volume, skip if room is not present.
  let volume = dimensions.width * dimensions.height * dimensions.depth;
  if (volume < Utils.ROOM_MIN_VOLUME) {
    return durations;
  }

  // Room surface area.
  let leftRightArea = dimensions.width * dimensions.height;
  let floorCeilingArea = dimensions.width * dimensions.depth;
  let frontBackArea = dimensions.depth * dimensions.height;
  let totalArea = 2 * (leftRightArea + floorCeilingArea + frontBackArea);
  for (let i = 0; i < Utils.NUMBER_REVERB_FREQUENCY_BANDS; i++) {
    // Effective absorptive area.
    let absorbtionArea =
      (coefficients.left[i] + coefficients.right[i]) * leftRightArea +
      (coefficients.down[i] + coefficients.up[i]) * floorCeilingArea +
      (coefficients.front[i] + coefficients.back[i]) * frontBackArea;
    let meanAbsorbtionArea = absorbtionArea / totalArea;

    // Compute reverberation using Eyring equation [1].
    // [1] Beranek, Leo L. "Analysis of Sabine and Eyring equations and their
    //     application to concert hall audience and chair absorption." The
    //     Journal of the Acoustical Society of America, Vol. 120, No. 3.
    //     (2006), pp. 1399-1399.
    durations[i] = Utils.ROOM_EYRING_CORRECTION_COEFFICIENT * k * volume /
      (-totalArea * Math.log(1 - meanAbsorbtionArea) + 4 *
      Utils.ROOM_AIR_ABSORPTION_COEFFICIENTS[i] * volume);
  }
  return durations;
}
557750
def sph_harmonic_tp(nrows, ncols, n, m):
    
    nuvec = pysphi.ynunm(n, m, n + 1)
    num = nuvec[1:] * (-1) ** m
    num = num[::-1]
    yvec = np.concatenate([num, nuvec])
    
    theta = np.linspace(0.0, np.pi, nrows)
    phi = np.linspace(0.0, 2.0 * np.pi - 2.0 * np.pi / ncols)
    nu = np.array(range(-n, n + 1), dtype=np.complex128)
    
    out = np.zeros([nrows, ncols], dtype=np.complex128)

    for nn in xrange(0, nrows):
        exp_vec = np.exp(1j * nu * theta[nn])
        for mm in xrange(0, ncols):          
            vl = np.sum(yvec * exp_vec)
            out[nn, mm] = 1j ** m * np.exp(1j * m * phi[mm]) * vl

    return out
445579
def Q(self):
        
        return np.array(list(self.center_frequencies)) \
               / np.array(list(self.bandwidths))
1056661
def calc_signal_power(params):

    '''
    

    '''

    for i, data_type in enumerate(['CSD','LFP','CSD_10_0', 'LFP_10_0']):
        if i % SIZE == RANK:

            # Load data
            if data_type in ['CSD','LFP']:
                fname=os.path.join(params.savefolder, data_type+'sum.h5')
            else:
                fname=os.path.join(params.populations_path, 'subsamples',
                                   str.split(data_type,'_')[0] + 'sum_' +
                                   str.split(data_type,'_')[1] + '_' +
                                   str.split(data_type,'_')[2] + '.h5')
            #open file
            f = h5py.File(fname)
            data = f['data'].value
            srate = f['srate'].value 
            tvec = np.arange(data.shape[1]) * 1000. / srate
        
            # slice
            slica = (tvec >= ana_params.transient)
            data = data[:,slica]
    
            # subtract mean
            dataT = data.T - data.mean(axis=1)
            data = dataT.T
            f.close()
    
            #extract PSD
            PSD=[]
            for i in np.arange(len(params.electrodeParams['z'])):
                if ana_params.mlab:
                    Pxx, freqs=plt.mlab.psd(data[i], NFFT=ana_params.NFFT,
                                        Fs=srate, noverlap=ana_params.noverlap,
                                        window=ana_params.window)
                else:
                    [freqs, Pxx] = hlp.powerspec([data[i]], tbin= 1.,
                                             Df=ana_params.Df, pointProcess=False)
                    mask = np.where(freqs >= 0.)
                    freqs = freqs[mask]
                    Pxx = Pxx.flatten()
                    Pxx = Pxx[mask]
                    Pxx = Pxx/tvec[tvec >= ana_params.transient].size**2
                PSD +=[Pxx.flatten()]
                
            PSD=np.array(PSD)
    
            # Save data
            f = h5py.File(os.path.join(params.savefolder, ana_params.analysis_folder,
                                       data_type + ana_params.fname_psd),'w')
            f['freqs']=freqs
            f['psd']=PSD
            f['transient']=ana_params.transient
            f['mlab']=ana_params.mlab
            f['NFFT']=ana_params.NFFT
            f['noverlap']=ana_params.noverlap
            f['window']=str(ana_params.window)
            f['Df']=str(ana_params.Df)
            f.close()

    return
403428
def bandpass_filter(rate=None, low=None, high=None, order=None):
    """"""
    assert low < high
    assert order >= 1
    return signal.butter(order,
                         (low / (rate / 2.), high / (rate / 2.)),
                         'pass')
551211
def _compute_soil_amplification(self, C, vs30, pga_bc, mean):
        
        # convert from base e (as defined in BA2008) to base 10 (as used in
        # AB2006)
        sal = np.log10(np.exp(self._get_site_amplification_linear(vs30, C)))
        sanl = np.log10(np.exp(
            self._get_site_amplification_non_linear(vs30, pga_bc, C)))

        idxs = vs30 < 2000.0
        mean[idxs] = mean[idxs] + sal[idxs] + sanl[idxs]
1366201
def _twofilter_smoothing_ON2(self, t, ti, info, phi, lwinfo):
        
        sp, sw = 0., 0.
        upb = lwinfo.max() + self.wgt[t].lw.max()
        if hasattr(self.model, 'upper_bound_trans'):
            upb += self.model.upper_bound_trans(t + 1)
        # Loop over n, to avoid having in memory a NxN matrix
        for n in range(self.N):
            omegan = np.exp(lwinfo + self.wgt[t].lw[n] - upb
                            + self.model.logpt(t + 1, self.X[t][n],
                                               info.hist.X[ti]))
            sp += np.sum(omegan * phi(self.X[t][n], info.hist.X[ti]))
            sw += np.sum(omegan)
        return sp / sw
574369
def lm_deltat(freqs, damping_times, modes):
    

    dt = {}
    for lmn in modes:
        l, m, nmodes = int(lmn[0]), int(lmn[1]), int(lmn[2])
        for n in range(nmodes):
            dt['%d%d%d' %(l,m,n)] = 1. / qnm_freq_decay(freqs['%d%d%d' %(l,m,n)],
                                       damping_times['%d%d%d' %(l,m,n)], 1./1000)

    delta_t = min(dt.values())
    if delta_t < min_dt:
        delta_t = min_dt

    return delta_t
399701
def getSoundFileDuration(fn):
    '''
    
    '''
    audiofile = wave.open(fn, "r")
    
    params = audiofile.getparams()
    framerate = params[2]
    nframes = params[3]
    
    duration = float(nframes) / framerate
    return duration
471088
def edf_greenhall_simple(alpha, d, m, S, F, N):
    
    L = m/F+m*d # length of filter applied to phase samples
    M = 1 + np.floor(S*(N-L) / m)
    J = min(M, (d+1)*S)
    inv_edf = (1.0/(pow(greenhall_sz(0, F, alpha, d), 2)*M))* \
               greenhall_BasicSum(J, M, S, F, alpha, d)
    return 1.0/inv_edf
1019298
private void mixChannelIntoBuffers(final int[] leftBuffer, final int[] rightBuffer, final int startIndex, final int endIndex, final ChannelMemory actMemo)
	{
		for (int i=startIndex; i<endIndex; i++)
		{
			// Retrieve the Sampledata for this point (interpolated, if necessary)
			long sample = actMemo.currentSample.getInterpolatedSample(doISP, actMemo.currentSamplePos, actMemo.currentTuningPos);

			// Resonance Filters
			if (actMemo.filterOn) sample = doResonance(actMemo, sample);

			// Volume Ramping (deClick) Left!
			long volL = actMemo.actRampVolLeft;
			if ((actMemo.deltaVolLeft>0 && volL>actMemo.actVolumeLeft) ||
				(actMemo.deltaVolLeft<0 && volL<actMemo.actVolumeLeft))
			{
				volL = actMemo.actRampVolLeft = actMemo.actVolumeLeft;
				actMemo.deltaVolLeft = 0;
			}
			else
			{
				actMemo.actRampVolLeft += actMemo.deltaVolLeft;
			}
			// Volume Ramping (deClick) Right!
			long volR = actMemo.actRampVolRight;
			if ((actMemo.deltaVolRight>0 && volR>actMemo.actVolumeRight) ||
				(actMemo.deltaVolRight<0 && volR<actMemo.actVolumeRight))
			{
				volR = actMemo.actRampVolRight = actMemo.actVolumeRight;
				actMemo.deltaVolRight = 0;
			}
			else
			{
				actMemo.actRampVolRight += actMemo.deltaVolRight;
			}
			
			// Fit into volume for the two channels
			leftBuffer[i] += (int)((sample * volL) >> Helpers.MAXVOLUMESHIFT);
			rightBuffer[i]+= (int)((sample * volR) >> Helpers.MAXVOLUMESHIFT);

			// Now fit the loops of the current Sample
			fitIntoLoops(actMemo);
			if (actMemo.instrumentFinished) break;
		}
	}
909150
def mel(sr, n_dft, n_mels=128, fmin=0.0, fmax=None, htk=False, norm=1):
    
    return librosa.filters.mel(sr=sr, n_fft=n_dft, n_mels=n_mels,
                               fmin=fmin, fmax=fmax,
                               htk=htk, norm=norm).astype(K.floatx())