Получить пиковое значение сигнала, если его частота лежит между двумя центрами бина

12

Пожалуйста, предположим следующее:

  • Частота основного сигнала была оценена с использованием БПФ и некоторых методов оценки частоты и находится между двумя центрами бинов
  • Частота дискретизации фиксирована
  • Вычислительные усилия не проблема

Зная частоту, каков наиболее точный способ оценить соответствующее пиковое значение фундаментальных сигналов?

Одним из способов может быть обнуление сигнала времени для увеличения разрешения БПФ, чтобы центр бина был ближе к расчетной частоте. В этом сценарии один момент, в котором я не уверен, заключается в том, могу ли я заполнять нулями столько, сколько я хочу, или в этом есть некоторые недостатки. Другой - какой центр ячеек я должен выбрать после заполнения нулями как тот, из которого я получаю пиковое значение (потому что каждый может точно не достичь интересующей частоты, даже после заполнения нулями).

Однако мне также интересно, есть ли другой метод, который может дать лучшие результаты, например, оценщик, который использует пиковые значения окружающих двух центров бинов для оценки пикового значения на интересующей частоте.

lR8n6i
источник
2
заполнение нулями до FFT является одним из способов. Другой способ - применить оконную функцию, которая подходит для ваших нужд. Окно с плоским верхом было разработано именно для этой цели. Конечно, если вы уже точно знаете частоту и интересуетесь только одним амплитудом, возможно, существуют более дешевые способы сделать это, чем БПФ.
Sellibitze
1
заполнение нулями не требуется: простая параболическая интерполяция (с 3 точками: imax-1, imax, imax + 1, где imaxпик FFT) даст вам точные результаты
Basj
Убедитесь, что функция интерполяции соответствует оконной функции. Плоская вершина тривиальна, в противном случае вам нужна подходящая пара (например, прямоугольное окно + интерполяция sinc, гауссовское окно + интерполяция
Гаусса
@CedronDawg этот вопрос и его ответы связаны (но не совпадают) с вашей точной формулой частоты. Может быть, вы можете найти это интересным.
Fat32

Ответы:

5

Первый алгоритм, который приходит на ум, - это алгоритм Гёртцела . Этот алгоритм обычно предполагает, что представляющая интерес частота является целым кратным основной частоты. Тем не менее, эта статья применяет (обобщенный) алгоритм к интересующему вас случаю.


Другая проблема заключается в том, что модель сигнала неверна. Это использует 2*%pi*(1:siglen)*(Fc/siglen). Это должно использовать 2*%pi*(0:siglen-1)*(Fc/siglen)для фазы, чтобы выйти правильно.

Я также думаю, что есть проблема с Fc=21.3очень низкой частотой . Низкочастотные реальные значения имеют тенденцию проявлять смещение, когда речь идет о проблемах оценки фазы / частоты.

Я также попытался выполнить грубый сеточный поиск для оценки фазы, и он дает тот же ответ, что и алгоритм Гёртцеля.

Ниже приведен график, который показывает смещение в обеих оценках (Гертцель: синий, Грубый: красный) для двух разных частот: Fc=21.3(сплошная) и Fc=210.3(пунктирная). Как видите, уклон для более высокой частоты намного меньше.

График оси представляет собой начальную фазу, изменяющуюся от 0 до .2 πx2π

введите описание изображения здесь

Питер К.
источник
Только что проверил код алгоритма Гёрзеля на основе статьи. Используя выходное значение DTFT, пик может быть получен очень точно. Тем не менее, есть коэффициент масштабирования ровно 1000. Таким образом, если исходный пик равен 1234, то после Гёрзеля он будет 1234. Кто-нибудь знает, откуда это может появиться?
lR8n6i
Сделал некоторое исследование в то же время. Вероятно, это связано с масштабированием амплитуды: масштабирование амплитуды во временной области = коэффициент частотной области * 2 / N, где N - длина сигнала. Это предположение верно?
lR8n6i
Здравствуй! Я только что узнал, что при использовании алгоритма Гёртцеля амплитуда при результирующем комплексном коэффициенте очень точная, но фаза совершенно неправильная. У кого-нибудь есть идея, откуда это может появиться? Под «фазой» я понимаю фазовую задержку, указанную в основной части исходного сигнала.
lR8n6i
1
@ Rickson1982 Фаза правильная. Вы просто не правильно интерпретируете это. :-) Помните: т.е. он на (90 градусов) превосходит то, что вы ожидаете. π/2sin(ω0t+ϕ)j2[ejϕδ~(ω+ω0+2πk)e+jϕδ~(ωω0+2πk)]π/2
Питер К.
4

Если вы хотите использовать несколько соседних бин FFT, а не только 2, то оконная интерполяция Sinc между результатами комплексного бина может дать очень точную оценку в зависимости от ширины окна.

Оконная интерполяция Sinc обычно встречается в высококачественных аудио апсэмплерах, поэтому статьи по этому вопросу будут иметь подходящие формулы интерполяции с анализом ошибок.

hotpaw2
источник
Спасибо за комментарий. Я тоже попробую этот подход.
lR8n6i
4

Если вы используете Flanagan [1], он вычисляется из разности фаз последовательных фазовых спектров Δϕ (мгновенная частота), и если вы восстанавливаете величину, используя правильный коэффициент (Instantaneous Magnitude) [2], используйте нормализованную функцию sinc: И в конце используйте параболическую интерполяцию вокруг пиковой величины, вы можете получить потрясающие результаты, сегодня я считаю, что это лучший способ, я его использовал, и результаты всегда очень солидно :-)

sin(πx)(πx)

[1] Дж. Л. Фланаган и Р. М. Голден, «Фазовый вокодер», Bell Systems Technical Journal, vol. 45, с. 1493–1509, 1966.

[2] K. Dressler, «Синусоидальное извлечение с использованием эффективной реализации БПФ с множественным разрешением», в Proc. 9th Int. Conf. о цифровых аудиоэффектах (DAFx-06), Монреаль, Канада, сентябрь 2006 г., стр. 247–252.

ederwander
источник
Здравствуй! Большое спасибо за все ваши комментарии. Я расширил свой код (см. Ниже), чтобы объединить фильтр Гертцеля с интерполяцией параболического пика, чтобы получить фазу. Однако результаты все еще не точны (+ - 3-4 градуса). Насколько это близко, или есть ошибки в понимании или кодировании?
lR8n6i
3

Один метод состоит в том, чтобы найти максимум и подогнать к нему параболу, а затем использовать максимум параболы в качестве оценки частоты и величины. Вы можете прочитать все о здесь: https://ccrma.stanford.edu/~jos/sasp/Sinusoidal_Peak_Interpolation.html

HerrLip
источник
3

У меня было много трудностей с этой проблемой пару лет назад.

Я разместил этот вопрос:

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

Я закончил вычисления с нуля и отправил ответ на свой вопрос.

Я удивлен, что мне не удалось найти подобную экспозицию в Интернете.

Я опубликую ответ снова здесь; Обратите внимание, что код разработан для сценария, в котором я перекрываю окно FFT в 4 раза.

π


Эта головоломка требует два ключа, чтобы разблокировать ее.

График 3.3:

введите описание изображения здесь

График 3.4:

введите описание изображения здесь

Код:

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;
}
Число Пи
источник
Вы интерполируете частоту, тогда как OP знает частоту и хочет интерполировать амплитуду.
Финнв
2

Этот код Python даст вам очень точный результат (я использовал его для большого количества музыкальных нот и получил ошибки менее 0,01% полутона) с параболической интерполяцией (метод, успешно используемый McAulay Quatieri, Serra и т. Д. В гармонике + невязка методы разделения)

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
источник
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

Частоты, с которыми вы работаете (частота дискретизации 21,3 Гц при 8 кГц) очень низкие. Поскольку это действительные сигналы, они будут демонстрировать смещение при оценке фазы для ** любой ** частоты.

На этом рисунке показан график зависимости смещения ( phase_est - phase_orig) для Fc = 210.3;(красным) от смещения для Fc = 21.3;. Как видите, смещение для 21.3дела гораздо значительнее .

Другой вариант - уменьшить частоту дискретизации. Зеленая кривая показывает смещение Fs = 800вместо 8000.

введите описание изображения здесь

lR8n6i
источник
1
Спасибо за обновление! Смотрите мой сюжет; Я все еще думаю, что любой фазовый оценщик будет иметь смещение для этой низкой частоты. Один из способов обойти это - использовать известную частоту (если она известна!) Для исправления смещения оценки фазы через справочную таблицу. Но вам нужно быть осторожным: смещение будет меняться с частотой. Еще один способ сделать это будет уменьшить частоту дискретизации.
Питер К.
1
И тебе спасибо! Однако, если вы используете Fs = 8000 Гц и Fc = 210 вместо 210,3, смещение выглядит еще хуже. Есть идеи, откуда это?
lR8n6i
1
Эрк! Без понятия. FWIW, оценщик Geortzel не имеет проблем: goertzel = atan(imag(y(2)),real(y(2)))*180/%pi + 90;. :-) Буду копать еще немного. Смотреть это пространство.
Петр К.
1
Параболическая интерполяция не делает то, что вы думаете, что делает. В частности, если вы замените расчет pс , 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;то вы получите гораздо лучшие ответы --- даже Fc=210. Я совсем не уверен, что текущая версия pдаст вам что-нибудь разумное. Формула интерполяции предназначена для интерполяции амплитуды параболы, но pинтерполирует фазу, которая просто ... странна.
Питер К.
1
Все это нормально, за исключением того, что местоположение пика ( p = (yp1 - ym1)/(2*(2*y0 - yp1 - ym1))) будет неправильным в некоторых случаях, если вы используете ФАЗЫ вместо амплитуд. Это связано с тем, что фазы могут пересекать границу +/- 180 градусов. Все, что нужно, чтобы исправить это для фазы, - это изменить эту строку на мой p2расчет выше.
Питер К.