我正在查看Goertzel算法，但它似乎没有处理相位问题？
#1 楼
在特定频率下使用DFT。然后从实部/成像部分计算振幅和相位。它为您提供了参考采样时间开始的相位。在“正常” FFT（或为所有N个谐波计算的DFT）中，通常使用f = k *（sample_rate ）/ N，其中k是整数。尽管它看起来似乎是牺牲品（特别是对于全整数的教会成员），但是在执行单个DFT时实际上可以使用k的非整数值。例如，假设您已经生成了（或获得）N = 27 Hz的正弦波的256个点。 （假设sample_rate = 200）。 256点FFT（或N点DFT）的“正常”频率将对应于：f = k *（sample_rate）/ N = k *（200）/ 256，其中k是整数。但是使用上面列出的参数，非整数“ k”为34.56将对应于27 Hz的频率。这就像创建一个DFT“箱”，该箱正好位于感兴趣的频率（27 Hz。）上。一些C ++代码（DevC ++编译器）可能如下所示： br />
上面的结果是-twopi / 4的相位，如生成的实点所示（放大器加倍以反映pos / neg频率）。
需要注意的几件事–我使用余弦生成测试波形并解释结果–您必须注意–相位参考时间= 0，这是您开始采样的时间（即：当您收集r [0]时） ），并且余弦是正确的解释。）
上面的代码既不优雅也不高效（例如：使用查找表获取sin / cos值等）。
使用较大的N时，您的结果将更加准确，并且由于上面的采样率和N不是彼此的倍数，因此存在一些误差。
当然，如果要更改采样率N或f，则必须更改代码和k的值。您可以在连续频率线上的任何位置放下DFT仓–只需确保您使用的k值对应于感兴趣的频率即可。
#2 楼
该问题可以表述为（非线性）最小二乘问题：$$ F（\ phi）= \ frac {1} {2} \ sum ^ {n} _ {i = 1} \ left [一个\ 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] $$
可以使用梯度下降法（一阶逼近），牛顿法，高斯-牛顿法或莱芬伯格-马夸特方法（二阶逼近-$ F''（\ phi）$
显然，上述目标函数由于周期性而具有多个最小值，因此可以添加一些惩罚项来区分其他最小值（例如，添加$ \ phi ^ {2 } $到模型方程式）。但是我认为优化将收敛到最近的最小值，您可以更新结果减去$ 2 \ pi k，k \ in N $。
评论
$ \ begingroup $
我不认为您需要因为周期性而受到处罚吗？您可以在收敛的相空间中取最小值，然后做一个模数$ 2 \ pi $，不是吗？
$ \ endgroup $
–太空
2012年9月6日在22:48
$ \ begingroup $
@Mohammad是的，但是某些优化技术可能会使用多个起点，这些起点应收敛到相同的值，或者使用单个全局最小化器假定一个凸函数，而该全局最小化器可以很好地用二次近似。另一个好处是，对于任何起点$ \ phi_ {0} $，我们都以相同的结果结束。
$ \ endgroup $
– Libor
2012年9月6日23:30
$ \ begingroup $
有趣。我是否可以邀请您也解决这个相关问题？ :-)
$ \ endgroup $
–太空
2012年9月7日下午14:16
$ \ begingroup $
@Mohammad好的，我在那里贡献了一点：）
$ \ endgroup $
– Libor
2012年9月7日15:30
$ \ begingroup $
函数fi（w）到哪里去了？ fi（w）不是常数，因此当您获取非常数的导数时，它如何变为零？
$ \ endgroup $
– SamFisher83
2012年9月10日下午16:13
#3 楼
Goertzel算法有几种不同的表达方式。提供2个状态变量（正交或接近）或复杂的状态变量（可能的输出）的变量通常可以用于参考Goertzel窗口中的某个点（例如中间点）来计算或估计相位。仅提供单个标量输出的输出通常不能。您还需要知道Goertzel窗口相对于时间轴的位置。
如果信号在Goertzel窗口中不是精确的整数周期，则相位估计在窗口中间的参考点周围进行定位可能比在相位的起点或终点处进行参考更为准确。
如果您知道信号的频率，则进行完整的FFT是过大的。另外，可以将Goertzel调谐到FFT长度中非周期性的频率，而FFT将需要附加插值或零填充以用于非周期性窗口内频率。
复杂的Goertzel等效于1箱DFT，它对余弦和正弦基向量或FFT旋转因子使用递归。
评论
$ \ begingroup $
相位估计在窗口内的任何地方都不会完全相同，因为您只需在窗口开始时将$ \ omega k $添加到相位估计中，即可计算出样本$ k的相位估计窗口中的$（$ k = 0 $是窗口的开始）？
$ \ endgroup $
–奥利·尼米塔洛（Olli Niemitalo）
19年1月8日在8:30
$ \ begingroup $
否，因为添加wk会在窗口结尾处产生与在非整数周期光圈正弦曲线中开始时不同的相位。但是1 bin DFT在同一点计算单个循环相位。因此，这三个值将全部不同。但是，无论f0是多少，中心相位始终与奇/偶函数的比率有关。
$ \ endgroup $
– hotpaw2
19年1月8日在13:19
$ \ begingroup $
尝试，但我不明白。
$ \ endgroup $
–奥利·尼米塔洛（Olli Niemitalo）
19年1月8日在16:11
$ \ begingroup $
使用余弦（在k = 0时为零相位），略微调整频率（通过微小的无理数，但不更改在k = 0时的相位）。 DFT报告阶段已更改！尝试对正好位于k = N / 2的余弦进行同样的处理。对于任何df，在k = N / 2时都没有变化。罪恶或任何混合都一样。将相位参考点居中显示，随着f0的变化，测得的相位变化较小。例如频率误差不会增加相位测量误差。
$ \ endgroup $
– hotpaw2
19年1月8日在16:56
$ \ begingroup $
是的，如果正弦波和Goertzel滤波器的频率不同，则相位估计误差在窗口中心处较小。在那种情况下，在窗口末端的相位估计值被一个常数所偏置，该常数是窗口中心与末端之间的距离与正弦波和Goertzel滤波器频率之差的乘积。减去该偏差会得到与中心估算值相同的尺寸误差，但它需要知道正弦波的频率。
$ \ endgroup $
–奥利·尼米塔洛（Olli Niemitalo）
19年1月9日在8:25
#4 楼
如果信号无噪声，则可以在两个信号中都识别出零交叉，并确定频率和相对相位。#5 楼
这取决于您对“快速”的定义是什么，您想要的估计准确度，是$ \ phi $还是相对于采样的相位以及函数和参考正弦波上有多少噪声。 />一种方法是仅对$ f（t）$进行FFT，然后查看最接近$ \ omega $的bin。但是，这取决于$ \ omega $接近bin中心频率。
所以：
“快速”是什么意思？
您需要估算的准确度如何？
您想要$ \ phi $（相对于参考的相位）还是相对于采样开始的相位？没关系吗？
每个信号的噪声电平是多少？
PS：我假设您的意思是$ f（t）= A \ sin（ωt+ ϕ）$，而不是比$ f（{\ Huge t}）= A \ sin（ω{\ Huge x} + ϕ）$。
#6 楼
起点：1）将信号和参考正弦波相乘：
$ F（t）$ =A⋅sin（ωt+ ϕ）⋅sin（ωt）=0.5⋅A⋅（cos（ϕ ）-cos（2⋅ωt+ ϕ））
2）找到周期$ T = \ pi / \ omega $：
$ I（\ phi）= \ int_0 ^ TF（t）dt的积分\ =0.5⋅A⋅cos（ϕ）\ cdot T $
3）您可以计算$ \ phi $：
$ cos（\ phi）= I（t）/（0.5 \ cdot A \ cdot T）$
考虑一下：
如何测量A？
如何在$ 0 ..（2 \ cdot \ pi）$区间确定$ \ phi $？ （想想“参考cos波”）
对于离散信号，请改变积分和，并仔细选择T！
#7 楼
您也可以这样做（以数字表示法）：np.arctan( (signal*cos).sum() / (signal*sin).sum() ))
其中信号是您的相移信号，cos和sin是参考信号，并生成一个近似值通过对两个乘积求和在一定时间内积分。
#8 楼
这是对@Kevin McGee的建议的改进，该建议使用具有分数bin索引的单频DFT。凯文（Kevin）的算法效果不佳：虽然在半垃圾箱和整个垃圾箱中非常精确，但也接近整体和一半，也相当不错，但是否则误差可能在5％以内，这对于大多数任务来说可能是不可接受的。我建议通过调整$ N $（即DFT窗口的长度）来改进Kevin的算法，以使$ k $尽可能接近整体。之所以可行，是因为与FFT不同，DFT不需要$ N $是2的幂。
下面的代码在Swift中，但应该直观清楚：
let f = 27.0 // frequency of the sinusoid we are going to generate
let S = 200.0 // sampling rate
let Nmax = 512 // max DFT window length
let twopi = 2 * Double.pi
// First, calculate k for Nmax, and then round it
var k = round(f * Double(Nmax) / S)
// The magic part: recalculate N to make k as close to whole as possible
// We also need to recalculate k once again due to rounding of N. This is important.
let N = Int(k * S / f)
k = f * Double(N) / S
// Generate the sinusoid
var r: [Double] = []
for i in 0..<N {
let t = Double(i) / S
r.append(sin(twopi * f * t))
}
// Compute single-frequency DFT
var R = 0.0, I = 0.0
let twopikn = twopi * k / Double(N)
for i in 0..<N {
let x = Double(i) * twopikn
R += r[i] * cos(x)
I += r[i] * sin(x)
}
R /= Double(N)
I /= Double(N)
let amp = 2 * sqrt(R * R + I * I)
let phase = atan2(I, R) / twopi
print(String(format: "k = %.2f R = %.8f I = %.8f A = %.8f φ/2π = %.8f", k, R, I, amp, phase))
评论
$ \ begingroup $
FFT只是有效地计算DFT的一种方法。对于现代图书馆，两个限制的力量不再存在。如果只需要一个或两个bin值，则最好像您一样直接计算它们。对于单个纯音（实数或复数），仅需要两个bin值即可精确计算频率，相位和幅度。参见dsprelated.com/showarticle/1284.php。数学是相当复杂的，但是有一些链接可以解释派生词。线性代数是真正理解的前提。
$ \ endgroup $
–雪铁龙·道格（Cedron Dawg）
19年11月20日在4:40
评论
$ \ begingroup $
通过调整N使k更接近整体，可以改进此方法。我发布了一个单独的答案，这提高了该算法的准确性。
$ \ endgroup $
– mojuba
19年11月15日在2:48