Uzyskaj wartość szczytową sygnału, jeśli jego częstotliwość leży między dwoma centrami bin

12

Przypuśćmy, że:

  • Częstotliwość podstawy sygnału została oszacowana za pomocą FFT i niektórych metod szacowania częstotliwości i leży pomiędzy dwoma centrami bin
  • Częstotliwość próbkowania jest stała
  • Wysiłek obliczeniowy nie stanowi problemu

Znając częstotliwość, jaki jest najdokładniejszy sposób oszacowania odpowiedniej wartości szczytowej fundamentalnych sygnałów?

Jednym ze sposobów może być zerowanie pola czasowego w celu zwiększenia rozdzielczości FFT, tak aby środek przedziału był bliższy oszacowanej częstotliwości. W tym scenariuszu jedną kwestią, o której nie jestem pewien, jest to, czy mogę zerować tak dużo, jak chcę, lub czy są w tym jakieś wady. Innym jest to, który środek przedziału powinienem wybrać po wypełnieniu zerowym, jako ten, od którego otrzymuję wartość szczytową (ponieważ nie można dokładnie trafić na częstotliwość będącą przedmiotem zainteresowania, nawet po wypełnieniu zerowym).

Zastanawiam się jednak również, czy istnieje inna metoda, która może zapewnić lepsze wyniki, na przykład estymator, który wykorzystuje wartości szczytowe otaczających dwóch centrów bin do oszacowania wartości szczytowej przy częstotliwości zainteresowania.

lR8n6i
źródło
2
zerowanie wypełniania przed FFT jest jednym ze sposobów. Kolejnym jest zastosowanie funkcji okna, która jest odpowiednia dla twoich neadów. Płaskie okno zostało zaprojektowane właśnie do tego celu. Oczywiście, jeśli znasz już dokładnie częstotliwość i interesuje Cię tylko jeden amputyd, prawdopodobnie są tańsze sposoby niż FFT.
sellibitze
1
wypełnienie zerowe nie jest wymagane: prosta interpolacja paraboliczna (z 3 punktami: imax-1, imax, imax + 1, gdzie imaxjest pik FFT) da dokładne wyniki
Basj,
Upewnij się, że funkcja interpolacji odpowiada funkcji okna. Flat-top jest trywialny, w przeciwnym razie potrzebujesz pasującej pary (np. Okno prostokątne + interpolacja sinc, okno gaussowskie + interpolacja gaussowska itp.)
finnw
@CedronDawg to pytanie i jego odpowiedzi są powiązane (ale nie takie same) z dokładną formułą częstotliwości. Może Cię to zainteresuje.
Fat32

Odpowiedzi:

5

Pierwszym algorytmem, który przychodzi mi do głowy, jest algorytm Goertzela . Algorytm ten zwykle zakłada, że ​​częstotliwość będąca przedmiotem zainteresowania jest całkowitą wielokrotnością częstotliwości podstawowej. Jednak ten artykuł stosuje (uogólniony) algorytm w przypadku, który Cię interesuje.


Innym problemem jest to, że model sygnału jest nieprawidłowy. Używa 2*%pi*(1:siglen)*(Fc/siglen). Należy użyć, 2*%pi*(0:siglen-1)*(Fc/siglen)aby faza wyszła poprawnie.

Myślę również, że istnieje problem z Fc=21.3bardzo niską częstotliwością . Sygnały o niskiej wartości rzeczywistej niskiej częstotliwości wykazują tendencyjność, jeśli chodzi o problemy z oszacowaniem fazy / częstotliwości.

Próbowałem również zgrubnego wyszukiwania siatki dla oszacowania fazy i daje to tę samą odpowiedź, co algorytm Goertzela.

Poniżej znajduje się wykres, który pokazuje odchylenie w obu oszacowaniach (Goertzel: niebieski, Gruboziarnisty: czerwony) dla dwóch różnych częstotliwości: Fc=21.3(stała) i Fc=210.3(przerywana). Jak widać odchylenie dla wyższej częstotliwości jest znacznie mniejsze.

Wykres osi jest początkową fazą zmieniającą się od 0 do .2 πx2)π

wprowadź opis zdjęcia tutaj

Peter K.
źródło
Właśnie przetestowałem kod algorytmu Goerzel na podstawie papieru. Wykorzystując wyjściową wartość DTFT, pik można uzyskać bardzo dokładnie. Istnieje jednak współczynnik skalowania wynoszący dokładnie 1000. Tak więc, jeśli pierwotny pik to 1,234, po Goerzelu będzie to 1234. Czy ktoś wie, skąd może pochodzić?
lR8n6i
W międzyczasie przeprowadziłem kilka badań. Prawdopodobnie ma to związek ze skalowaniem amplitudy: skalowanie amplituda w dziedzinie czasu = współczynnik w dziedzinie częstotliwości * 2 / N, gdzie N jest długością sygnału. Czy to założenie jest słuszne?
lR8n6i
Cześć! Właśnie odkryłem, że przy użyciu algorytmu Goertzela amplituda przy wynikowym złożonym współczynniku jest bardzo dokładna, ale faza jest całkowicie błędna. Czy ktoś ma pomysł, skąd to może pochodzić? Przez „fazę” rozumiem opóźnienie fazowe określone w podstawie pierwotnego sygnału.
lR8n6i
1
@ Rickson1982 Faza jest poprawna. Po prostu nie interpretujesz go poprawnie. :-) Pamiętaj: tj. jest o (90 stopni) od tego, czego oczekujesz. π/2grzech(ω0t+ϕ)jot2)[mi-jotϕδ~(ω+ω0+2)πk)-mi+jotϕδ~(ω-ω0+2)πk)]π/2)
Peter K.
4

Jeśli chcesz używać wielu sąsiednich pojemników FFT, a nie tylko 2, to interpolacja Sinc w okienku między złożonymi wynikami bin może dać bardzo dokładne oszacowanie, w zależności od szerokości okna.

Interpolacja okienkowana jest powszechnie spotykana w wysokiej jakości próbnikach audio, więc artykuły na ten temat będą miały odpowiednie formuły interpolacyjne z analizą błędów.

hotpaw2
źródło
Dziękuję za komentarz. Spróbuję również tego podejścia.
lR8n6i
4

Jeśli użyjesz Flanagana [1], jest on obliczany na podstawie różnicy faz kolejnych widm fazowych ϕϕ (Częstotliwość chwilowa), a jeśli rekonstruujesz wielkość za pomocą właściwego współczynnika (Natychmiastowa wielkość) [2], użyj znormalizowanej funkcji sinusa: A na koniec użyj interpolacji parabolicznej wokół wielkości szczytowej, możesz uzyskać niesamowite rezultaty, dziś uważam, że jest to najlepszy sposób, wykorzystałem go, a wyniki są zawsze bardzo solidny :-)

grzech(πx)(πx)

[1] JL Flanagan i RM Golden, „Phaser vocoder”, Bell Systems Technical Journal, vol. 45, s. 1493–1509, 1966.

[2] K. Dressler, „Ekstrakcja sinusoidalna z wykorzystaniem efektywnej implementacji FFT o wielu rozdzielczościach”, w Proc. 9th Int. Konf. w sprawie cyfrowych efektów dźwiękowych (DAFx-06), Montreal, Kanada, wrzesień 2006, s. 247–252.

Ederwander
źródło
Cześć! Wielkie dzięki za wszystkie komentarze. Rozszerzyłem swój kod (patrz poniżej), aby połączyć filtr Goertzela z paraboliczną interpolacją pików, aby uzyskać fazę. Jednak wyniki nadal nie są dokładne (+ - 3-4 stopni). Czy jest to tak bliskie, jak to możliwe, czy występują błędy w zrozumieniu lub kodowaniu?
lR8n6i
3

Kilka lat temu miałem wiele problemów z tym właśnie problemem.

Zadałem to pytanie:

/programming/4633203/extracting-precise-frequencies-from-fft-bins-using-phase-change-between-frames

Skończyłem robić obliczenia od zera i opublikowałem odpowiedź na moje własne pytanie.

Dziwi mnie, że nie udało mi się znaleźć podobnej ekspozycji w Internecie.

Ponownie opublikuję odpowiedź tutaj; zwróć uwagę, że kod jest przeznaczony na scenariusz, w którym nakładam okno FFT 4x.

π


Ta łamigłówka wymaga dwóch kluczy, aby ją odblokować.

Wykres 3.3:

wprowadź opis zdjęcia tutaj

Wykres 3.4:

wprowadź opis zdjęcia tutaj

Kod:

for (int k = 0; k <= fftFrameSize/2; k++) 
{
    // compute magnitude and phase 
    bins[k].mag = 2.*sqrt(fftBins[k].real*fftBins[k].real + fftBins[k].imag*fftBins[k].imag);
    bins[k].phase = atan2(fftBins[k].imag, fftBins[k].real);

    // Compute phase difference Δϕ fo bin[k]
    double deltaPhase;
    {
        double measuredPhaseDiff = bins[k].phase - gLastPhase[k];
        gLastPhase[k] = bins[k].phase;

        // Subtract expected phase difference <-- FIRST KEY
        // Think of a single wave in a 1024 float frame, with osamp = 4
        //   if the first sample catches it at phase = 0, the next will 
        //   catch it at pi/2 ie 1/4 * 2pi
        double binPhaseExpectedDiscrepancy = M_TWOPI * (double)k / (double)osamp;
        deltaPhase = measuredPhaseDiff - binPhaseExpectedDiscrepancy;

        // Wrap delta phase into [-Pi, Pi) interval 
        deltaPhase -= M_TWOPI * floor(deltaPhase / M_TWOPI + .5);
    }

    // say sampleRate = 40K samps/sec, fftFrameSize = 1024 samps in FFT giving bin[0] thru bin[512]
    // then bin[1] holds one whole wave in the frame, ie 44 waves in 1s ie 44Hz ie sampleRate / fftFrameSize
    double bin1Freq = (double)sampleRate / (double)fftFrameSize;
    bins[k].idealFreq = (double)k * bin1Freq;

    // Consider Δϕ for bin[k] between hops.
    // write as 2π / m.
    // so after m hops, Δϕ = 2π, ie 1 extra cycle has occurred   <-- SECOND KEY
    double m = M_TWOPI / deltaPhase;

    // so, m hops should have bin[k].idealFreq * t_mHops cycles.  plus this extra 1.
    // 
    // bin[k].idealFreq * t_mHops + 1 cycles in t_mHops seconds 
    //   => bins[k].actualFreq = bin[k].idealFreq + 1 / t_mHops
    double tFrame = fftFrameSize / sampleRate;
    double tHop = tFrame / osamp;
    double t_mHops = m * tHop;

    bins[k].freq = bins[k].idealFreq + 1. / t_mHops;
}
Liczba Pi
źródło
Interpolujesz częstotliwość, podczas gdy OP zna częstotliwość i chce interpolować amplitudę.
finnw
2

Ten kod python daje bardzo dokładny wynik (użyłem go do wielu nut i uzyskałem błędy poniżej 0,01% półtonu) z interpolacją paraboliczną (metoda z powodzeniem stosowana przez McAulay Quatieri, Serra itp. W harmonicznej + resztkowa techniki separacji)

import matplotlib.pyplot as plt
import numpy as np
from scipy.io.wavfile import read
from scipy.fftpack import fft, ifft
import math

(fs, x) = read('test.wav')
if (len(x.shape) == 2):    # if stereo we keep left channel only
 x = x[:,1]

n=x.size
freq = np.arange(n)*1.0/n*fs 
xfft = abs(fft(x))

imax=np.argmax(xfft)  
p=1.0/2*(xfft[imax-1]/xfft[imax]-xfft[imax+1]/xfft[imax])/(xfft[imax-1]/xfft[imax]-2+xfft[imax+1]/xfft[imax])   # parabolic interpolation 
print 'Frequence detectee avec interpolation parabolique :',(imax+p)*1.0/n*fs, 'Hz'
Basj
źródło
1
clear all
clc

for phase_orig = 0:pi/18:pi,

%% Specify and generate signal
Amp = 1;                     % Amplitude of signal
Fs = 8000;                   % samples per second
dt = 1/Fs;                   % seconds per sample
Fc = 21.3;                   % Hz
StopTime = 0.25;             % seconds
t = (0:dt:StopTime-dt)';     % seconds

siglen = length(t);
sig = Amp * 1.5 * sin(2*pi*(0:siglen-1)*(Fc/siglen) + phase_orig) + 1.5 * Amp * sin(2*pi*(0:siglen-1)*(Fc/siglen) * 3) ...
  + 1.5 * Amp * sin(2*pi*(0:siglen-1)*(Fc/siglen) * 5)+ 0.3 * Amp * sin(2*pi*(0:siglen-1)*(Fc/siglen) * 7) ...
  + 1.3 * Amp * sin(2*pi*(0:siglen-1)*(Fc/siglen) * 9)+ 1.4 * Amp * sin(2*pi*(0:siglen-1)*(Fc/siglen) * 11);

%% Estimate the peak value of the signals fundamental using Goertzel algorithm
peak = 0;
indvec = [Fc-1 Fc Fc+1];

% Check the input data
if ~isvector(sig) || isempty(sig)
  error('X must be a nonempty vector')
end

if ~isvector(indvec) || isempty(indvec)
  error('INDVEC must be a nonempty vector')
end
if ~isreal(indvec)
  error('INDVEC must contain real numbers')
end

% forcing x to be column
sig = reshape(sig,siglen,1);

% initialization
no_freq = length(indvec); %number of frequencies to compute
y = zeros(no_freq,1); %memory allocation for the output coefficients

% Computation via second-order system
% loop over the particular frequencies
for cnt_freq = 1:no_freq
  %for a single frequency:
  %a/ precompute the constants
  pik_term = 2*pi*(indvec(cnt_freq))/(siglen);
  cos_pik_term2 = cos(pik_term) * 2;
  cc = exp(-1i*pik_term); % complex constant
  %b/ state variables
  s0 = 0;
  s1 = 0;
  s2 = 0;
  %c/ 'main' loop
  for ind = 1:siglen-1 %number of iterations is (by one) less than the length of signal
    %new state
    s0 = sig(ind) + cos_pik_term2 * s1 - s2;  % (*)
    %shifting the state variables
    s2 = s1;
    s1 = s0;
  end
  %d/ final computations
  s0 = sig(siglen) + cos_pik_term2 * s1 - s2; %correspond to one extra performing of (*)
  y(cnt_freq) = s0 - s1*cc; %resultant complex coefficient

  %complex multiplication substituting the last iterationA
  %and correcting the phase for (potentially) non-integer valued
  %frequencies at the same time
  y(cnt_freq) = y(cnt_freq) * exp(-1i*pik_term*(siglen-1));
end

  % perfom amplitude scaling
  peak = abs(y(2)) * 2 / siglen

% perform parabolic interpolation to get the phase estimate
phase_orig=phase_orig*180/pi
ym1 = angle(unwrap(y(1)));
y0 = angle(unwrap(y(2)));
yp1 = angle(unwrap(y(3)));

p = (yp1 - ym1)/(2*(2*y0 - yp1 - ym1)); 
phase = y0 - 0.25*(ym1-yp1)*p;
phase_est = phase * 180/pi + 90;
phase_est = mod(phase_est+180,360)-180
end

Częstotliwości, z którymi mamy do czynienia (21,3 Hz próbkowanych przy 8 kHz) są bardzo niskie. Ponieważ są to sygnały o wartościach rzeczywistych, będą wykazywać błąd w szacowaniu faz dla ** dowolnej ** częstotliwości.

To zdjęcie pokazuje wykres odchylenia ( phase_est - phase_orig) dla Fc = 210.3;(na czerwono) w stosunku do odchylenia dla Fc = 21.3;. Jak widać, przesunięcie jest znacznie większe w 21.3przypadku.

Inną opcją jest zmniejszenie częstotliwości próbkowania. Zielona krzywa pokazuje odchylenie Fs = 800zamiast 8000.

wprowadź opis zdjęcia tutaj

lR8n6i
źródło
1
Dziękuję za aktualizację! Zobacz moją fabułę; Nadal uważam, że każdy estymator fazy będzie miał tendencję do tak niskiej częstotliwości. Jednym ze sposobów obejścia tego jest użycie znanej częstotliwości (jeśli jest znana!), Aby skorygować błąd oszacowania fazy za pomocą tabeli przeglądowej. Ale musisz być ostrożny: obciążenie będzie zmieniać się wraz z częstotliwością. Innym sposobem na to będzie zmniejszenie częstotliwości próbkowania.
Peter K.
1
Dziękuję też! Jeśli jednak używasz Fs = 8000 Hz i Fc = 210 zamiast 210,3, odchylenie wygląda jeszcze gorzej. Masz pomysł, skąd to się bierze?
lR8n6i
1
Erk! Brak pomysłu. FWIW, estymator Geortzel nie ma problemów: goertzel = atan(imag(y(2)),real(y(2)))*180/%pi + 90;. :-) Kopie trochę więcej. Obserwuj tą przestrzeń.
Peter K.
1
Interpolacja paraboliczna nie robi tego, co myślisz. W szczególności, jeśli zastąpi Obliczanie pprzy p2 = (abs(y(3)) - abs(y(1)))/(2*(2*abs(y(2)) - abs(y(3)) - abs(y(1)))); phase2 = y0 - 0.25*(ym1-yp1)*p2;czym można uzyskać znacznie lepsze odpowiedzi --- nawet Fc=210. Nie jestem wcale pewien, czy obecna wersja pda ci coś sensownego. Formuła interpolacji służy do interpolacji AMPLITUDY paraboli, ale pinterpoluje fazę, która jest po prostu ... dziwna.
Peter K.
1
Wszystko to jest OK, Z WYJĄTKIEM, że położenie piku ( p = (yp1 - ym1)/(2*(2*y0 - yp1 - ym1))) będzie przez pewien czas niepoprawne, jeśli użyjesz FAZY zamiast amplitud. Jest tak, ponieważ fazy mogą przeskakiwać wokół granicy +/- 180 stopni. Wszystko, co jest potrzebne, aby to naprawić dla fazy, to zmienić tę linię na moje p2obliczenia powyżej.
Peter K.