Scipy:对一些选定频率进行傅立叶变换

Question

我在信号上使用scipy.fft，通过移动窗口绘制随时间变化的频率振幅(这里是一个示例，时间在X上，频率在Y上，振幅是颜色)。

然而，只有少数几个频率让我感兴趣(大约只有3，4个频率)。使用FFT，我似乎无法只选择我想要的频率(因为频率范围显然是由算法决定的)，所以我计算了很多无用的东西，如果信号太长，我的程序甚至会用MemoryError崩溃。

我该怎么办？我必须使用自定义傅立叶变换吗？在这种情况下，欢迎良好实现的链接？或者有scipy方法吗？

---

编辑

在@jfaller回答后，我决定(尝试)实现Goertzel算法。我想到了这个：https://gist.github.com/4128537但它不起作用(频率440没有出现，更不用说峰值了，我没有费心应用合适的窗口)。有人帮忙吗！？我不擅长DSP。

Answer 1

您确实希望使用Goertzel算法：http://en.wikipedia.org/wiki/Goertzel_algorithm.基本上，它是一个单点FFT，如果你只需要信号中有限数量的频率，它是有效的。如果你在从维基百科中提取算法时遇到困难，请ping回来，我会帮你的。此外，如果你在谷歌上搜索一些资源，就会发现存在用python编写的DTMF解码器(按键电话解码器)。你可以看看他们是怎么做的。

Answer 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()