Python signal.kaiser函数代码示例

def taper(flag, m, samples):
    # m = samples for taper
    # samples = total samples
    window = kaiser(2*m+1, beta=14)

    if flag == 'front':
        # cut and replace the second half of window with 1s
        ones = np.ones(samples-m-1)
        window = window[0:(m+1)]
        window = np.concatenate([window, ones])

    elif flag == 'end':
        # cut and replace the first half of window with 1s
        ones = np.ones(samples-m-1)
        window = window[(m+1):]
        window = np.concatenate([ones, window])

    elif flag == 'all':
        ones = np.ones(samples-2*m-1)
        window = np.concatenate([window[0:(m+1)], ones, window[(m+1):]])

    # avoid concatenate error
    if window.size < samples:
        window = np.append(window, 1)

    if window.size != samples:
        print(window.size)
        print(samples)
        print("[ERROR]: taper and data do not have the same number of samples.")
        window = np.ones(samples)

    return window

def apply_kaiserbessel_window(X, alpha=6.5):
    """
    Apply a Kaiser-Bessel window to X.

    Parameters
    ----------
    X : ndarray, shape=(n_samples, n_features)
        Input array of samples

    alpha : float, optional (default=6.5)
        Tuning parameter for Kaiser-Bessel function. alpha=6.5 should make
        perfect reconstruction possible for MDCT.

    Returns
    -------
    X_windowed : ndarray, shape=(n_samples, n_features)
        Windowed version of X.
    """
    beta = np.pi * alpha
    win = sg.kaiser(X.shape[1], beta)
    row_stride = 0
    col_stride = win.itemsize
    strided_win = as_strided(win, shape=X.shape,
                             strides=(row_stride, col_stride))
    return X * strided_win

def firwin_kaiser_lpf(f_pass, f_stop, d_stop, fs = 1.0, N_bump=0):
    """
    Design an FIR lowpass filter using the sinc() kernel and
    a Kaiser window. The filter order is determined based on 
    f_pass Hz, f_stop Hz, and the desired stopband attenuation
    d_stop in dB, all relative to a sampling rate of fs Hz.
    Note: the passband ripple cannot be set independent of the
    stopband attenuation.

    Mark Wickert October 2016
    """
    wc = 2*np.pi*(f_pass + f_stop)/2/fs
    delta_w = 2*np.pi*(f_stop - f_pass)/fs
    # Find the filter order
    M = np.ceil((d_stop - 8)/(2.285*delta_w))
    # Adjust filter order up or down as needed
    M += N_bump
    N_taps = M + 1
    # Obtain the Kaiser window
    beta = signal.kaiser_beta(d_stop)
    w_k = signal.kaiser(N_taps,beta)
    n = np.arange(N_taps)
    b_k = wc/np.pi*np.sinc(wc/np.pi*(n-M/2)) * w_k
    b_k /= np.sum(b_k)
    print('Kaiser Win filter taps = %d.' % N_taps)
    return b_k

 def test_spectrogram_energy_conservation(self):
     """Test the energy conservation property of the spectrogram."""
     signal, _ = fmlin(128, 0.1, 0.4)
     window = kaiser(17, 3 * np.pi)
     tfr, ts, freqs = cohen.Spectrogram(signal, n_fbins=64, fwindow=window).run()
     e_sig = (np.abs(signal) ** 2).sum()
     self.assertAlmostEqual(tfr.sum().sum() / 64, e_sig)

def THD(signal, sample_rate):
    """Measure the THD for a signal
    
    This function is not yet trustworthy.
    
    Returns the estimated fundamental frequency and the measured THD,
    calculated by finding peaks in the spectrum.
    
    There are two definitions for THD, a power ratio or an amplitude ratio
    When finished, this will list both
    
    """
    # Get rid of DC and window the signal
    signal -= mean(signal) # TODO: Do this in the frequency domain, and take any skirts with it?
    windowed = signal * kaiser(len(signal), 100)
    
    # Find the peak of the frequency spectrum (fundamental frequency)
    f = rfft(windowed)
    i = argmax(abs(f))
    true_i = parabolic(log(abs(f)), i)[0]
    print 'Frequency: %f Hz' % (sample_rate * (true_i / len(windowed)))
    
    print 'fundamental amplitude: %.3f' % abs(f[i])
    
    # Find the values for the first 15 harmonics.  Includes harmonic peaks only, by definition
    # TODO: Should peak-find near each one, not just assume that fundamental was perfectly estimated.
    # Instead of limited to 15, figure out how many fit based on f0 and sampling rate and report this "4 harmonics" and list the strength of each
    for x in range(2, 15):
        print '%.3f' % abs(f[i * x]),
   
    THD = sum([abs(f[i*x]) for x in range(2,15)]) / abs(f[i])
    print '\nTHD: %f%%' % (THD * 100)
    return

def freq_from_hps(signal, fs):
    """
    Estimate frequency using harmonic product spectrum

    Low frequency noise piles up and overwhelms the desired peaks

    Doesn't work well if signal doesn't have harmonics
    """
    signal = asarray(signal) + 0.0

    N = len(signal)
    signal -= mean(signal)  # Remove DC offset

    # Compute Fourier transform of windowed signal
    windowed = signal * kaiser(N, 100)

    # Get spectrum
    X = log(abs(rfft(windowed)))

    # Remove mean of spectrum (so sum is not increasingly offset
    # only in overlap region)
    X -= mean(X)

    # Downsample sum logs of spectra instead of multiplying
    hps = copy(X)
    for h in range(2, 9):  # TODO: choose a smarter upper limit
        dec = decimate(X, h, zero_phase=True)
        hps[:len(dec)] += dec

    # Find the peak and interpolate to get a more accurate peak
    i_peak = argmax(hps[:len(dec)])
    i_interp = parabolic(hps, i_peak)[0]

    # Convert to equivalent frequency
    return fs * i_interp / N  # Hz

def freq_from_fft(signal, fs):
    """
    Estimate frequency from peak of FFT

    Pros: Accurate, usually even more so than zero crossing counter
    (1000.000004 Hz for 1000 Hz, for instance).  Due to parabolic
    interpolation being a very good fit for windowed log FFT peaks?
    https://ccrma.stanford.edu/~jos/sasp/Quadratic_Interpolation_Spectral_Peaks.html
    Accuracy also increases with signal length

    Cons: Doesn't find the right value if harmonics are stronger than
    fundamental, which is common.
    """
    signal = asarray(signal)

    N = len(signal)

    # Compute Fourier transform of windowed signal
    windowed = signal * kaiser(N, 100)
    f = rfft(windowed)

    # Find the peak and interpolate to get a more accurate peak
    i_peak = argmax(abs(f))  # Just use this value for less-accurate result
    i_interp = parabolic(log(abs(f)), i_peak)[0]

    # Convert to equivalent frequency
    return fs * i_interp / N  # Hz

def firwin_kaiser_hpf(f_stop, f_pass, d_stop, fs = 1.0, N_bump=0):
    """
    Design an FIR highpass filter using the sinc() kernel and
    a Kaiser window. The filter order is determined based on 
    f_pass Hz, f_stop Hz, and the desired stopband attenuation
    d_stop in dB, all relative to a sampling rate of fs Hz.
    Note: the passband ripple cannot be set independent of the
    stopband attenuation.

    Mark Wickert October 2016
    """
    # Transform HPF critical frequencies to lowpass equivalent
    f_pass_eq = fs/2. - f_pass
    f_stop_eq = fs/2. - f_stop
    # Design LPF equivalent
    wc = 2*np.pi*(f_pass_eq + f_stop_eq)/2/fs
    delta_w = 2*np.pi*(f_stop_eq - f_pass_eq)/fs
    # Find the filter order
    M = np.ceil((d_stop - 8)/(2.285*delta_w))
    # Adjust filter order up or down as needed
    M += N_bump
    N_taps = M + 1
    # Obtain the Kaiser window
    beta = signal.kaiser_beta(d_stop)
    w_k = signal.kaiser(N_taps,beta)
    n = np.arange(N_taps)
    b_k = wc/np.pi*np.sinc(wc/np.pi*(n-M/2)) * w_k
    b_k /= np.sum(b_k)
    # Transform LPF equivalent to HPF
    n = np.arange(len(b_k))
    b_k *= (-1)**n
    print('Kaiser Win filter taps = %d.' % N_taps)
    return b_k

def freq_from_hps(signal, fs):
    """Estimate frequency using harmonic product spectrum
    
    Low frequency noise piles up and overwhelms the desired peaks
    """
    N = len(signal)
    signal -= mean(signal) # Remove DC offset
    
    # Compute Fourier transform of windowed signal
    windowed = signal * kaiser(N, 100)
    
    # Get spectrum
    X = log(abs(rfft(windowed)))
    
    # Downsample sum logs of spectra instead of multiplying
    hps = copy(X)
    for h in arange(2, 9): # TODO: choose a smarter upper limit
        dec = decimate(X, h)
        hps[:len(dec)] += dec
    
    # Find the peak and interpolate to get a more accurate peak
    i_peak = argmax(hps[:len(dec)])
    i_interp = parabolic(hps, i_peak)[0]
    
    # Convert to equivalent frequency
    return fs * i_interp / N # Hz

 def test_spectrogram_linearity(self):
     signal, _ = fmlin(128, 0.1, 0.4)
     window = kaiser(17, 3 * np.pi)
     tfr1, _, _ = cohen.spectrogram(signal, n_fbins=64, window=window)
     tfr2, _, _ = cohen.spectrogram(signal * 2, n_fbins=64, window=window)
     x = np.sum(np.sum(tfr2))
     y = np.sum(np.sum(tfr1))
     self.assertEqual(x / y, 4)

def calc_freq_fft(signal):
  N = len(signal)
  # Compute Fourier transform of windowed signal
  windowed = signal * ss.kaiser(N, 100)
  f = np.fft.rfft(windowed)
  # Find the peak and interpolate to get a more accurate peak
  i_peak = np.argmax(abs(f)) # Just use this value for less-accurate result
  i_interp = parabolic(np.log(abs(f)), i_peak)[0]
  # Convert to equivalent frequency
  return sample_rate * i_interp / N # Hz

 def test_spectrogram_linearity(self):
     """Test the linearity property of the spectrogram."""
     signal, _ = fmlin(128, 0.1, 0.4)
     window = kaiser(17, 3 * np.pi)
     tfr1, _, _ = cohen.Spectrogram(signal, n_fbins=64,
                                    fwindow=window).run()
     tfr2, _, _ = cohen.Spectrogram(signal * 2, n_fbins=64,
                                    fwindow=window).run()
     x = np.sum(np.sum(tfr2))
     y = np.sum(np.sum(tfr1))
     self.assertEqual(x / y, 4)

 def test_basic(self):
     assert_allclose(signal.kaiser(6, 0.5),
                     [0.9403061933191572, 0.9782962393705389,
                      0.9975765035372042, 0.9975765035372042,
                      0.9782962393705389, 0.9403061933191572])
     assert_allclose(signal.kaiser(7, 0.5),
                     [0.9403061933191572, 0.9732402256999829,
                      0.9932754654413773, 1.0, 0.9932754654413773,
                      0.9732402256999829, 0.9403061933191572])
     assert_allclose(signal.kaiser(6, 2.7),
                     [0.2603047507678832, 0.6648106293528054,
                      0.9582099802511439, 0.9582099802511439,
                      0.6648106293528054, 0.2603047507678832])
     assert_allclose(signal.kaiser(7, 2.7),
                     [0.2603047507678832, 0.5985765418119844,
                      0.8868495172060835, 1.0, 0.8868495172060835,
                      0.5985765418119844, 0.2603047507678832])
     assert_allclose(signal.kaiser(6, 2.7, False),
                     [0.2603047507678832, 0.5985765418119844,
                      0.8868495172060835, 1.0, 0.8868495172060835,
                      0.5985765418119844])

def CircularDeconvolution(x,y,dt,fmax,beta=None):
    
    from numpy.fft import rfft, irfft, fftshift, fft, ifft, fftfreq, ifftshift
    from numpy import floor,zeros,hstack,conj,pi,exp,linspace,ones,real,pi
    from scipy.signal import blackmanharris,kaiser,hann
    from matplotlib.pyplot import plot, show
    
    # mx = abs(x).argmax()
    # my = abs(y).argmax()
    # Nx = len(x)
   #  Ny = len(y)
    Nmin = len(x)+len(y)-1
    
   
    N = FFTLengthPower2((len(x)+len(y)-1))
    
    

    # X = rfft(hstack((x[mx::],zeros(abs(Ny-Nx)),x[0:mx-1])))
    X = fft(x,n=N)
    Sxx = fftshift(X*conj(X))
    # Y = rfft(y[m::],n=Ny)
    # Y = rfft(hstack((y[mx::],y[0:mx-1])))
    Y = fft(y,n=N)
    Syx = fftshift(Y*conj(X))
    
    f = fftshift(fftfreq(N,dt))
    
    fpass = [(f>=-fmax)&(f<=fmax)]
    
    H = Syx[fpass]/Sxx[fpass]
    
    
    
    if beta is None:
    
        H = hann(len(H))*H
        
    else:
        
        H = kaiser(len(H),beta)*H
        

    HH = zeros(N)+0.0*1j
    
    HH[fpass] = H.copy()
    
    H = ifftshift(HH)
  
    h = real(ifft(H))

    return h[0:Nmin], H

def kaiser_derived(M, beta):
    """ Return a Kaiser-Bessel derived window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero or less, an empty
        array is returned.
    beta : float
        Kaiser-Bessel window shape parameter.

    Returns
    -------
    w : ndarray
        The window, normalized to fulfil the Princen-Bradley condition.

    Notes
    -----
    This window is only defined for an even number of taps.

    References
    ----------
    .. [1] Wikipedia, "Kaiser window",
           https://en.wikipedia.org/wiki/Kaiser_window

    """
    M = int(M)
    try:
        from scipy.signal import kaiser_derived as scipy_kd
        return scipy_kd(M, beta)
    except ImportError:
        pass

    if M < 1:
        return np.array([])

    if M % 2:
        raise ValueError(
            "Kaiser Bessel Derived windows are only defined for even number "
            "of taps"
        )

    w = np.zeros(M)
    kaiserw = kaiser(M // 2 + 1, beta)
    csum = np.cumsum(kaiserw)
    halfw = np.sqrt(csum[:-1] / csum[-1])
    w[:M//2] = halfw
    w[-M//2:] = halfw[::-1]

    return w