具体来说，我想生成信号：
$$ y（t）= \ sin（2 \ pi ft）$$
其中$ t = n / f_s $为样本编号$ n $。
一个想法是保持跟踪每个时钟周期旋转一个角度$ \ Delta \ phi = 2 \ pi f / f_s $的向量$（x，y）$。
作为Matlab伪代码实现（实数在C中实现）：
%% Initialization code
f_s = 32768; % sample rate [Hz]
f = 19.875; % some constant frequency [Hz]
v = [1 0]; % initial condition
d_phi = 2*pi * f / f_s; % change in angle per clock cycle
% initialize the rotation matrix (only once):
R = [cos(d_phi), -sin(d_phi) ; ...
sin(d_phi), cos(d_phi)]
然后，在每个时钟周期，我们将向量旋转一点点： />
这样就可以在每个周期仅用4个乘法来计算振荡器。但是，我担心相位误差和幅度稳定性。 （在简单的测试中，我感到惊讶的是振幅没有立即消失或爆炸-也许
sincos指令保证$ \ sin ^ 2 + \ cos ^ 2 = 1 $？）。什么是正确的方法吗？
#1 楼
没错，严格的递归方法很容易随着迭代次数的增加而累积错误。通常要完成的一种更可靠的方法是使用数控振荡器（NCO）。基本上，您有一个可跟踪振荡器瞬时相位的累加器，其更新如下：$$
\ delta = \ frac {2 \ pi f} {f_s}
$$
$$
\ phi [n] =（\ phi [n-1] + \ delta）\ mod 2 \ pi
$$$
然后，每次您都需要将NCO中的累积相位转换为所需的正弦输出。如何执行此操作取决于您对计算复杂性，准确性等的要求。一种显而易见的方法是将输出计算为
$$
x_c [n] = \ cos（\ phi [n]）
$$
$$
x_s [n] = \ sin（\ phi [n]）
$$$
使用现有的任何正弦/余弦实现。在高通量和/或嵌入式系统中，相位到正弦/余弦值的映射通常是通过查找表完成的。查找表的大小（即，您对正弦和余弦的相位参数进行的量化的数量）可以用作内存消耗和近似误差之间的折衷。令人高兴的是，所需的计算量通常与表的大小无关。另外，如果需要，可以利用余弦和正弦函数固有的对称性来限制LUT的大小。您只需要真正存储采样周期的四分之一即可。 （例如，线性或三次插值）。
这种方法的另一个好处是，将频率或相位调制与这种结构结合起来是微不足道的。可以通过相应地改变$ \ delta $来调制输出的频率，并且可以通过直接加到$ \ phi [n] $来实现相位调制。
#2 楼
您所拥有的是一个非常有效的振荡器。潜在的数值漂移问题实际上可以解决。您的状态变量v有两部分，一方面是实际部分，另一方面是虚部。让我们叫r和i。我们知道r ^ 2 + i ^ 2 =1。随着时间的流逝，它可能会上下漂移，但是可以通过乘以增益校正因子（例如$$ g = \ frac {1} {\ sqrt { r ^ {2} + i ^ {2}}} $$显然这很昂贵，但是我们知道增益校正非常接近于单位，我们可以用一个简单的泰勒近似扩展到
$$ g = \ frac {1} {\ sqrt {r ^ {2} + i ^ {2}}} \\约\ frac {1} {2} \ cdot \ left（3-（ r ^ {2} + i ^ {2}）\ right）$$
此外，我们不需要对每个样本都执行此操作，但是每100或1000个样本就足够了保持稳定。如果您执行基于帧的处理，这将特别有用。每帧更新一次就可以了。这是一个快速的Matlab计算10,000,000个样本。
%% seed the oscillator
% set parameters
f0 = single(100); % say 100 Hz
fs = single(44100); % sample rate = 44100;
nf = 1024; % frame size
% initialize phasor and state
ph = single(exp(-j*2*pi*f0/fs));
state = single(1 + 0i); % real part 1, imaginary part 0
% try it
x = zeros(nf,1,'single');
testRuns = 10000;
for k = 1:testRuns
% overall frames
% sample: loop
for i= 1:nf
% phasor multiply
state = state *ph;
% take real part for cosine, or imaginary for sine
x(i) = real(state);
end
% amplitude corrections through a taylor exansion aroud
% abs(state) very close to 1
g = single(.5)*(single(3)-real(state)*real(state)-imag(state)*imag(state) );
state = state*g;
end
fprintf('Deviation from unity amplitude = %f\n',g-1);
评论
$ \ begingroup $
Hilmar在另一个问题中进一步解释了此答案：dsp.stackexchange.com/a/1087/34576
$ \ endgroup $
–西罗林顿
19/12/9在19:47
#3 楼
如果不让它递归更新向量v，就可以避免幅度波动的不稳定。相反，将原型向量v旋转到当前输出相位。这仍然需要一些触发功能，但每个缓冲区仅需要一次。没有幅度漂移和任意频率
伪代码：
如果可以承受量化的频率转换，则可以省去乘法，cexp所需的trig函数以及超过2pi的模数余量。例如fs / 1024（用于1024个采样相量缓冲区）。
评论
$ \ begingroup $
感谢您的回答。 sincos的执行时间与少数乘法相比如何？使用mod操作是否需要注意一些陷阱？
$ \ endgroup $
– nibot
11年8月29日在19:09
$ \ begingroup $
吸引人的是，系统中的所有振荡器都可以使用相同的相幅LUT。
$ \ endgroup $
– nibot
11年8月29日在19:12
$ \ begingroup $
mod 2pi的作用是什么？我也看到过执行mod 1.0的实现。您能扩展一下模运算的作用吗？
$ \ endgroup $
–BigBrownBear00
17-2-22在14:53
$ \ begingroup $
@ BigBrownBear00：模运算使$ \ phi [n] $保持在可管理的范围内。实际上，如果您没有模数，则随着时间的推移，它将增长为非常大的正数或负数（累计相的总量）。这可能是不好的，原因有几个，包括最终溢出或算术精度损失，以及余弦和正弦函数评估的性能降低。如果不需要首先将参数减少到$ [0，2 \ pi）$范围内，则典型的实现会更快。
$ \ endgroup $
–Jason R
17年2月22日在16:11
$ \ begingroup $
$ 2 \ pi $与1.0的因数是实现细节。这取决于平台三角函数的作用域。如果他们期望值在$ [0，1.0）$范围内（即角度以周期测量），则将调整$ \ phi [n] $的方程式以反映该不同的单位。上面答案中的解释假定使用弧度的典型角度单位。
$ \ endgroup $
–Jason R
17-2-22在16:13