# Hilbert Transform in C provides possibly strange results

I am calculating some signals that depend on the Hilbert Transform and, after following various guides online and SO, my function looks like this:

void Hilbert_Transf(COMPLEX *cplx_data, const int LENGTH)
{
float Imag_Val = 0.0f;
for(int i = 0; i < LENGTH; i++)
{
Imag_Val = cplx_data[i].y;
cplx_data[i].y = (Imag_Val > 0 ? -PI / 2 : (Imag_Val < 0 ? PI / 2 : 0));
}
}


It is essentially testing if the imaginary value is positive, so it is replaced by -PI / 2, or by PI / 2 if it is negative, and by 0 if 0. To test the function I used a sine wave, computed by:

const int LENGTH = 1024;
float input[LENGTH];
for(int i = 0; i < LENGTH; i++)
input[i] = sin(i * 2 * PI * 0.01);


The FFT process is throught Intel MKL, so in the function above I just considered COMPLEX *cplx_data to be a complex in the form MKL takes and the member .y being its imaginary part. Then I write to files the input data (sine wave) and output from the HT, and plot them. So you can see the input sine in blue and the HT of it in red, as well as the "serrated" peaks, troughs. I wasn't expecting this "stepped" behavior. Can you guys tell me what is wrong/missing and how it should be? If you need more information, just let me know.

========== EDIT ===========

Matt L., following your script, here is the result of a plot which looks more decent: Blue is the input and red, the output. I will of course find a way to reproduce this in C, but I will ask a couple more questions on the comments so we can close it.

• Whatever you are computing here, it's not a Hilbert Transform. Recommendation: start with the math, make sure that this is correct (and exactly what you want), than implement the code, and then verify the code against the math. – Hilmar Jan 20 '19 at 15:39
• yeah, what Hilmar said. – robert bristow-johnson Jan 21 '19 at 11:30

The function in your question doesn't compute the Hilbert transform.

The Hilbert transform can be implemented either in the time domain by filtering the input signal with a filter kernel that approximates the ideal kernel of the Hilbert transform (there are approximations for FIR and IIR Hilbert transformers), or it can be implemented in the frequency domain.

An example of a frequency domain implementation is the Matlab/Octave function hilbert. But note that this function actually computes the analytic signal, so the Hilbert transform is the imaginary part of the result.

A simple direct frequency domain implementation of the Hilbert transform is shown in the Matlab/Octave script below. Note that it only uses real (I)FFTs.

% compute Hilbert transform of vector x

N = length(x);
X = fft(x);

if ~rem(N,2),       % N is even

n = N/2;

% set values at DC and Nyquist to zero
X(1) = 0; X(n+1) = 0;

% multiply with -1i * sign(w)
X(2:n) = -1i * X(2:n);
X(n+2:N) = 1i * X(n+2:N);

else                % N is odd

n = (N+1)/2;

% set value at DC to zero
X(1) = 0;

% multiply with -1i * sign(w)
X(2:n) = -1i * X(2:n);
X(n+1:N) = 1i * X(n+1:N);

endif

y = real(ifft(X));

• @Matt L. Thanks for your reply. You can checking here: link , at minute 6:00. But if you see any mistake in the explanation, I am all ears. As I mentioned, I am using C, not Matlab, and I intend to apply Euler's formula to rotate the phase of a signal by X degrees once I get the HT right. – JayY Jan 20 '19 at 16:29
• @JayY: He's talking about a phase shift of $-\pi/2$, not about replacing the imaginary part by $\pm\pi/2$. I know you're using C, but the reason why I linked to Matlab's hilbert function is because in the documentation the algorithm is explained, and it's pretty straightforward. – Matt L. Jan 20 '19 at 18:20
• @Matt L. Again, thanks for your comment. I have inspected Matlab/Octave's Hilbert script and reproduced it in C (it zeroes half of the real/imag arrays, and multiplies by 2 the other half). When I apply the inverse FFT to have the real signal back so I can proceed to Euler's formula, it really just changes the magnitude. As you said, the HT is the imaginary part of the result, but (FFT_Length == Length / 2 + 1) in the FFT packages I use, so I can't multiply without doing iFFT. What do you suggest? – JayY Jan 20 '19 at 18:51
• @JayY: I added a short code fragment to my answer explaining how the Hilbert transform can be implemented using real FFTs. – Matt L. Jan 21 '19 at 11:30
• @Matt L. Thanks for the sample and please see above the result. I'd like to ask just for a couple more clarifications: 1 - why do you set the DC and the Nyquist freq to 0? 2 - This sample code requires N to be even, but I don't know beforehand if my input will be even. In case it is odd, then it is when these "limits" calculated by hilbert.m come into play (not demonstrated here)? – JayY Jan 21 '19 at 12:09

I am not a signal processing expert, but I have made it a good practice not to mix concepts. There is something called the Hilbert Transform, and there is an analytic signal.

Here is what I do to compute the Hilbert Transform.

• Starting with amplitude sampled values, equally spaced samples, the time-domain signal is placed in an array of $$N$$ locations of type float.
• Take the FFT of the time-domain array, and the result is a frequency domain array of $$N$$ points of type complex.
• Visit every bin location from $$0$$ to $$N-1$$, copy the real value to a temp, copy the imag value to the real value, then copy the temp to the imag value, and finally multiply the real value part by $$-1$$.
• Do this for every bin location (no shortcuts).
• Take the inverse FFT of the frequency domain array, and the result is the original time-domain signal shifted in phase by 90 degrees.

The phase-shifted signal can then be used for signal processing.

Here is what I do to form the analytic signal:

• Starting with amplitude sampled values, equally spaced samples, the time domain signal is placed in an array of $$N$$ locations of type float.
• Take the FFT of the time-domain array, and the result is a frequency domain array of $$N$$ points of type complex.
• Zero the value in real and imag at location bin.
• Visit every bin location from $$1$$ to $$N/2$$, where bin $$N/2$$ is the Nyquist frequency), and multiply the real value by $$2$$ and multiply the imag value by $$2$$.
• Now visit every bin location from $$N/2 + 1$$ to $$N-1$$ (the last location of the array) and zero the real and imag values of each bin location; that means "zero the negative frequencies" of the FFT frequency domain array.
• Take the inverse FFT of the frequency domain array, and the result is the analytic signal that can be used for signal processing.