信息处理 - 如何实现数字振荡器？ - 吾爱随笔录

如何实现数字振荡器？

信息处理软件实现振荡器

2022-01-18 00:31:16

我有一个浮点数字信号处理系统，它以固定的采样率运行 $f_s = 32768$ 使用 x86-64 处理器实现的每秒采样数。假设 DSP 系统同步锁定到任何事情，在某个频率下实现数字振荡器的最佳方法是什么 $f$ ?

具体来说，我想生成信号：

y (t) = \sin (2 π f t)

$y(t) = \sin(2\pi f t)$ 在哪里

t = n / f_{s}

$t=n/f_s$ 样品编号

n

$n$ .

一个想法是跟踪向量 $(x,y)$ 我们旋转一个角度 $\Delta\phi = 2\pi f/f_s$ 在每个时钟周期。

作为一个 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)]

然后，在每个时钟周期，我们将向量旋转一点：

%% in-loop code

while (forever),
  v = R*v;        % rotate the vector by d_phi
  y = v(1);       % this is the sine wave we're generating
  output(y);
end

这允许每个周期仅用 4 次乘法来计算振荡器。但是，我会担心相位误差和幅度稳定性。（在简单的测试中，我很惊讶振幅没有立即消失或爆炸——也许sincos指令保证 $\sin^2+\cos^2=1$ ？）。

这样做的正确方法是什么？

3个回答

你是对的，严格递归方法很容易随着迭代次数的增加而累积错误。通常这样做的一种更稳健的方式是使用数控振荡器 (NCO)。基本上，您有一个累加器来跟踪振荡器的瞬时相位，更新如下：

δ = \frac{2 π f}{f_{s}}

$\delta = \frac{2\pi f}{f_s}$

ϕ [n] = (ϕ [n - 1] + δ) \mod 2 π

$\phi[n] = (\phi[n-1] + \delta) \mod 2\pi$

然后，在每个时刻，您都需要将 NCO 中的累积相位转换为所需的正弦输出。如何做到这一点取决于您对计算复杂性、准确性等的要求。一种明显的方法是将输出计算为

x_{c} [n] = \cos (ϕ [n])

$x_c[n] = \cos(\phi[n])$

x_{s} [n] = \sin (ϕ [n])

$x_s[n] = \sin(\phi[n])$

使用您可用的任何正弦/余弦实现。在高吞吐量和/或嵌入式系统中，从相位到正弦/余弦值的映射通常通过查找表完成。查找表的大小（即您对正弦和余弦的相位参数所做的量化量）可以用作内存消耗和近似误差之间的权衡。好消息是所需的计算量通常与表大小无关。此外，如果需要，您可以利用余弦和正弦函数中固有的对称性来限制 LUT 大小；您只需要存储采样正弦曲线周期的四分之一。

如果您需要比合理尺寸的 LUT 更高的精度，那么您始终可以查看表格样本之间的插值（例如线性或三次插值）。

这种方法的另一个好处是，将频率或相位调制与这种结构结合起来很简单。输出的频率可以通过相应地改变来调制，相位调制可以通过简单地直接添加到来实现。 $\delta$ $\phi[n]$

你所拥有的是一个非常好的和高效的振荡器。潜在的数值漂移问题实际上是可以解决的。您的状态变量 v 有两个部分，一个本质上是实部，另一个是虚部。让我们调用 r 和 i。我们知道 r^2+i^2 = 1。随着时间的推移，这可能会上下漂移，但是可以通过乘以像这样的增益校正因子来轻松校正

g = \frac{1}{\sqrt{r^{2} + i^{2}}}

$g = \frac{1}{\sqrt{r^{2}+i^{2}}}$

显然这是非常昂贵的，但是我们知道增益校正非常接近于单位，我们可以通过简单的泰勒展开来近似这个

g = \frac{1}{\sqrt{r^{2} + i^{2}}} \approx \frac{1}{2} \cdot (3 - (r^{2} + i^{2}))

$g = \frac{1}{\sqrt{r^{2}+i^{2}}}\approx \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);

如果不使其递归地更新向量 v，则可以避免不稳定的幅度漂移。相反，将原型向量 v 旋转到当前输出相位。这仍然需要一些三角函数，但每个缓冲区只需要一次。

无幅度漂移和任意频率

伪代码：

init(freq)
  precompute Nphasor samples in phasor
  phase=0

gen(Nsamps)
    done=0
    while done < Nsamps:
       ndo = min(Nsamps -done, Nphasor)
       append to output : multiply buf[done:done+ndo) by cexp( j*phase )
       phase = rem( phase + ndo * 2*pi*freq/fs,2*pi)
       done = done+ndo

如果您可以容忍量化的频率转换，您可以取消乘法、cexp 所需的三角函数以及超过 2pi 的模余数。例如 fs/1024 用于 1024 采样相量缓冲器。

其它你可能感兴趣的问题

上一篇高阶滤波器的级联双二阶部分如何工作？下一篇如何从图像中去除眩光/亮度限制？