利用帧间相位变化从FFT桶中提取精确频率
利用帧间相位变化从FFT桶中提取精确频率
提问于 2011-01-08 09:03:56
我一直在浏览这篇精彩的文章：http://blogs.zynaptiq.com/bernsee/pitch-shifting-using-the-ft/
虽然很棒，但它是极其艰难和沉重的。这种材料真让我紧张。
我从Stefan的代码模块中提取了数学，该模块计算给定的垃圾桶的确切频率。但我不明白最后的计算。最后有人能给我解释一下数学结构吗？
在深入研究代码之前，让我设置一下场景：
- 假设我们设置了fftFrameSize = 1024，所以我们处理的是512+1回收箱
- 例如，Bin1 1的理想频率适用于帧内的单个波。当采样率为40 40Hz时，tOneFrame = 1024/40K秒=1/40，所以Bin1最好是收集40 so信号。
- 设置osamp (overSample) = 4，我们沿着输入信号前进256个步骤。因此，第一个分析检查字节零到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;
}
}我只是看不清楚，尽管它似乎在盯着我的脸。请有人从零开始，一步一步地解释这个过程好吗？
回答 6
Stack Overflow用户
回答已采纳
发布于 2011-02-08 09:10:35
最后，我发现了这一点；实际上，我不得不从零开始推导它。我知道会有一些简单的方法来推导它，我(通常)的错误是试图遵循别人的逻辑，而不是使用我自己的常识。
这个难题需要两个键才能解锁。
- 第一个关键是了解过采样是如何在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;
}Stack Overflow用户
发布于 2011-01-08 09:19:27
基本原理很简单。如果一个给定的分量完全匹配一个bin频率，那么它的相位就不会从一个FT变化到另一个FT。但是，如果频率与bin频率不完全对应，则在连续FTs之间会发生相位变化。频率增量只是：
delta_freq = delta_phase / delta_time然后，对该成分的频率的精确估计将是：
freq_est = bin_freq + delta_freqStack Overflow用户
发布于 2011-02-06 11:14:33
我已经为表现良好自己实现了这个算法。当你一次采取另一个FFT偏移，你期望相位变化根据偏移，即两个FFT采取的256个样本应有一个相位差为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;校正的工作方式是假设一个垃圾箱中的信号具有本中心频率，然后计算该信号的期望相移。这个预期的移位是从实际的移位中减去的，留下了错误。取余数(模2π) (-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请注意，许多相邻的回收箱通常被更正为相同的频率，因为增量校正可以达到0.5 * FFT_N / FFT_STEP桶，因此您使用的FFT_STEP越小，就越有可能进行更正(但这会增加所需的处理能力，以及由于不准确而导致的不精确)。
我希望这有帮助:)
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/4633203
复制相关文章
相似问题
腾讯云开发者
