fftw computing analytic signal from r2c (real to complex) data

Question

I have a real signal which I want to fourier transform, compute the analytic signal in the frequency domain and then transform back. The problem with

fftw_plan fftw_plan_dft_r2c_1d(int n, double *in, fftw_complex *out,
                           unsigned flags);

is that it only has the positive frequency components, while I need to zero the negative frequency components (fftw_complex(f)=0 for f < 0).

The obvious solution would be to zero-pad fftw_complex and transform back with complex to complex, however, for performance/ memory reasons I would like to avoid this.

@Rearden The values of the negative frequency components are the conjugate complex values of the positive frequency components. This means: If you have computed: f(1) = a_1 + j b_1, f(2) = a_2 +j b_2 , ... Then you can construct the values for the negative frequencies by: f(-1) = a_1 - j b_1, f(-2) = a_2 - j b_2. — Heinz M.
The solution is to do what you want to avoid, use a full length complex IFFT. — hotpaw2

francis francis · Accepted Answer · 2017-09-20T22:18:03

So the aim is to optimize the computation of the analytic signal corresponding to a real signal u(). As you wrote, the obvious solution consists in:

Compute the forward real-to-complex DFT using r2c FFTW function. The Fourier coefficients corresponding to negative frequencies are not computed, as they are complex conjugates of the Fourier coefficients corresponding to positive frequencies.
Pad the complex array with zeros, so as to zero all negative frequencies.
Compute the complex-to-complex backward DFT to obtain the complex analytical signal using c2c FFTW.

First, since the analytical signal of u() is complex in the time domain, it is twice as big as u() and saving memory is going to be hard. Since the real part of the analytic signal is already known (that is u() itself), computing the analytical signal is rigorously equivalent to applying the Hilbert tranform to u(), since its result is the imaginary part of the analytical signal.

To apply the Hilbert transform using FFTW, please refer to the algorithms in [1]:

Apply r2c FFTW transform.
For all (that is positive) frequencies, switch the real part and the complex part. Take the opposite while switching. This boils down to multiplying each complex coefficient by -j, where j^2=-1.
Apply c2r FFTW transform. The resulting real array is the imaginary part of the analytical signal.

As a result, the imaginary part of the analytical signal can be computed by using a c2r transform instead of a c2c transform. Nevertheless, the imaginary part is obtained as a second array: the arrays must be zipped shut to obtain a result perfectly equivalent to the first method. Using advanced real-data dfts and strides of 2 could help removing that operation.

Even further,still in [1], it can be noticed that multiplying by -j.sign(w) in the frequency domain corresponds to a convolution in the time domain (See Discrete Hilbert transform). The discrete convolution kernel h is the backward DFT of -j.sign(w). As a result, it writes:

$h_i=\sum_{k=0}^{n-1}-j.sign(w_k).e^{2\pi \frac{k.i}{n} j}$

If n is even, indexes k>0 strictly below n/2 correspond to positive frequencies, index n/2 must be silenced as it corresponds to both n/2 and -n/2 frequencies and indexes higher than n/2 correspond to negative frequencies. Index 0 can be silenced. As a result, the sum writes:

$h_i=\sum_{k=1}^{n/2-1}-j.e^{2\pi \frac{k.i}{n} j}+\sum_{k=1}^{n/2-1}j.e^{-2\pi \frac{k.i}{n} j}$

The geometric sequences are easily summed and the magic operates: all even terms of the convolution kernel h are null!

i%2==0 => h_i=0

(On 2017/09/21, wikipedia is wrong regarding the value of h_i for odd i, see Todoran et. al. [1])

As a consequence even terms of u() are never mixed with odd terms of u() as the convolution with h is performed. Indeed, odd terms of h(u)=u*h depend on even terms of u() and even terms of h(u) depend on odd terms of u(). This seems a promicing track but let's quote Todoran et. al. [1] about it:

This algorithm seems to be computed in a shorter time. In fact, it requires a longer time than the algorithms that use the discrete Fourier transform. This is explained by the fact that for DFT were developed fast calculus algorithms (FFT -Fast Fourier Transform).

Short conclusion: FFTW wins again...

[1] Gheorghe TODORAN, Rodica HOLONEC and Ciprian IAKAB, Discrete Hilbert Transform. Numerical Algorithm.ACTA ELECTROTEHNICA, 2008, 49, 4, 485-490

fftw computing analytic signal from r2c (real to complex) data

1 Answers