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

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.

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 $x$$2\pi$

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.

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 :-)

[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.

Jedną z metod jest znalezienie maksimum i dopasowanie do niego paraboli, a następnie użycie maksimum paraboli jako oszacowania częstotliwości i wielkości. Możesz przeczytać wszystko o tutaj: https://ccrma.stanford.edu/~jos/sasp/Sinusoidal_Peak_Interpolation.html

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ć.

• Pierwszym kluczem jest zrozumienie, w jaki sposób nakładanie się okna FFT wprowadza obrót w fazie bin.

• Drugi klucz pochodzi z wykresu 3.3 i 3.4 tutaj (dzięki Stephan Bernsee za zgodę na kopiowanie zdjęć tutaj).

Wykres 3.3:

Wykres 3.4:

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;
}

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'

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.