使用帧之间的相位变化从 FFT Bins 中提取精确频率
Posted
技术标签:
【中文标题】使用帧之间的相位变化从 FFT Bins 中提取精确频率【英文标题】:Extracting precise frequencies from FFT Bins using phase change between frames 【发布时间】:2011-06-05 17:09:03 【问题描述】:我一直在浏览这篇精彩的文章:http://blogs.zynaptiq.com/bernsee/pitch-shifting-using-the-ft/
虽然很棒,但它非常艰难和沉重。这种材料真的让我很紧张。
我从 Stefan 的代码模块中提取了数学,该模块计算给定 bin 的确切频率。但我不明白最后的计算。有人能解释一下最后的数学结构吗?
在深入研究代码之前,让我先设置一下场景:
假设我们设置 fftFrameSize = 1024,所以我们要处理 512+1 个 bin
例如,Bin[1] 的理想频率适合帧中的单个波。在 40KHz 的采样率下,tOneFrame = 1024/40K 秒 = 1/40s,因此 Bin[1] 理想情况下会收集 40Hz 的信号。
设置 osamp (overSample) = 4,我们以 256 步长沿输入信号前进。因此,第一次分析检查字节 0 到 1023,然后检查字节 256 到 1279,依此类推。注意每个浮点数被处理 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;
我只是看不清楚,尽管它似乎在盯着脸看。有人可以从头开始一步一步地解释这个过程吗?
【问题讨论】:
如何获得BIN * bins
它代表什么?
【参考方案1】:
基本原理很简单。如果给定的组件与 bin 频率完全匹配,则其相位不会从一个 FT 变为下一个 FT。然而,如果频率不完全对应于 bin 频率,则在连续 FT 之间会有相位变化。频率增量只是:
delta_freq = delta_phase / delta_time
然后对组件频率的精确估计将是:
freq_est = bin_freq + delta_freq
【讨论】:
抱歉我很笨,但我还是不明白为什么这是真的。使用这个数学我仍然觉得很不靠谱。 如果 2 个 FFT 的偏移量与正弦波的一个周期不同,那么即使正弦波频率以 bin 为中心,也会发生相位变化。 知道频率的一个定义是相位变化率,即f = dϕ/dt
。
我会冒险有人嫉妒你的 l33tDSPsk1llz :p 好吧,不是我。我非常感谢您和 HotPaw 提供了全新的视角。现在我真的可以理解这一点了——终于!!!
@Ohmu:很高兴听到你取得了进展——如果你打算做更多这类事情,我建议阅读一本好的介绍性 DSP 书——Richard Lyons 的书,了解数字信号处理,非常好,比大多数人实用得多。【参考方案2】:
这是相位声码器方法使用的频率估计技术。
如果您及时查看(固定频率和固定幅度)正弦波上的单个点,相位将随时间提前与频率成正比的量。或者你也可以反过来:如果你测量一个正弦曲线的相位在任何单位时间内变化了多少,你就可以计算出那个正弦曲线的频率。
相位声码器使用两个 FFT 参考两个 FFT 窗口估计相位,两个 FFT 的偏移量是 2 个相位测量值之间的时间距离。从那时起,您就有了对该 FFT bin 的频率估计(FFT bin 大致是一个滤波器,用于隔离正弦分量或其他适合该 bin 的足够窄带信号)。
要使此方法起作用,使用的 FFT bin 附近的频谱必须相当稳定,例如频率不变等。这是相位声码器需要的假设。
【讨论】:
【参考方案3】:恰好落在 bin 频率上的信号频率将 bin 相位提前 2π 的整数倍。由于 FFT 的周期性,对应于 bin 频率的 bin 相位是 2π 的倍数,因此在这种情况下没有相位变化。你提到的文章也解释了这一点。
【讨论】:
如果 FFT 步长与 FFT 大小相同,那将是正确的。然而,这里的步长变得更小(osamp 因子),然后即使对于中心频率,相位也不再保持不变。例如。仅考虑一个样本的 FFT 步骤。对于较低频率,基本上不会有相移,而对于非常高的频率,可能会有高达 PI 的相位差。 我已经回答了我自己的问题。但是,如果我对我的答案给予赏金,它将丢失。由于他很棒的开源项目(Performous),我打算把它给 Tronic,但他有很多积分!所以......享受;)【参考方案4】:我自己为Performous 实现了这个算法。当您在一个时间偏移处进行另一个 FFT 时,您希望相位根据偏移量而变化,即相隔 256 个样本的两个 FFT 对于信号中存在的所有频率应该具有 256 个样本的相位差(这假设信号本身是稳定的,这对于像 256 个样本这样的短期内是一个很好的假设)。
现在,您从 FFT 获得的实际相位值不是样本,而是相位角,因此它们会因频率而异。在下面的代码中,phaseStep 值是每个 bin 所需的转换因子,即对于 bin x 对应的频率,相移将为 x * phaseStep。对于 bin 中心频率,x 将是一个整数(bin 编号),但对于实际检测到的频率,它可以是任何实数。
const double freqPerBin = SAMPLE_RATE / FFT_N;
const double phaseStep = 2.0 * M_PI * FFT_STEP / FFT_N;
通过假设 bin 中的信号具有 bin 中心频率,然后计算预期相移来进行校正。从实际班次中减去这个预期班次,留下误差。取余数(模 2 pi)(-pi 到 pi 范围)并使用 bin 中心 + 校正计算最终频率。
// 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
请注意,许多相邻的 bin 通常最终被校正为相同的频率,因为 delta 校正可以高达 0.5 * FFT_N / FFT_STEP bin所需的处理能力以及由于不准确而导致的不精确性)。
我希望这会有所帮助:)
【讨论】:
我现在有一些“论文风格”的基本原理可供参考。但我不够聪明,无法从这些解释中自己制定数学。我经过一些解释,逐行生成数学。数学证明。 也许这会有所帮助? sengpielaudio.com/calculator-timedelayphase.htm(时间延迟以毫秒为单位,但我想您可以将 256 个样本转换为适当的时间)【参考方案5】:我终于想通了;真的我不得不从头开始。我知道会有一些简单的方法来推导它,我的(通常)错误是试图遵循别人的逻辑而不是使用我自己的常识。
这个谜题需要两把钥匙来解锁。
第一个关键是了解过采样如何在 bin 相位上引入旋转。
第二个键来自此处的图 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;
【讨论】:
编辑:查看我在math.stackexchange.com/questions/9416/… 的回答以了解垃圾箱轮换【参考方案6】:也许这会有所帮助。将 FFT 箱视为指定小时钟或转子,每个都以箱的频率旋转。对于稳定的信号,转子的(理论)下一个位置可以使用您没有得到的位中的数学来预测。 针对这个“应该”(理想)位置,您可以计算几个有用的东西:(1) 与相邻帧的 bin 中的相位差,phase vocoder 使用它来更好地bin 频率的估计,或者 (2) 更一般地说 相位偏差,这是音频中音符开始或其他事件的积极指标。
【讨论】:
以上是关于使用帧之间的相位变化从 FFT Bins 中提取精确频率的主要内容,如果未能解决你的问题,请参考以下文章