我正在使用R并试图通过将快速傅立叶变换应用于大量声波（1000s音频文件）来恢复频率（实际上，只是接近实际频率的一个数字）每个文件并确定每个文件最高幅度的频率。我希望能够尽快恢复这些峰值频率。 FFT方法是我最近了解到的一种方法，我认为它应该适用于此任务，但我愿意接受不依赖于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赫兹）的实际频率最快的人。

更多的上下文

我有许多文件（几个GBS），各包含被调制几次第二，有时信号实际上完全消失，从而有只是沉默的基调。我想识别未调制音调的频率。我知道他们都应该在6000赫兹以下的地方，但我不知道比这更精确。如果（大，如果）我理解正确，我在这里有一个好的方法，这只是一个让它更快的问题。就我而言，我以前没有数字信号处理的经验，所以我非常感谢与数学/方法相关的任何提示和指示，以及关于如何更好地以编程方式处理这些问题的建议。

来源

2014-01-21 Jota

使用FFT的问题是它假定输入是周期性的。信号快照中存在的大多数频率通常不是这种情况。 –

@MatthewLundberg我的理解是我有一个音调，它是一个相对固定的频率，比如800±50Hz，有时候会进行调制，但是是信号中存在的主要频率。这将被视为是周期性的，正确的，应该通过这种方法来识别？如果没有，为什么不;我误解了什么？ – Jota

我的意思是周期性的是，FFT认为给定的信号是在两个方向上背对背重复播放到无穷大。除了几个频率之外，这会引入边缘效应。这些边缘效应可能会或可能不会为您的结果着色。 –

在更好地理解了这个任务和一些涉及的术语后，我遇到了一些其他的方法，我将在这里介绍。这些额外的方法允许窗口函数和更多，真的，而我的问题中最快的方法不是。我也只是通过分配的一些函数的结果为对象和索引对象，而不是再次运行该功能

#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对象中，tuneR包中有一个S4对象。

library(Rmpfr) # to avoid an integer overflow problem in `WHICH.MOD` 
microbenchmark(
    WHICH.MOD = which((mfft<-Mod(fft([email protected])))[1:(length([email protected])/2)] == max(abs(mfft[1:(length([email protected])/2)]))) * mpfr(s.rate,100)/length([email protected]), 
    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。我将这个标记为答案，因为我没有任何其他答案，但随时可以提供其他答案。如果它能提供更好的解决方案，我会投票赞成，也可以选择它作为答案。

来源

2014-01-24 00:08:28 Jota

从FFT中有效提取信号的频率

更多的上下文

回答

相关问题