P i*_*P i 25 c math signal-processing fft phase
我一直在浏览这篇精彩的文章:http://blogs.zynaptiq.com/bernsee/pitch-shifting-using-the-ft/
虽然太棒了,但它非常艰难而且沉重.这种材料真的让我感到舒服.
我从Stefan的代码模块中提取了数学,该代码模块计算给定bin的确切频率.但我不明白最后的计算.有人能告诉我最后的数学结构吗？
在深入研究代码之前,让我设置一下场景:
假设我们设置fftFrameSize = 1024,所以我们处理512 + 1个bin
例如,Bin [1]的理想频率适合帧中的单个波.在40KHz的采样率下,tOneFrame = 1024/40K秒= 1/40秒,因此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;
}
}
我只是看不清楚,即使它似乎在盯着脸.有人可以一步一步地从头开始解释这个过程吗？
Pau*_*l R 12
基本原理很简单.如果给定的组件与bin频率完全匹配,则其相位将不会从一个FT变为下一个FT.但是,如果频率与bin频率不完全对应,那么在连续的FT之间将存在相位变化.频率增量只是:
delta_freq = delta_phase / delta_time
然后,对组件频率的精确估计将是:
freq_est = bin_freq + delta_freq
Tro*_*nic 11
我已经为Performous自己实现了这个算法.当您在时间偏移处采用另一个FFT时,您希望相位根据偏移而改变,即,相隔256个采样的两个FFT应该对信号中存在的所有频率具有256个采样的相位差(这假设信号本身是稳定的,这对短期如256个样本是一个很好的假设.
现在,您从FFT获得的实际相位值不是采样,而是采用相位角,因此它们将根据频率而不同.在下面的代码中,phaseStep值是每个bin所需的转换因子,即对于与bin 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范围),并用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越小,校正就越远(但这会增加处理能力)由于不准确而需要和不精确).
我希望这有帮助 :)
这是相位声码器方法使用的频率估计技术.
如果你及时观察(固定频率和固定幅度)正弦波上的单个点,相位将随时间推移一个与频率成比例的量.或者您可以进行相反的操作:如果您测量正弦曲线的相位在任何时间单位内的变化程度,您可以计算该正弦波的频率.
相位声码器使用两个FFT来参考两个FFT窗口来估计相位,并且两个FFT的偏移是两个相位测量之间的时间距离.从那时起,您就可以对该FFT区域进行频率估计(FFT区间大致是一个滤波器,用于隔离正弦分量或适合该区间的其他足够窄带信号).
为了使这种方法起作用,使用中的FFT箱附近的频谱必须相当稳定,例如频率不变等.这是相位声码器所需的假设.
最后我想出了这个; 我真的必须从头开始.我知道会有一些简单的方法来推导它,我(通常)的错误是试图遵循其他人的逻辑而不是使用我自己的常识.
这个拼图需要两把钥匙来解锁它.
第一个关键是要了解过采样如何在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;
}