How can I effectively calculate the phase angle of a complex number that is (essentially) equal to zero?

Question

I'm writing a C++ program that takes the FFT of a real input signal containing double values and returns a vector X containing std::complex<double> values. Once I have the vector of results I then attempt to calculate the magnitude and phase of the result.

I am running into an issue with calculating the phase angle when one of the outputs is "zero". Zero is in quotes because when a calculation that results in 0 returns a double, the returned value will be very near zero, but not quite exactly zero.

For example, at index 3 my output array has the calculated "zero" value:

X[3] = 3.0531133177191805e-16 - i*5.5511151231257827e-17

I am trying to use the standard library std::arg function that is supposed to return the phase angle of a complex number. std::arg(X[3])

While X[3] is essentially 0, it is not EXACTLY 0 and the way phase is calculated this causes a problem because the calculation uses the ratio of the imaginary part divided by the ratio of the real part which is far from 0!

Doing the actual calculation results in a far from desirable result.

How can I make C++ realize that the result is really 0 so I can get the correct phase angle?

I'm looking for a more elegant solution than using an arbitrary hard-coded "epsilon" value to compare the double to, but so far searching online I haven't had any luck coming up with something better.

Answer 1

It was brought to my attention that my original answer will not work when the input to the FFT has been scaled. I believe I have an actual valid solution now... The original answer is kept below so that the comments still make sense.

From the comments on this answer and others, I've gathered that calculating the exact rounding error in the language may technically be possible but it is definitely not practical. The best practical solution seems to be to allow the user to provide their own noise threshold (in dB) and ignore any data points whose power level falls below that threshold. It would be impossible to come up with a generic threshold for all situations, but the user can provide a reasonable threshold based on the signal-to-noise ratio of the signal being analyzed and pass that in.

A generic phase calculation function is shown below that calculates the phase angles for a vector of complex data points.

std::vector<double> Phase(std::vector<std::complex<double>> X, double threshold, double amplitude)
{
    size_t N = X.size();
    std::vector<double> X_phase(N);

    std::transform(X.begin(), X.end(), X_phase.begin(), [threshold, amplitude](const std::complex<double>& value) {
        double level = 10.0 * std::log10(std::norm(value) / std::pow(amplitude, 2.0));
        return level > threshold ? std::arg(value) : 0.0;
    });

    return X_phase;
}

This function takes 3 arguments:

1. The vector of complex signal data you want to calculate the phase of.
2. A sensible threshold -- Can be calculated from the signal-to-noise ratio of whatever measurement device was used to capture the signal. If your signal contains no noise other than the rounding errors of the language itself you can set this to some arbitrary really low value, like -120dB.
3. The maximum possible amplitude of your input signal. If your signal is calculated, this should simply be set to the amplitude of your signal. If your signal is measured, this should be set to the maximum amplitude the measuring device is capable of measuring (If your signal comes from reading an audio file, often its data will be normalized between -1.0 and 1.0. In this case you would just set the amplitude value to 1.0).

This new implementation still provides me with the correct results, but is much more robust. By leaving the threshold calculation to the user they can set the most sensible value themselves based on the characteristics of the measurement device used to measure their input signal.

Please let me know if you notice any errors or any ways I can further improve the design!

---

Original Answer

I found a solution that seems generic enough.

In the #include <limits> header, there is a constant value for std::numeric_limits<double>::digits10.

According the the documentation:

The value of std::numeric_limits<T>::digits10 is the number of base-10 digits that can be represented by the type T without change, that is, any number with this many significant decimal digits can be converted to a value of type T and back to decimal form, without change due to rounding or overflow.

Using this I can filter out any output values that have a magnitude lower than this limit:

Calculate the phase of X[3]:

int N = X.size();

auto tmp = std::abs(X[3])/N > std::pow(10, -std::numeric_limits<double>::digits10)
         ? value
         : 0.0

double phase = std::arg(tmp);

This effectively filters out any values that are not precisely zero due to rounding errors within the C++ language itself. It will NOT however filter out garbage data caused by noise in the input signal.

After adding this to my phase calculation I get the expected results.