【问题标题】：Efficiently extract frequency of signal from FFT从 FFT 中有效提取信号频率
【发布时间】：2014-01-21 03:12:08
【问题描述】：

我正在使用 R 并尝试通过对每个声波（1000 个音频文件）应用快速傅里叶变换并识别频率来恢复频率（实际上，只是一个接近实际频率的数字）每个文件的幅度最高。我希望能够尽快恢复这些峰值频率。 FFT 方法是我最近了解的一种方法，我认为它应该适用于这项任务，但我对不依赖 FFT 的答案持开放态度。我已经尝试了几种应用 FFT 并获得最大幅度频率的方法，并且自从我使用第一种方法以来，我已经看到了显着的性能提升，但如果可能的话，我想进一步加快执行时间。

这里是示例数据：

s.rate<-44100                        # sampling frequency
t <- 2                               # seconds, for my situation, I've got 1000s of 1 - 5 minute files to go through
ind <- seq(s.rate*t)/s.rate          # time indices for each step
                                     # let's add two sin waves together to make the sound wave
f1 <- 600                            # Hz: freq of sound wave 1
y <- 100*sin(2*pi*f1*ind)            # sine wave 1
f2 <- 1500                           # Hz: freq of sound wave 2
z <- 500*sin(2*pi*f2*ind+1)          # sine wave 2
s <- y+z                             # the sound wave: my data isn't this nice, but I think this is an OK example

我尝试的第一种方法是使用 seewave 包中的 fpeaks 和 spec 函数，它似乎有效。但是，它的速度非常慢。

library(seewave)
fpeaks(spec(s, f=s.rate), nmax=1, plot=F) * 1000  # *1000 in order to recover freq in Hz
[1] 1494
# pretty close, quite slow

在阅读了更多内容后，我尝试了下一种方法，其中

spec(s, f=s.rate, plot=F)[which(spec(s, f=s.rate, plot=F)[,2]==max(spec(s, f=s.rate, plot=F)[,2])),1] * 1000 # again need to *1000 to get Hz
   x 
1494 
# pretty close, definitely faster

多看几圈后，我发现这种方法效果不错。

which(Mod(fft(s)) == max(abs(Mod(fft(s))))) * s.rate / length(s)
[1] 1500  
# recovered the exact frequency, and quickly!

以下是一些性能数据：

library(microbenchmark)
microbenchmark(
  WHICH.MOD = which(Mod(fft(s))==max(abs(Mod(fft(s))))) * s.rate / length(s),
  SPEC.WHICH = spec(s,f=s.rate,plot=F)[which(spec(s,f=s.rate,plot=F)[,2] == max(spec(s,f=s.rate,plot=F)[,2])),1] * 1000,   # this is spec from the seewave package
  # to recover a number around 1500, you have to multiply by 1000
  FPEAKS.SPEC = fpeaks(spec(s,f=s.rate),nmax=1,plot=F)[,1] * 1000, # fpeaks is from the seewave package... again, need to multiply by 1000
  times=10)

Unit: milliseconds
        expr       min        lq    median        uq       max neval
   WHICH.MOD     10.78     10.81     11.07     11.43     12.33    10
  SPEC.WHICH     64.68     65.83     66.66     67.18     78.74    10
 FPEAKS.SPEC 100297.52 100648.50 101056.05 101737.56 102927.06    10

好的解决方案是最快恢复接近（± 10 Hz）真实频率的频率。

更多上下文

我有很多文件（几个 GB），每个文件都包含一个每秒被调制几次的音调，有时信号实际上完全消失了，所以只有沉默。我想确定未调制音调的频率。我知道它们都应该低于 6000 Hz，但我不知道比这更准确。如果（大如果）我理解正确，我在这里有一个好的方法，这只是让它更快的问题。仅供参考，我以前没有数字信号处理方面的经验，因此除了有关如何更好地以编程方式处理此问题的建议外，我还欣赏与数学/方法相关的任何提示和指示。

【问题讨论】：

使用 FFT 的问题在于它假设输入是周期性的。对于信号快照中存在的大多数频率，通常情况并非如此。
@MatthewLundberg 我的理解是我有一个音调，它是一个相对恒定的频率，比如 800 ± 50 Hz ，它有时会被调制，但它是信号中存在的主要频率。这将被认为是周期性的，正确的，并且应该可以通过这种方法识别？如果没有，为什么不；我误会了什么？
我所说的周期性是指 FFT 认为给定的信号是背靠背重放的，在两个方向上都是无穷大。这会为除少数频率之外的所有频率引入边缘效应。这些边缘效应可能会或可能不会影响您的结果。
@MatthewLundberg 感谢您的提示；我将不得不目视检查多个文件的频谱图，以查看该方法是否有效或有问题。对于处理边缘效应或替代方法有什么建议吗？
在 FT 之前在每一端应用平滑窗口可能会有所帮助。小波变换试图解决这个问题，但是您在频域中失去了准确性（并且在时域中获得了准确性；请注意，FT 不提供时域信息）。我很好奇您在使用窗口或非窗口 FT 的真实数据中发现了什么。除了“Windowed FFT”之外，另一个可能有助于查找相关信息的搜索词是“Spectral Leakage”。

标签： r performance audio signal-processing fft

【解决方案1】：

在对这项任务和所涉及的一些术语有了更好的理解后，我发现了一些其他方法，我将在这里介绍这些方法。这些额外的方法允许窗口函数和更多，实际上，我的问题中最快的方法没有。我也只是通过将一些函数的结果分配给一个对象并索引该对象而不是再次运行该函数来加快速度

#i.e.
(ms<-meanspec(s,f=s.rate,wl=1024,plot=F))[which.max(ms[,2]),1]*1000 
# instead of 
meanspec(s,f=s.rate,wl=1024,plot=F)[which.max(meanspec(s,f=s.rate,wl=1024,plot=F)[,2]),1]*1000

我有我最喜欢的方法，但我欢迎建设性的警告、反馈和意见。

microbenchmark(
  WHICH.MOD = which((mfft<-Mod(fft(s)))[1:(length(s)/2)] == max(abs(mfft[1:(length(s)/2)]))) * s.rate / length(s),
  MEANSPEC = (ms<-meanspec(s,f=s.rate,wl=1024,plot=F))[which.max(ms[,2]),1]*1000,
  DFREQ.HIST = (h<-hist(dfreq(s,f=s.rate,wl=1024,plot=F)[,2],200,plot=F))$mids[which.max(h$density)]*1000,
  DFREQ.DENS = (dens <- density(dfreq(s,f=s.rate,wl=1024,plot=F)[,2],na.rm=T))$x[which.max(dens$y)]*1000,
  FPEAKS.MSPEC = fpeaks(meanspec(s,f=s.rate,wl=1024,plot=F),nmax=1,plot=F)[,1]*1000 , 
  times=100)

Unit: milliseconds
         expr       min        lq    median        uq      max neval
    WHICH.MOD  8.119499  8.394254  8.513992  8.631377 10.81916   100
     MEANSPEC  7.748739  7.985650  8.069466  8.211654 10.03744   100
   DFREQ.HIST  9.720990 10.186257 10.299152 10.492016 12.07640   100
   DFREQ.DENS 10.086190 10.413116 10.555305 10.721014 12.48137   100
 FPEAKS.MSPEC 33.848135 35.441716 36.302971 37.089605 76.45978   100

DFREQ.DENS 返回一个离实际值最远的频率值。其他方法返回接近实际值的值。

对于我的一个音频文件（即真实数据），性能结果有点不同（见下文）。上面使用的数据和下面用于性能数据的真实数据之间的一个潜在相关区别是，上面的数据只是一个数字向量，而我的真实数据存储在 Wave 对象中，这是来自 @ 的 S4 对象987654325@包。

library(Rmpfr) # to avoid an integer overflow problem in `WHICH.MOD`
microbenchmark(
  WHICH.MOD = which((mfft<-Mod(fft(d@left)))[1:(length(d@left)/2)] == max(abs(mfft[1:(length(d@left)/2)]))) * mpfr(s.rate,100) / length(d@left),
  MEANSPEC = (ms<-meanspec(d,f=s.rate,wl=1024,plot=F))[which.max(ms[,2]),1]*1000,
  DFREQ.HIST = (h<-hist(dfreq(d,f=s.rate,wl=1024,plot=F)[,2],200,plot=F))$mids[which.max(h$density)]*1000,
  DFREQ.DENS = (dens <- density(dfreq(d,f=s.rate,wl=1024,plot=F)[,2],na.rm=T))$x[which.max(dens$y)]*1000,
  FPEAKS.MSPEC = fpeaks(meanspec(d,f=s.rate,wl=1024,plot=F),nmax=1,plot=F)[,1]*1000 , 
  times=25) 

Unit: seconds
         expr      min       lq   median       uq      max neval
    WHICH.MOD 3.249395 3.320995 3.361160 3.421977 3.768885    25
     MEANSPEC 1.180119 1.234359 1.263213 1.286397 1.315912    25
   DFREQ.HIST 1.468117 1.519957 1.534353 1.563132 1.726012    25
   DFREQ.DENS 1.432193 1.489323 1.514968 1.553121 1.713296    25
 FPEAKS.MSPEC 1.207205 1.260006 1.277846 1.308961 1.390722    25

WHICH.MOD 实际上必须运行两次才能考虑左右音频通道（即我的数据是立体声），因此它需要的时间比输出指示的要长。注意：我需要使用Rmpfr 库，以便WHICH.MOD 方法处理我的真实数据，因为我遇到了整数溢出问题。

有趣的是，FPEAKS.MSPEC 在我的数据中表现得非常好，而且它似乎返回了一个非常准确的频率（基于我对频谱图的视觉检查）。 DFREQ.HIST 和 DFREQ.DENS 速度很快，但输出频率与我判断的实际值不太接近，而且都是相对丑陋的解决方案。到目前为止，我最喜欢的解决方案 MEANSPEC 使用 meanspec 和 which.max。我将其标记为答案，因为我没有任何其他答案，但请随时提供其他答案。我会投票给它，如果它提供更好的解决方案，我可能会选择它作为答案。

【讨论】：