spectrum.py 文件源码

def blackman_tukey(x,
                   M,
                   L,
                   y=None,
                   window='boxcar',
                   window_args=[],
                   d=1,
                   full=False):
    """
    Compute the Blackman-Tukey cross power spectral density (PSD)
    estimate between the time-domain signals *x* and *y* (must be the
    same length as *x*). If *y* is not given, compute the power
    spectral density estimate of *x*.  Use the spectral window with
    identifier *window* (see the options in
    :func:scipy.`signal.get_window`, e.g., a tuple can be used to pass
    arguments to the window function) and length *M* (i.e., the
    maximum auto-correlation lag to include in the estimate). Compute
    the estimate at *L* uniformly spaced frequency samples where *d*
    is the time domain sample interval. If not *full*, return the
    tuple containing the length *L* PSD estimate and length *L*
    corresponding frequencies. If *full*, also return the estimated
    cross correlation and window function (i.e., a tuple with four
    elements).
    """
    N = len(x)
    assert M <= N
    if y is None:
        y = x
    else:
        assert len(y) == N
    Rxy = scipy.signal.correlate(x, y) / N
    Rxy_window = Rxy[(N - 1) - M:(N - 1) + M + 1]
    window = scipy.signal.get_window(window, 2*M + 1, fftbins=False)
    k_range = NP.arange(0, L)
    shift = NP.exp(2j * NP.pi * k_range * M / L)
    Sxy = NP.fft.fft(window * Rxy_window, n=L) * shift * d
    f = NP.fft.fftfreq(L, d=d)
    if full:
        return (Sxy, f, Rxy, window)
    else:
        return (Sxy, f)