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 

I는 seewave 패키지에서 fpeaksspec 및 기능을 사용하려고 한 첫 번째 방법은, 작동하는 것으로 보인다. 그러나, 그것은 매우 어렵습니다.

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 

을 시도, 나는 합리적으로 잘 작동하는이 방법을 발견했다.

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에서) 복구하는 사람이 될 것입니다 : 여기

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

약간의 성능 데이터입니다. 각각 여러 번에게 두 번째 변조, 그냥 침묵이되도록 가끔 신호가 실제로 완전히 사라됩니다 톤을 포함

나는 많은 파일 (몇 GB를)를 가지고

더 많은 상황. 나는 변조되지 않은 음색의 주파수를 확인하고 싶습니다. 나는 그들이 모두 6000 Hz 미만의 어딘가에 있어야한다는 것을 알고 있지만, 그것보다 정확하게는 모른다. (큰 경우) 제가 올바르게 이해한다면, 여기에 OK 접근법이 있습니다. 단지 더 빨리 만드는 문제 일뿐입니다. 필자는 디지털 신호 처리에 대한 이전 경험이 없으므로 프로그래밍 방식으로 접근하는 것이 더 나은 방법에 대한 조언 외에도 수학/방법과 관련된 팁과 포인터에 감사드립니다.

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 

을 기능을 실행하여 물건을 조금 빨라,하지만 건설적인 경고를 환영 피드백 , 그리고 의견.

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 실제 값과 주파수 값 먼을 반환합니다. 다른 접근법은 실제 값에 가까운 값을 반환합니다.

내 오디오 파일 중 하나 (즉, 실제 데이터)와 성능 결과가 약간 다릅니다 (아래 참조).위에서 사용 된 데이터와 아래 성능 데이터에 사용 된 실제 데이터 사이의 잠재적 인 관련성의 차이점 중 하나는 데이터가 숫자의 벡터 일 뿐이며 실제 데이터는 tuneR 패키지의 S4 객체 인 Wave 개체에 저장된다는 것입니다.

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

실제로는 더 이상 출력이 표시보다 소요되므로 (즉, 내 데이터가 스테레오) 왼쪽과 오른쪽 오디오 채널을 고려하여 두 번 실행해야합니다. 참고 : 정수 넘침 문제가있어서 WHICH.MOD 접근 방식을 실제 데이터와 함께 사용하려면 Rmpfr 라이브러리를 사용해야했습니다.

흥미롭게도 FPEAKS.MSPEC은 내 데이터와 함께 매우 잘 수행되었으며 꽤 정확한 주파수를 반환하는 것으로 보입니다 (스펙트로 그램의 시각적 검사를 기반으로 함). DFREQ.HIST 및 DFREQ.DENS은 빠르지 만 출력 빈도가 실제 가치라고 판단한 것과 거의 같지 않으며 둘 다 비교적 못생긴 해결책입니다. 지금까지 내가 가장 좋아하는 솔루션 인 MEANSPEC은 meanspec과 which.max을 사용합니다. 다른 답변이 없기 때문에 이것을 답으로 표시 하겠지만 다른 대답은 자유롭게 제공하십시오. 나는 더 나은 해결책을 제공한다면 그것을 투표하고 대답으로 선택할 것입니다.