0
$\begingroup$

Background

To generate a sine wave the straight forward way is to call the sin method for each data point:

Iteration (pseudo code):

sample = sin(angular_freq * t + phase)
t = t+1

An alternative to calling sin for each iteration is to use complex numbers. Multiplication by a unit complex number is equivalent to rotation in the complex plane.

Initialize:

omega = cos(angular_freq) + i * sin(angular_freq)
phasor = 1 + i*0

Iteration:

phasor = phasor * omega
sample = imag(phasor)

(Note: this method also requires periodic normalization. There is an approximation to scale by (3 - a*a - b*b) / 2 which avoids square roots.)

Question

I've run benchmarks with my code locally and the complex generation method is slightly quicker (4.8 ns/op vs 3.6 ns/op). Probably not enough to matter for my purposes, but it got me curious. If I wanted to generate an FM signal with this efficiently how would that be done? It looks like I'd need to recalculate the rotation complex number (omega) for each frequency change. Is there a way to avoid that recalculation and only need complex math?

$\endgroup$
1
  • $\begingroup$ what language you are using? what is the platform? Running C on cortex-m is different than Python on a mac $\endgroup$ May 25, 2022 at 6:12

1 Answer 1

2
$\begingroup$

First off: Your "slow" pseudocode isn't correct; sin(freq*time + phase) needs to save freq*time+phase for the next iteration, and only change the +phase term over time (as phase=freq/f_sample), if you want to be able to change freq over time. Otherwise changing the frequency will make jumps in your phase, and you can't have that. It doesn't change anyting about the complexity of that approach (you need to do one multiplication of the instantaneous frequency with a constant 1/f_sample, and an addition per output sample).


Is there a way to avoid that recalculation and only need complex math?

Generally, when thinking about this, you'll notice that you really need to compute the sin and cos of a value; and, since standard libraries (of whatever language you're using, unless it's very experimental) will be relatively cleverly implemented, if there was a lower-effort way of doing that through complex math, with the same quality, they'd do it already.

What you can of course do is reduce the precision you calculate $\omega$ with; it's not like having 20 significant digits will do you any good when encoding audio for FM transmission over a noisy narrowband channel! So, if your programming language lets you choose, opt for sincosf instead of sincos (or sinf instead of sin, etc), to use the single-precision instead of double-precision methods.

Note that you're indeed at risk of accumulating errors in the saved phasor (in both the corrected "slow" and the "fast" method), but this accumulation doesn't matter to an DC-blocking audio receiver anyways.

The next step would be reducing the precision of your trigonometric functions further. You don't need high precision, usually, for FM. And even if you do, you might be able to exclude a lot of special cases yourself, so that the computation might be faster.

For example, if you'd implement

initialize:
samp_rate_factor = 1 / samp_rate * modulation_index_normalized
lastphase = 0
freq = get_accumulated_message_samples()

for each natural n:
  delta     = freq[n] * samp_rate_factor
  newphase  = lastphase + delta
  out[n]    = sin(newphase)
  lastphase = newphase

That would constitute a frequency modulation. Now, if you chose samp_rate_factor appropriately, then delta will always be $|\text{delta}|< \pi$. That means that newphase can never differ by more than $\pi$ from lastphase. That means you can do a simple bounds check in every iteration:

for each natural n:
  delta     = freq[n] * samp_rate_factor
  newphase  = lastphase + delta
  if newphase < -pi:
    newphase += 2*pi
  else if newphase > pi:
    newphase -= 2*pi
  out[n]    = sin(newphase)
  lastphase = newphase

Now, that looks like unnecessary additional computation, after all, sin deals with arbitrary real numbers, right? It does, but by internally wrapping them to a $2\pi$-wide interval (or $\pi$, or even $\pi/2$ wide) and then exploiting symmetry and periodicity. That wrapping (modf) is more expensive for "any" possible value than just our two comparisons above. Great, now we have much easier requirements for sin!

Over that interval, a rough polynomial approximation might do. If we choose to approximate $\sin(x)$ using a polynomial $P_N(x) = \sum_{k=0}^{N-1} a_k \cdot x^{2k+1}$ with $N$, exclusively odd, exponents, then we can try how low we can get $N$ before the quality of our FM reduces significantly. A different answer about calculation of trigonometric functions pointed to Olli W's website, which nicely investigates this.

Now, looping back to your original question: What about doing it as a multiplication of complex phasors? Well, that's totally possible, as you demonstrate, and multiplication of two complex numbers is indeed something that needs four multiplications only! It does also have multiplicative error propagation, so my guess is that, since you, no matter what, need some form of trigonometric function, you're better off directly just summing phase and phase difference, and doing one trigonometric function on the result for each sample. Look at the above polynomial:

To calculate the value of $P_N(x)$, you need to do:

  1. calculate $x^1, x^3,\ldots, x^{2(N-1)+1}$, but that's just calculating $x\cdot x$ once, and multiplying $x$ with it $N-1$ times, so, $N$ multiplications in total, and
  2. calculate the product $a_k\cdot x^{2k+1}$, so $N$ multiplications, and, finally.
  3. the $N$-element sum, so $N-1$ additions.

In total, you get ca. $3N$ real-valued operations! Not bad, hm!

So, let's go with $N=4$. The optimal parameter vector $\mathbf\theta=(a_0, a_1, a_2, a_{3})^T$ can be found using numerical optimization quickly¹.

\begin{align} P_4(x) &= a_{ 0 } \cdot x^{ 1 } + a_{ 1 } \cdot x^{ 3 } + a_{ 2 } \cdot x^{ 5 } + a_{ 3 } \cdot x^{ 7 }\\ &= 0.99977253 \cdot x^{ 1 } + -0.16613475 \cdot x^{ 3 } + 0.0080648643 \cdot x^{ 5 } + -0.00015190283 \cdot x^{ 7 } \end{align}


¹ Code to do that, Python3, scipy, numpy, happy fun times:

# Copyright 2022 Marcus Müller
# SPDX-License-Identifier: LGPL-3.0-or-later

import numpy
from scipy import optimize

def odd_polynomial_residual_to_sin(coeff_vec: np_t.ArrayLike) -> float:
  """
  This function calculates the sum of squared differences at 1024
  points on [-\pi;\pi] between the odd-only polynomial using the
  coefficients passed as arguments and the numpy sine.
  We're making sure the approximate math is already "rounded" down
  to float32, so that we don't approximate values we can't really use.
  """
  coeffs = numpy.asarray(coeff_vec, dtype=numpy.float32)
  x = numpy.linspace(-numpy.pi, numpy.pi, 1024, dtype=numpy.float32)
  sine = numpy.sin(x)
  poly_approx = sum((x**(idx*2+1)).astype(numpy.float32) * coeff for idx, coeff in enumerate(coeffs))
  return sum((poly_approx - sine)**2)

# Use [0.5, 0.5, 0.5, 0.5] as initial guess, i.e., try to improve
# the sine approximation that 0.5x¹ + 0.5x³ + 0.5x⁵ + 0.5x⁷ is
result = optimize.minimize(odd_polynomial_residual_to_sin, [0.5]*4)

# print the result
print(" + ".join(f"{coeff:.8g} \\cdot x^{{ {2*k+1} }}" for k, coeff in enumerate(result.x)))
# 0.99977253 \cdot x^{ 1 } + -0.16613475 \cdot x^{ 3 } + 0.0080648643 \cdot x^{ 5 } + -0.00015190283 \cdot x^{ 7 }

The function looks like this:

sine vs polynomial approximation

and the difference between the sin and the polynomial approximation like this:

difference

if you constrain your optimization to the range $[-\pi/2; \pi/2]$ (i.e., change linspace(-pi, pi, …) to linspace(-pi/2, pi/2,…) in odd_polynomial_residual_to_sin), you get much lower error, even, as that region, containing a single critical point, is much "nicer" to approximate with polynomials (note the scaling of the vertical axis!):

difference when optimizing for the small range

You might not even need that kind of precision. So, reducing from 4 to 3 coefficients (replace [0.5] * 4 with [0.5] * 3):

difference: only three coefficients, optimizing for the +- pi/2 range

Note that this sine approximation uses a whopping 6 multiplications and 2 additions, and any reasonably optimizing compiler targeting a CPU that has pipelined arithmetic units and a 4-vectorized addition will gladly make it so that this might even have a throughput of 1 sine every 6 (or even 5) cycles:

// Copyright 2022 Marcus Müller
// SPDX-License-Identifier: LGPL-3.0-or-later
#include <array>

constexpr float sine_approx(float x) {
  const std::array<float, 3> c = {0.99977088, -0.16582583,  0.00757373};
  const float x_sq = x * x;
    float x_pwr = x;
    float sum = x_pwr * c[0];
    x_pwr *= x_sq;
    sum += x_pwr * c[1];
    x_pwr *= x_sq;
    sum += x_pwr * c[2];
    // x_pwr *= x_sq;
    // sum += x_pwr * c[3];
    return sum;
}

In my benchmarking, that's only roughly 1.5 times faster than std::sin(float) (which is libc's sinf under the hood, pretty much); I was expecting more! Also, thanks to SSE, it seems, that the $P_3$ and $P_4$ variants are nearly identically fast. So, you get "precision for free".

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.