r - 从 FFT 中有效提取信号频率

Question

我正在使用 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

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

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，但我不知道比这更准确。如果(大如果)我理解正确，我在这里有一个好的方法，这只是让它更快的问题。仅供引用，我以前没有数字信号处理方面的经验，因此除了有关如何更好地以编程方式处理此问题的建议外，我还欣赏与数学/方法相关的任何提示和指示。

Answer 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

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

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 中。对象，来自 tuneR 的 S4 对象包裹。

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实际上必须运行两次才能考虑左右音频 channel (即我的数据是立体声)，因此它需要的时间比输出指示的要长。注意:我需要使用 Rmpfr图书馆为了 WHICH.MOD处理我的真实数据的方法，因为我遇到整数溢出问题。

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