该问题可以表述为（非线性）最小二乘问题：

$$F(\phi) = \frac{1}{2}\sum^{n}_{i=1}\left[ A \cdot \sin(\omega i+\phi) - f_{i}(\omega) \right]^{2}$$

其中是相对于最小化的目标函数。$F(\phi)$$\phi$

导数非常简单：

$$F'(\phi)=\sum^{n}_{i=1} A\cdot \cos(\omega i+\phi) \left[A \cdot \sin(\omega i+\phi) - f_{i}(\omega)\right]$$

上述目标函数可以使用梯度下降法（一阶近似）、牛顿法、高斯-牛顿法或Levenberg-Marquardt 法（二阶近似 -需要在这些方法中提供）迭代地最小化。$F''(\phi)$

显然，上述目标函数由于周期性存在多个最小值，因此可以添加一些惩罚项来区分其他最小值（例如，在模型方程但我认为优化只会收敛到最接近的最小值，您可以更新结果减去。$\phi^{2}$$2\pi k, k\in N$

在特定频率下使用 DFT。然后从实部/图像部分计算幅度和相位。它为您提供参考采样时间开始的相位。

在“正常”FFT（或为所有 N 个谐波计算的 DFT）中，您通常使用 f = k*(sample_rate)/N 计算频率，其中 k 是整数。尽管这看起来有点亵渎神明（尤其是对于全整数教会的成员而言），但在进行单个 DFT 时，您实际上可以使用 k 的非整数值。

例如，假设您已经生成（或获得）N = 256 个 27 Hz 正弦波点。（比如说，sample_rate = 200）。256 点 FFT（或 N 点 DFT）的“正常”频率对应于：f = k*(sample_rate)/N = k*(200)/256，其中 k 是整数。但是使用上面列出的参数，34.56 的非整数“k”对应于 27 Hz 的频率。这就像创建一个完全以感兴趣的频率（27 Hz）为中心的 DFT 'bin'。一些 C++ 代码（DevC++ 编译器）可能如下所示：

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <cmath>
using namespace std;

// arguments in main needed for Dev-C++ I/O
int main (int nNumberofArgs, char* pszArgs[ ] ) {
const long N = 256 ;
double sample_rate = 200., amp, phase, t, C, S, twopi = 6.2831853071795865; 
double  r[N] = {0.}, i[N] = {0.}, R = 0., I = 0. ;
long n ;

// k need not be integer
double k = 34.56;

// generate real points
for (n = 0; n < N; n++) {
    t =  n/sample_rate;
    r[n] = 10.*cos(twopi*27.*t - twopi/4.);
}  // end for

// compute one DFT
for (n = 0; n < N; n++) {
    C = cos(twopi*n*k/N); S = sin(twopi*n*k/N);
    R = R + r[n]*C + i[n]*S;
    I = I + i[n]*C - r[n]*S;
} // end for

cout<<"\n\ndft results for N = " << N << "\n";
cout<<"\nindex k     real          imaginary       amplitude         phase\n";

amp = 2*sqrt( (R/N)*(R/N) + (I/N)*(I/N) ) ;
phase = atan2( I, R ) ;
// printed R and I are scaled
printf("%4.2f\t%11.8f\t%11.8f\t%11.8f\t%11.8f\n",k,R/N,I/N,amp,phase);

cout << "\n\n";
system ("PAUSE");
return 0;
} // end main

//**** end program

（PS：我希望上面的内容能很好地转化为stackoverflow——其中一些可能会环绕）

上面的结果是-twopi/4的相位，如生成的实点所示（amp加倍以反映pos/neg频率）。

有几点需要注意——我使用余弦来生成测试波形并解释结果——你必须小心——相位以时间 = 0 为参考，这是你开始采样的时间（即：当你收集 r[0] ），余弦是正确的解释）。

上面的代码既不优雅也不高效（例如：对 sin/cos 值使用查找表等）。

当您使用较大的 N 时，您的结果会变得更准确，并且由于采样率和上面的 N 不是彼此的倍数，因此会有一点误差。

当然，如果要更改采样率 N 或 f，则必须更改代码和 k 的值。您可以在连续频率线上的任何位置插入 DFT 箱 - 只需确保您使用的 k 值对应于感兴趣的频率。

Goertzel 算法有几种不同的公式。提供 2 个状态变量（正交或接近）或复杂状态变量作为可能输出的那些通常可用于参考 Goertzel 窗口中的某个点（例如中间）来计算或估计相位。单独提供单个标量输出的通常不能。

您还需要知道 Goertzel 窗口相对于时间轴的位置。

如果您的信号在 Goertzel 窗口中不完全是整数周期，则窗口中间参考点周围的相位估计可能比参考开始或结束的相位更准确。

如果您知道信号的频率，则完整的 FFT 是多余的。此外，可以将 Goertzel 调谐到 FFT 长度中非周期性的频率，而 FFT 将需要对非周期性窗口频率进行额外的插值或零填充。

复数 Goertzel 相当于使用余弦和正弦基向量或 FFT 旋转因子的递归的 DFT 的 1 个 bin。

如果您的信号是无噪声的，您可以识别两者中的零交叉并确定频率和相对相位。