Извлечение точных частот из FFT Bins с использованием изменения фазы между кадрами
Я просматривал эту фантастическую статью: http://blogs.zynaptiq.com/bernsee/pitch-shifting-using-the-ft/
Будучи фантастическим, он очень тяжелый и тяжелый. Этот материал действительно меня растягивает.
Я извлек математику из модуля кода Стефана, который вычисляет точную частоту для данного бина. Но я не понимаю последний расчет. Может кто-нибудь объяснить мне математическую конструкцию в конце?
Прежде чем копаться в коде, позвольте мне установить сцену:
Допустим, мы установили fftFrameSize = 1024, поэтому имеем дело с 512+1 бинами
Например, идеальная частота Bin[1] соответствует одной волне в кадре. При частоте дискретизации 40 кГц tOneFrame = 1024/40K секунд = 1/40 с, поэтому Bin[1] в идеале собирал бы сигнал 40 Гц.
Установив osamp (overSample) = 4, мы продвигаемся по нашему входному сигналу с шагом 256. Таким образом, первый анализ исследует байты с нуля до 1023, затем с 256 по 1279 и т. Д. Обратите внимание, что каждый float обрабатывается 4 раза.
...
void calcBins(
long fftFrameSize,
long osamp,
float sampleRate,
float * floats,
BIN * bins
)
{
/* initialize our static arrays */
static float gFFTworksp[2*MAX_FRAME_LENGTH];
static float gLastPhase[MAX_FRAME_LENGTH/2+1];
static long gInit = 0;
if (! gInit)
{
memset(gFFTworksp, 0, 2*MAX_FRAME_LENGTH*sizeof(float));
memset(gLastPhase, 0, (MAX_FRAME_LENGTH/2+1)*sizeof(float));
gInit = 1;
}
/* do windowing and re,im interleave */
for (long k = 0; k < fftFrameSize; k++)
{
double window = -.5*cos(2.*M_PI*(double)k/(double)fftFrameSize)+.5;
gFFTworksp[2*k] = floats[k] * window;
printf("sinValue: %f", gFFTworksp[2*k]);
gFFTworksp[2*k+1] = 0.;
}
/* do transform */
smbFft(gFFTworksp, fftFrameSize, -1);
printf("\n");
/* this is the analysis step */
for (long k = 0; k <= fftFrameSize/2; k++)
{
/* de-interlace FFT buffer */
double real = gFFTworksp[2*k];
double imag = gFFTworksp[2*k+1];
/* compute magnitude and phase */
double magn = 2.*sqrt(real*real + imag*imag);
double phase = atan2(imag,real);
/* compute phase difference */
double phaseDiff = phase - gLastPhase[k];
gLastPhase[k] = phase;
/* subtract expected phase difference */
double binPhaseOffset = M_TWOPI * (double)k / (double)osamp;
double deltaPhase = phaseDiff - binPhaseOffset;
/* map delta phase into [-Pi, Pi) interval */
// better, but obfuscatory...
// deltaPhase -= M_TWOPI * floor(deltaPhase / M_TWOPI + .5);
while (deltaPhase >= M_PI)
deltaPhase -= M_TWOPI;
while (deltaPhase < -M_PI)
deltaPhase += M_TWOPI;
(РЕДАКТИРОВАТЬ:) Теперь немного, я не понимаю:
// Get deviation from bin frequency from the +/- Pi interval
// Compute the k-th partials' true frequency
// Start with bin's ideal frequency
double bin0Freq = (double)sampleRate / (double)fftFrameSize;
bins[k].idealFreq = (double)k * bin0Freq;
// Add deltaFreq
double sampleTime = 1. / (double)sampleRate;
double samplesInStep = (double)fftFrameSize / (double)osamp;
double stepTime = sampleTime * samplesInStep;
double deltaTime = stepTime;
// Definition of frequency is rate of change of phase, i.e. f = dϕ/dt
// double deltaPhaseUnit = deltaPhase / M_TWOPI; // range [-.5, .5)
double freqAdjust = (1. / M_TWOPI) * deltaPhase / deltaTime;
// Actual freq <-- WHY ???
bins[k].freq = bins[k].idealFreq + freqAdjust;
}
}
Я просто не вижу этого ясно, хотя кажется, что он смотрит в лицо. Может ли кто-нибудь объяснить этот процесс с нуля, шаг за шагом?
6 ответов
Наконец я понял это; действительно я должен был вывести это с нуля. Я знал, что будет какой-то простой способ получить его, моя (обычная) ошибка заключалась в том, чтобы попытаться следовать логике других людей, а не использовать свой собственный здравый смысл.
Эта головоломка требует два ключа, чтобы разблокировать ее.
Первый ключ заключается в том, чтобы понять, как передискретизация приводит к вращению фазы бина.
Второй ключ взят из графиков 3.3 и 3.4 здесь: http://www.dspdimension.com/admin/pitch-shifting-using-the-ft/
...
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 bin0Freq = (double)sampleRate / (double)fftFrameSize;
bins[k].idealFreq = (double)k * bin0Freq;
// 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;
}
Основной принцип очень прост. Если данный компонент точно соответствует частоте бина, то его фаза не изменится с одного FT на следующий. Однако, если частота не соответствует точно частоте ячейки, то между последовательными FT будет изменение фазы. Частота дельта всего:
delta_freq = delta_phase / delta_time
и уточненная оценка частоты компонента будет тогда:
freq_est = bin_freq + delta_freq
Я сам реализовал этот алгоритм для Performous. Когда вы берете другое БПФ со смещением по времени, вы ожидаете, что фаза изменится в соответствии со смещением, то есть два БПФ, взятые с интервалом 256 выборок, должны иметь разность фаз 256 выборок для всех частот, присутствующих в сигнале (это предполагает, что сами сигналы устойчивы, что является хорошим предположением для коротких периодов, таких как 256 выборок).
Теперь фактические значения фазы, которые вы получаете из БПФ, находятся не в выборках, а в фазовом угле, поэтому они будут отличаться в зависимости от частоты. В следующем коде значение phaseStep - это коэффициент преобразования, необходимый для каждого элемента, то есть для частоты, соответствующей элементу x, сдвиг фазы будет равен x * phaseStep. Для центральных частот ячейки x будет целым числом (номером ячейки), но для фактически обнаруженных частот это может быть любое действительное число.
const double freqPerBin = SAMPLE_RATE / FFT_N;
const double phaseStep = 2.0 * M_PI * FFT_STEP / FFT_N;
Коррекция работает, предполагая, что сигнал в ячейке имеет центральную частоту ячейки, а затем вычисляя ожидаемый сдвиг фазы для этого. Этот ожидаемый сдвиг вычитается из фактического сдвига, оставляя ошибку. Остаток (по модулю 2 pi) берется (диапазон от -pi до pi), и окончательная частота вычисляется с помощью центра бина + коррекция.
// process phase difference
double delta = phase - m_fftLastPhase[k];
m_fftLastPhase[k] = phase;
delta -= k * phaseStep; // subtract expected phase difference
delta = remainder(delta, 2.0 * M_PI); // map delta phase into +/- M_PI interval
delta /= phaseStep; // calculate diff from bin center frequency
double freq = (k + delta) * freqPerBin; // calculate the true frequency
Обратите внимание, что многие смежные элементы часто оказываются скорректированными на одну и ту же частоту, поскольку дельта-коррекция может составлять до 0,5 * FFT_N / FFT_STEP в любом случае, так что чем меньше FFT_STEP, который вы используете, тем дальше будут возможны корректировки (но это увеличивает мощность обработки) нужна, а также неточность из-за неточностей).
Надеюсь, это поможет:)
Это метод оценки частоты, используемый методами фазового вокодера.
Если вы посмотрите на одну точку синусоидальной волны (с фиксированной частотой и с фиксированной амплитудой) во времени, фаза будет продвигаться со временем на величину, пропорциональную частоте. Или вы можете сделать обратное: если вы измеряете, насколько фаза синусоиды изменяется в течение любой единицы времени, вы можете вычислить частоту этой синусоиды.
Фазный вокодер использует два БПФ для оценки фазы со ссылкой на два окна БПФ, а смещение двух БПФ представляет собой расстояние между двумя фазовыми измерениями во времени. Исходя из этого, у вас есть ваша оценка частоты для этого элемента FFT (элемент FFT представляет собой примерно фильтр для изоляции синусоидального компонента или другого достаточно узкополосного сигнала, который помещается в этот элемент).
Чтобы этот метод работал, спектр вблизи используемой ячейки БПФ должен быть достаточно стационарным, например, не изменяться по частоте и т. Д. Это предположение, которое требуется для фазового вокодера.
Может быть, это поможет. Представьте, что ячейки FFT задают маленькие часы или роторы, каждый из которых вращается с частотой ячейки. Для стабильного сигнала (теоретическая) следующая позиция ротора может быть предсказана с использованием математики в бите, который вы не получите. В противовес этому "идеальному" (идеальному) положению вы можете вычислить несколько полезных вещей: (1) разницу с фазой в бине смежного кадра, которая используется фазовым вокодером для лучшей оценки частоты бина, или (2) в более общем смысле отклонение фазы, которое является положительным индикатором появления ноты или какого-либо другого события в аудио.
Частоты сигнала, которые попадают точно в фазу дискретизации, увеличиваются на целое число, кратное 2π. Поскольку фазы бина, которые соответствуют частотам бина, кратны 2π из-за периодической природы БПФ, в этом случае не происходит никакого изменения фазы. Статья, которую вы упоминаете, также объясняет это.