Dapatkan nilai puncak sinyal jika frekuensinya terletak di antara dua pusat bin

Silakan anggap sebagai berikut:

• Frekuensi fundamental sinyal telah diperkirakan menggunakan FFT dan beberapa metode estimasi frekuensi dan terletak di antara dua pusat bin
• Frekuensi pengambilan sampel ditetapkan
• Upaya komputasi tidak menjadi masalah

Mengetahui frekuensi, apa cara paling akurat untuk memperkirakan nilai puncak sinyal fundamental yang sesuai?

Salah satu cara mungkin dengan mem-pad sinyal waktu untuk meningkatkan resolusi FFT sehingga pusat bin akan lebih dekat ke frekuensi yang diperkirakan. Dalam skenario ini, satu hal yang saya tidak yakin tentang adalah apakah saya dapat nol-pad sebanyak yang saya inginkan atau jika ada beberapa kelemahan dalam melakukannya. Yang lain adalah tempat bin center yang harus saya pilih setelah zero padding sebagai yang saya dapatkan nilai puncaknya (karena orang mungkin tidak mengenai frekuensi yang diinginkan, bahkan setelah zero padding).

Namun, saya juga bertanya-tanya apakah ada metode lain yang dapat memberikan hasil yang lebih baik, yaitu penduga yang menggunakan nilai puncak dari dua pusat bin di sekitarnya untuk memperkirakan nilai puncak pada frekuensi yang diinginkan.

Algoritma pertama yang muncul dalam pikiran adalah Algoritma Goertzel . Algoritma itu biasanya mengasumsikan bahwa frekuensi yang menarik adalah kelipatan bilangan bulat dari frekuensi dasar. Namun, makalah ini menerapkan algoritma (umum) untuk kasus yang Anda minati.

Masalah lain adalah bahwa model sinyal tidak benar. Itu menggunakan 2*%pi*(1:siglen)*(Fc/siglen). Itu harus digunakan 2*%pi*(0:siglen-1)*(Fc/siglen)untuk fase untuk keluar dengan benar.

Saya juga berpikir ada masalah dengan frekuensinya Fc=21.3sangat rendah. Sinyal bernilai real frekuensi rendah cenderung menunjukkan bias ketika datang ke masalah estimasi fase / frekuensi.

Saya juga mencoba pencarian grid kasar untuk estimasi fase, dan memberikan jawaban yang sama dengan algoritma Goertzel.

Di bawah ini adalah plot yang menunjukkan bias pada kedua taksiran (Goertzel: biru, Kasar: merah) untuk dua frekuensi berbeda: Fc=21.3(padat) dan Fc=210.3(putus-putus). Seperti yang Anda lihat bias untuk frekuensi yang lebih tinggi jauh lebih sedikit.

$x$$2\pi$

Jika Anda bersedia menggunakan beberapa nampan FFT yang bersebelahan, bukan hanya 2, maka interpolasi Sinc berjendela antara hasil nampan kompleks dapat menghasilkan perkiraan yang sangat akurat, tergantung pada lebar jendela.

Interpolasi Windowed Sinc umumnya ditemukan dalam upampler audio berkualitas tinggi, sehingga makalah tentang subjek tersebut akan memiliki formula interpolasi yang sesuai dengan analisis kesalahan.

[1] JL Flanagan dan RM Golden, "Phase vocoder," Bell Systems Technical Journal, vol. 45, hlm. 1493–1509, 1966.

[2] K. Dressler, "Ekstraksi sinusoid menggunakan implementasi FFT multi-resolusi yang efisien," di Proc. Int 9. Conf. tentang Efek Audio Digital (DAFx-06), Montreal, Kanada, September 2006, hlm. 247–252.

Salah satu metode adalah menemukan maksimum dan sesuai parabola tentang hal itu, dan kemudian menggunakan maksimum parabola sebagai perkiraan frekuensi dan besarnya. Anda dapat membaca semua tentang di sini: https://ccrma.stanford.edu/~jos/sasp/Sinusoidal_Peak_Interpolation.html

Saya memiliki banyak kesulitan dengan masalah yang tepat ini beberapa tahun yang lalu.

Saya memposting pertanyaan ini:

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

Saya akhirnya melakukan perhitungan dari awal, dan memposting jawaban untuk pertanyaan saya sendiri.

Saya terkejut bahwa saya tidak dapat menemukan penjelasan serupa di Internet.

Saya akan mengirim jawabannya lagi di sini; perhatikan bahwa kode ini dirancang untuk skenario di mana saya tumpang tindih jendela FFT saya dengan 4x.

π

Teka-teki ini membutuhkan dua kunci untuk membukanya.

• Kunci pertama adalah untuk memahami bagaimana tumpang tindih jendela FFT memperkenalkan rotasi pada fase bin.

• Kunci kedua berasal dari Grafik 3.3 & 3.4 di sini (terima kasih kepada Stephan Bernsee untuk izin untuk menyalin foto di sini).

Grafik 3.3:

Grafik 3.4:

Kode:

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

Kode python ini akan memberikan Anda hasil yang sangat akurat (saya menggunakannya untuk banyak not musik dan mendapatkan kesalahan semitone kurang dari 0,01%) dengan interpolasi parabola (metode yang berhasil digunakan oleh McAulay Quatieri, Serra, dll. Dalam harmonik + residual teknik pemisahan)

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

Frekuensi yang Anda hadapi (21,3Hz sampel pada 8kHz) sangat rendah. Karena ini adalah sinyal bernilai nyata, mereka akan menunjukkan bias dalam estimasi fase untuk ** frekuensi ** apa pun.

Gambar ini menunjukkan plot bias ( phase_est - phase_orig) untuk Fc = 210.3;(berwarna merah) versus bias untuk Fc = 21.3;. Seperti yang Anda lihat, offset jauh lebih penting untuk 21.3kasing.

Pilihan lain adalah mengurangi laju sampling Anda. Kurva hijau menunjukkan bias untuk Fs = 800bukan 8000.