Python utils.get_freqs函数代码示例

    def spectrum_fourier(self):
        """

        This is the spectrum estimated as the FFT of the time-series

        Returns
        -------
        (f,spectrum): f is an array with the frequencies and spectrum is the
        complex-valued FFT.

        """

        data = self.input.data
        sampling_rate = self.input.sampling_rate

        fft = fftpack.fft
        if np.any(np.iscomplex(data)):
            # Get negative frequencies, as well as positive:
            f = np.linspace(-sampling_rate/2., sampling_rate/2., data.shape[-1])
            spectrum_fourier = np.fft.fftshift(fft(data))
        else:
            f = tsu.get_freqs(sampling_rate, data.shape[-1])
            spectrum_fourier = fft(data)[..., :f.shape[0]]

        return f, spectrum_fourier

def test_cached_coherence():
    """Testing the cached coherence functions """
    NFFT = 64 #This is the default behavior
    n_freqs = NFFT//2 + 1
    ij = [(0,1),(1,0)]
    ts = np.loadtxt(os.path.join(test_dir_path,'tseries12.txt'))
    freqs,cache = tsa.cache_fft(ts,ij)

    #Are the frequencies the right ones?
    yield npt.assert_equal,freqs,ut.get_freqs(2*np.pi,NFFT)
                     
    #Check that the fft of the first window is what we expect:
    hann = mlab.window_hanning(np.ones(NFFT))
    w_ts = ts[0][:NFFT]*hann
    w_ft = np.fft.fft(w_ts)[0:n_freqs]

    #This is the result of the function:
    first_window_fft = cache['FFT_slices'][0][0]
    
    yield npt.assert_equal,w_ft,first_window_fft
    
    coh_cached = tsa.cache_to_coherency(cache,ij)[0,1]
    f,c = tsa.coherency(ts)
    coh_direct = c[0,1]

    yield npt.assert_almost_equal,coh_direct,coh_cached

    def filtered_fourier(self):

        """Filter the time-series by passing it to the Fourier domain and null
        out the frequency bands outside of the range [lb,ub] """

        freqs = tsu.get_freqs(self.sampling_rate,self.data.shape[-1])

        if self.ub is None:
            self.ub = freqs[-1]
            
        power = np.fft.fft(self.data)
        idx_0 = np.hstack([np.where(freqs<self.lb)[0],
                           np.where(freqs>self.ub)[0]])

        #Make sure that you keep the DC component:
        keep_dc = np.copy(power[...,0])
        power[...,idx_0] = 0
        power[...,-1*idx_0] = 0 #Take care of the negative frequencies
        power[...,0] = keep_dc #And put the DC back in when you're done:
        
        data_out = np.fft.ifft(power)

        data_out = np.real(data_out) #In order to make sure that you are not
                                      #left with float-precision residual
                                      #complex parts

        return ts.TimeSeries(data=data_out,
                                 sampling_rate=self.sampling_rate,
                                 time_unit=self.time_unit) 

def test_cached_coherence():
    """Testing the cached coherence functions """
    NFFT = 64  # This is the default behavior
    n_freqs = NFFT // 2 + 1
    ij = [(0, 1), (1, 0)]
    ts = np.loadtxt(os.path.join(test_dir_path, 'tseries12.txt'))
    freqs, cache = tsa.cache_fft(ts, ij)

    # Are the frequencies the right ones?
    npt.assert_equal(freqs, utils.get_freqs(2 * np.pi, NFFT))

    # Check that the fft of the first window is what we expect:
    hann = mlab.window_hanning(np.ones(NFFT))
    w_ts = ts[0][:NFFT] * hann
    w_ft = fftpack.fft(w_ts)[0:n_freqs]

    # This is the result of the function:
    first_window_fft = cache['FFT_slices'][0][0]

    npt.assert_equal(w_ft, first_window_fft)

    coh_cached = tsa.cache_to_coherency(cache, ij)[0, 1]
    f, c = tsa.coherency(ts)
    coh_direct = c[0, 1]

    npt.assert_almost_equal(coh_direct, coh_cached)

    # Only welch PSD works and an error is thrown otherwise. This tests that
    # the error is thrown:
    with pytest.raises(ValueError) as e_info:
        tsa.cache_fft(ts, ij, method=methods[2])

    # Take the method in which the window is defined on input:
    freqs, cache1 = tsa.cache_fft(ts, ij, method=methods[3])
    # And compare it to the method in which it isn't:
    freqs, cache2 = tsa.cache_fft(ts, ij, method=methods[4])
    npt.assert_equal(cache1, cache2)

    # Do the same, while setting scale_by_freq to False:
    freqs, cache1 = tsa.cache_fft(ts, ij, method=methods[3],
                                  scale_by_freq=False)
    freqs, cache2 = tsa.cache_fft(ts, ij, method=methods[4],
                                  scale_by_freq=False)
    npt.assert_equal(cache1, cache2)

    # Test cache_to_psd:
    psd1 = tsa.cache_to_psd(cache, ij)[0]
    # Against the standard get_spectra:
    f, c = tsa.get_spectra(ts)
    psd2 = c[0][0]

    npt.assert_almost_equal(psd1, psd2)

    # Test that prefer_speed_over_memory doesn't change anything:
    freqs, cache1 = tsa.cache_fft(ts, ij)
    freqs, cache2 = tsa.cache_fft(ts, ij, prefer_speed_over_memory=True)
    psd1 = tsa.cache_to_psd(cache1, ij)[0]
    psd2 = tsa.cache_to_psd(cache2, ij)[0]
    npt.assert_almost_equal(psd1, psd2)

 def __init__(self,time_series,lb=0,ub=None,boxcar_iterations=2):
     self.data = time_series.data 
     self.sampling_rate = time_series.sampling_rate
     self.freqs = tsu.get_freqs(self.sampling_rate,self.data.shape[-1])
     self.ub=ub
     self.lb=lb
     self.time_unit=time_series.time_unit
     self._boxcar_iterations=boxcar_iterations

    def spectrum_fourier(self):
        """ Simply the non-normalized Fourier transform for a real signal"""

        data = self.input.data
        sampling_rate = self.input.sampling_rate
        
        fft = np.fft.fft
        f = tsu.get_freqs(sampling_rate,data.shape[-1])
        spectrum_fourier = fft(data)[...,:f.shape[0]]
        return f,spectrum_fourier 

def correlation_spectrum(x1, x2, Fs=2 * np.pi, norm=False):
    """
    Calculate the spectral decomposition of the correlation.

    Parameters
    ----------
    x1,x2: ndarray
       Two arrays to be correlated. Same dimensions

    Fs: float, optional
       Sampling rate in Hz. If provided, an array of
       frequencies will be returned.Defaults to 2

    norm: bool, optional
       When this is true, the spectrum is normalized to sum to 1

    Returns
    -------
    f: ndarray
       ndarray with the frequencies

    ccn: ndarray
       The spectral decomposition of the correlation

    Notes
    -----

    This method is described in full in [Cordes2000]_

    .. [Cordes2000] D Cordes, V M Haughton, K Arfanakis, G J Wendt, P A Turski,
    C H Moritz, M A Quigley, M E Meyerand (2000). Mapping functionally related
    regions of brain with functional connectivity MR imaging. AJNR American
    journal of neuroradiology 21:1636-44

    """

    x1 = x1 - np.mean(x1)
    x2 = x2 - np.mean(x2)
    x1_f = np.fft.fft(x1)
    x2_f = np.fft.fft(x2)
    D = np.sqrt(np.sum(x1 ** 2) * np.sum(x2 ** 2))
    n = x1.shape[0]

    ccn = ((np.real(x1_f) * np.real(x2_f) +
           np.imag(x1_f) * np.imag(x2_f)) /
           (D * n))

    if norm:
        ccn = ccn / np.sum(ccn) * 2  # Only half of the sum is sent back
                                     # because of the freq domain symmetry.
                                     # XXX Does normalization make this
                                     # strictly positive?

    f = utils.get_freqs(Fs, n)
    return f, ccn[0:(n / 2 + 1)]

    def frequencies(self):
        """Get the central frequencies for the frequency bands, given the
           method of estimating the spectrum """

        self.method['Fs'] = self.method.get('Fs',self.input.sampling_rate)
        NFFT = self.method.get('NFFT',64)
        Fs = self.method.get('Fs')
        freqs = tsu.get_freqs(Fs,NFFT)
        lb_idx,ub_idx = tsu.get_bounds(freqs,self.lb,self.ub)
        
        return freqs[lb_idx:ub_idx]

    def frequencies(self):
        """Get the central frequencies for the frequency bands, given the
           method of estimating the spectrum """

        # Get the sampling rate from the seed time-series:
        self.method["Fs"] = self.method.get("Fs", self.seed.sampling_rate)
        NFFT = self.method.get("NFFT", 64)
        Fs = self.method.get("Fs")
        freqs = tsu.get_freqs(Fs, NFFT)
        lb_idx, ub_idx = tsu.get_bounds(freqs, self.lb, self.ub)

        return freqs[lb_idx:ub_idx]

def test_coherence_welch():
    """Tests that the code runs and that the resulting matrix is symmetric """  

    t = np.linspace(0,16*np.pi,1024)
    x = np.sin(t) + np.sin(2*t) + np.sin(3*t) + np.random.rand(t.shape[-1])
    y = x + np.random.rand(t.shape[-1])

    method = {"this_method":'welch',
              "NFFT":256,
              "Fs":2*np.pi}
    
    f,c = tsa.coherence(np.vstack([x,y]),csd_method=method)
    np.testing.assert_array_almost_equal(c[0,1],c[1,0])

    f_theoretical = ut.get_freqs(method['Fs'],method['NFFT'])
    npt.assert_array_almost_equal(f,f_theoretical)

def test_coherence_partial():
    """ Test partial coherence"""

    t = np.linspace(0,16*np.pi,1024)
    x = np.sin(t) + np.sin(2*t) + np.sin(3*t) + np.random.rand(t.shape[-1])
    y = x + np.random.rand(t.shape[-1])
    z = x + np.random.rand(t.shape[-1])

    method = {"this_method":'welch',
              "NFFT":256,
              "Fs":2*np.pi}
    f,c = tsa.coherence_partial(np.vstack([x,y]),z,csd_method=method)

    f_theoretical = ut.get_freqs(method['Fs'],method['NFFT'])
    npt.assert_array_almost_equal(f,f_theoretical)
    npt.assert_array_almost_equal(c[0,1],c[1,0])

def test_coherency():
    """
    Tests that the coherency algorithm runs smoothly, using the different
    csd routines, that the resulting matrix is symmetric and for the welch
    method, that the frequency bands in the output make sense
    """

    for method in methods:
        f, c = tsa.coherency(tseries, csd_method=method)

        npt.assert_array_almost_equal(c[0, 1], c[1, 0].conjugate())
        npt.assert_array_almost_equal(c[0, 0], np.ones(f.shape))

        if method is not None and method['this_method'] != "multi_taper_csd":
            f_theoretical = utils.get_freqs(method['Fs'], method['NFFT'])
            npt.assert_array_almost_equal(f, f_theoretical)

def test_coherency_welch():
    """Tests that the coherency algorithm runs smoothly, using the welch csd
    routine, that the resulting matrix is symmetric and that the frequency bands
    in the output make sense"""
    
    t = np.linspace(0,16*np.pi,1024)
    x = np.sin(t) + np.sin(2*t) + np.sin(3*t) + np.random.rand(t.shape[-1])
    y = x + np.random.rand(t.shape[-1])

    method = {"this_method":'welch',
              "NFFT":256,
              "Fs":2*np.pi}

    f,c = tsa.coherency(np.vstack([x,y]),csd_method=method)

    npt.assert_array_almost_equal(c[0,1],c[1,0].conjugate())
    npt.assert_array_almost_equal(c[0,0],np.ones(f.shape))
    f_theoretical = ut.get_freqs(method['Fs'],method['NFFT'])
    npt.assert_array_almost_equal(f,f_theoretical)

    def spectrum_fourier(self):
        """

        This is the spectrum estimated as the FFT of the time-series

        Returns
        -------
        (f,spectrum): f is an array with the frequencies and spectrum is the
        complex-valued FFT.

        """

        data = self.input.data
        sampling_rate = self.input.sampling_rate

        fft = fftpack.fft
        f = tsu.get_freqs(sampling_rate, data.shape[-1])
        spectrum_fourier = fft(data)[..., :f.shape[0]]
        return f, spectrum_fourier

def cache_fft(time_series, ij, lb=0, ub=None,
                  method=None, prefer_speed_over_memory=False,
                  scale_by_freq=True):
    """compute and cache the windowed FFTs of the time_series, in such a way
    that computing the psd and csd of any combination of them can be done
    quickly.

    Parameters
    ----------

    time_series : float array
       An ndarray with time-series, where time is the last dimension

    ij: list of tuples
      Each tuple in this variable should contain a pair of
      indices of the form (i,j). The resulting cache will contain the fft of
      time-series in the rows indexed by the unique elements of the union of i
      and j

    lb,ub: float
       Define a frequency band of interest, for which the fft will be cached

    method: dict, optional
        See :func:`get_spectra` for details on how this is used. For this set
        of functions, 'this_method' has to be 'welch'


    Returns
    -------
    freqs, cache

        where: cache =
             {'FFT_slices':FFT_slices,'FFT_conj_slices':FFT_conj_slices,
             'norm_val':norm_val}

    Notes
    -----

    - For these functions, only the Welch windowed periodogram ('welch') is
      available.

    - Detrending the input is not an option here, in order to save
      time on an empty function call.

    """
    if method is None:
        method = {'this_method': 'welch'}  # The default

    this_method = method.get('this_method', 'welch')

    if this_method == 'welch':
        NFFT = method.get('NFFT', 64)
        Fs = method.get('Fs', 2 * np.pi)
        window = method.get('window', mlab.window_hanning)
        n_overlap = method.get('n_overlap', int(np.ceil(NFFT / 2.0)))
    else:
        e_s = "For cache_fft, spectral estimation method must be welch"
        raise ValueError(e_s)
    time_series = utils.zero_pad(time_series, NFFT)

    #The shape of the zero-padded version:
    n_channels, n_time_points = time_series.shape

    # get all the unique channels in time_series that we are interested in by
    # checking the ij tuples
    all_channels = set()
    for i, j in ij:
        all_channels.add(i)
        all_channels.add(j)

    # for real time_series, ignore the negative frequencies
    if np.iscomplexobj(time_series):
        n_freqs = NFFT
    else:
        n_freqs = NFFT // 2 + 1

    #Which frequencies
    freqs = utils.get_freqs(Fs, NFFT)

    #If there are bounds, limit the calculation to within that band,
    #potentially include the DC component:
    lb_idx, ub_idx = utils.get_bounds(freqs, lb, ub)

    n_freqs = ub_idx - lb_idx
    #Make the window:
    if mlab.cbook.iterable(window):
        assert(len(window) == NFFT)
        window_vals = window
    else:
        window_vals = window(np.ones(NFFT, time_series.dtype))

    #Each fft needs to be normalized by the square of the norm of the window
    #and, for consistency with newer versions of mlab.csd (which, in turn, are
    #consistent with Matlab), normalize also by the sampling rate:

    if scale_by_freq:
        #This is the normalization factor for one-sided estimation, taking into
        #account the sampling rate. This makes the PSD a density function, with
        #units of dB/Hz, so that integrating over frequencies gives you the RMS
        #(XXX this should be in the tests!).
#.........这里部分代码省略.........