Scipy：几个选定频率的傅立叶变换

Posted 2023-03-12

【中文标题】Scipy：几个选定频率的傅立叶变换

【英文标题】：Scipy : fourier transform of a few selected frequencies

【发布时间】：2012-11-10 02:14:13

【问题描述】：

我在一个信号上使用scipy.fft，带有一个移动窗口来绘制随时间变化的频率幅度（这里是an example，时间在X，频率在Y，幅度是颜色）。

但是，我只对少数几个频率感兴趣（大约 3、4 个频率）。使用 FFT 似乎我不能只选择我想要的频率（因为显然频率范围是由算法决定的），所以我计算了很多无用的东西，如果信号太长了。

我该怎么办？我是否必须使用自定义傅立叶变换 - 在这种情况下，欢迎提供良好实现的链接 - 还是有 scipy 方式？

---

编辑

@jfaller 回答后，我决定（尝试）实现 Goertzel 算法。我想出了这个：https://gist.github.com/4128537，但它不起作用（频率 440 没有出现，没关系峰值，我没有费心应用适当的窗口）。有什么帮助！？我不擅长 DSP。

【问题讨论】：

其实，短时傅里叶变换，看这个问题***.com/questions/2459295/stft-and-istft-in-python可能是我需要的...

这不是你说的吗？即滑动窗口的 DFT。

哦...没错，我正在使用scipy.fft，而那家伙也在使用相同的。所以有同样的问题。

信号需要在窗口上方大致静止。这是窗口大小的上限。下限是所需的频率分辨率。希望它们是一致的。

一般来说，简单地计算一个频率上的 DFT 是一个 N 阶运算。如果要计算 N 个点上 M 个频率的 DFT 系数，并且 M k n)) 用于您感兴趣的 k。

【参考方案1】：

您确实希望使用 Goertzel 算法：http://en.wikipedia.org/wiki/Goertzel_algorithm。基本上，它是单点的 FFT，如果您只需要信号中有限数量的频率，它就很有效。如果您在从 Wikipedia 中分离算法时遇到问题，请回复，我会帮助您。此外，如果您搜索一些资源，则存在用 python 编写的 DTMF 解码器（按键式电话解码器）。你可以看看他们是怎么做的。

【讨论】：

Arrr ...我不能说我没有尝试过 (gist.github.com/4128537)。但肯定有问题。我不太擅长 DSP :( 任何提示都会很棒！

好的，我晚上要登出，但我今晚或明天会看看，并更新我的评论。同时，从***链接的C代码还不错：netwerkt.wordpress.com/2011/08/25/goertzel-filter

@sebpiq：您是否考虑过使用scipy.signal.lfilter 和lfiltic？这将是应用级联 IIR、FIR 滤波器的快速方法。

jfaller：这篇文章看起来不错。我正在看一看。 @eryksun ：感谢您的提示，我将对其进行优化。

【参考方案2】：

在@jfaller、@eryksun 的大力帮助下……我在 Python 中实现了一个简单的 Goertzel 算法：https://gist.github.com/4128537。我会在这里重新粘贴。正如@eryksun 提到的，可以使用scipy.signal.lfilter 和lfiltic 对其进行优化，但目前我会坚持使用它，因为我可能想用另一种语言实现它：

import math

def goertzel(samples, sample_rate, *freqs):
    """
    Implementation of the Goertzel algorithm, useful for calculating individual
    terms of a discrete Fourier transform.

    `samples` is a windowed one-dimensional signal originally sampled at `sample_rate`.

    The function returns 2 arrays, one containing the actual frequencies calculated,
    the second the coefficients `(real part, imag part, power)` for each of those frequencies.
    For simple spectral analysis, the power is usually enough.

    Example of usage : 

        # calculating frequencies in ranges [400, 500] and [1000, 1100]
        # of a windowed signal sampled at 44100 Hz

        freqs, results = goertzel(some_samples, 44100, (400, 500), (1000, 1100))
    """
    window_size = len(samples)
    f_step = sample_rate / float(window_size)
    f_step_normalized = 1.0 / window_size

    # Calculate all the DFT bins we have to compute to include frequencies
    # in `freqs`.
    bins = set()
    for f_range in freqs:
        f_start, f_end = f_range
        k_start = int(math.floor(f_start / f_step))
        k_end = int(math.ceil(f_end / f_step))

        if k_end > window_size - 1: raise ValueError('frequency out of range %s' % k_end)
        bins = bins.union(range(k_start, k_end))

    # For all the bins, calculate the DFT term
    n_range = range(0, window_size)
    freqs = []
    results = []
    for k in bins:

        # Bin frequency and coefficients for the computation
        f = k * f_step_normalized
        w_real = 2.0 * math.cos(2.0 * math.pi * f)
        w_imag = math.sin(2.0 * math.pi * f)

        # Doing the calculation on the whole sample
        d1, d2 = 0.0, 0.0
        for n in n_range:
            y  = samples[n] + w_real * d1 - d2
            d2, d1 = d1, y

        # Storing results `(real part, imag part, power)`
        results.append((
            0.5 * w_real * d1 - d2, w_imag * d1,
            d2**2 + d1**2 - w_real * d1 * d2)
        )
        freqs.append(f * sample_rate)
    return freqs, results


if __name__ == '__main__':
    # quick test
    import numpy as np
    import pylab

    # generating test signals
    SAMPLE_RATE = 44100
    WINDOW_SIZE = 1024
    t = np.linspace(0, 1, SAMPLE_RATE)[:WINDOW_SIZE]
    sine_wave = np.sin(2*np.pi*440*t) + np.sin(2*np.pi*1020*t)
    sine_wave = sine_wave * np.hamming(WINDOW_SIZE)
    sine_wave2 = np.sin(2*np.pi*880*t) + np.sin(2*np.pi*1500*t)
    sine_wave2 = sine_wave2 * np.hamming(WINDOW_SIZE)

    # applying Goertzel on those signals, and plotting results
    freqs, results = goertzel(sine_wave, SAMPLE_RATE, (400, 500),  (1000, 1100))

    pylab.subplot(2, 2, 1)
    pylab.title('(1) Sine wave 440Hz + 1020Hz')
    pylab.plot(t, sine_wave)

    pylab.subplot(2, 2, 3)
    pylab.title('(1) Goertzel Algo, freqency ranges : [400, 500] and [1000, 1100]')
    pylab.plot(freqs, np.array(results)[:,2], 'o')
    pylab.ylim([0,100000])

    freqs, results = goertzel(sine_wave2, SAMPLE_RATE, (400, 500),  (1000, 1100))

    pylab.subplot(2, 2, 2)
    pylab.title('(2) Sine wave 660Hz + 1200Hz')
    pylab.plot(t, sine_wave2)

    pylab.subplot(2, 2, 4)
    pylab.title('(2) Goertzel Algo, freqency ranges : [400, 500] and [1000, 1100]')
    pylab.plot(freqs, np.array(results)[:,2], 'o')
    pylab.ylim([0,100000])

    pylab.show()