4

First, please excuse my ignorance in this field, I'm a programmer by trade but have been stuck in a situation a little beyond my expertise (in math and signals processing).

I have a Matlab script that I need to port to a C++ program (without compiling the matlab code into a DLL). It uses the hilbert() function with one argument. I'm trying to find a way to implement the same thing in C++ (i.e. have a function that also takes only one argument, and returns the same values).

I have read up on ways of using FFT and IFFT to build it, but can't seem to get anything as simple as the Matlab version. The main thing is that I need it to work on a 128*2000 matrix, and nothing I've found in my search has showed me how to do that.

I would be OK with either a complex value returned, or just the absolute value. The simpler it is to integrate into the code, the better.

Thank you.

Hristo Iliev
  • 72,659
  • 12
  • 135
  • 186
Jordan
  • 3,998
  • 9
  • 45
  • 81
  • Although technically a programming question, this one is better suited for the [DSP](http://dsp.stackexchange.com/) site. – Hristo Iliev Aug 13 '12 at 21:58

4 Answers4

9

The MatLab function hilbert() does actually not compute the Hilbert transform directly but instead it computes the analytical signal, which is the thing one needs in most cases. It does it by taking the FFT, deleting the negative frequencies (setting the upper half of the array to zero) and applying the inverse FFT. It would be straight forward in C/C++ (three lines of code) if you've got a decent FFT implementation.

GregC
  • 7,737
  • 2
  • 53
  • 67
André Bergner
  • 1,429
  • 10
  • 10
  • Yes, I found the algorithm on the Matlab page (http://www.mathworks.com/help/toolbox/signal/ref/hilbert.html). I installed the FFTW library (http://www.fftw.org/) and then just followed the algorithm steps. FFTW's scaling on the inverse transforms gave me a bit of confusion, but it all works now. – Jordan Aug 14 '12 at 20:16
  • I have it working perfectly, but the problem is that it's SLOW. I mean, very slow, slower than MATLAB was. I'm using FFTW's advanced interface to do many transforms with the same plan, but it's taking 20 seconds to do 9088 hilbert transforms (the three steps) of 2078 data points each. MATLAB, on the same computer, takes 6 seconds. Any ideas how to do it faster? – Jordan Aug 15 '12 at 18:07
  • 1
    The FFT (FAST Fourier transform) factorizes the data set into smaller sets (of prime number length) and performs a DFT on each of those. It is only fast if it can factorize the data as much as possible, ideally powers of 2 (the smalles prime number). You need to zero pad your data to a length where it can be factorized easily (take 4096 for instance). Your number 2078 factors into 2 * 1039, which is responsible for the poor performance. – André Bergner Aug 15 '12 at 20:16
  • Thank you, this is helpful. I changed my number of samples to 2048, and achieved a 40% time decrease. Unfortunately, it's still significantly slower than Matlab. Any other cool tips? – Jordan Aug 16 '12 at 20:59
  • Also, FYI, I found that it was much faster to go with a base prime of 3 for 2187 (only 109 padded zeros) then to go with a base of 2 for 4096 (2018 padded zeros). – Jordan Aug 16 '12 at 22:01
  • Nice, that's interesting. Btw, most people believe it only works with base 2 and must be that, but the secret are the small primes. – André Bergner Aug 17 '12 at 06:47
2

This looks pretty good, as long as you can deal with the GPL license. Part of a much larger numerical computing resource.

Matt Phillips
  • 9,465
  • 8
  • 44
  • 75
  • I had a look at that code and cannot recommend it. It computes the hilbert transform by direct convolution. There are two problems to that approach. Firstly, a direct convolution is very slow opposed to a FFT based fast convolution, especially since the kernel of the Hilbert transform is decaying very slowly. Secondly, the HT kernel has a singulary at zero and the numerical computation of the convolution integral is ill-posed. – André Bergner Aug 17 '12 at 06:55
  • @AndréBergner Thank you for your input. But as for speed, what about the time taken to perform the FFT, IFFTs required for the frequency domain convolution? Wouldn't the time *in total* still be larger this way? – Matt Phillips Aug 17 '12 at 13:48
  • No, this technique is called fast convolution and very popular in audio processing, for instance to apply a long reverb to music in real time. And, as I wrote, there is also the numeric problem with the singularity. In that linked software they avoid it, by moving the kernel a bit, so that the singularity is between two samples. But that adds a fractional delay to the output, which is probably not desired as it distorts the result when constructing the analytic signal. Btw, there is also another fast approximation method based on fast recursive (IIR) filters – so called Hilbert transformers. – André Bergner Aug 18 '12 at 16:21
0

Simple code below. (Note: this was part of a bigger project). The value for L is based on the your determination of your order, N. With N = 2L-1. Round N to an odd number. xbar below is based on the signal you define as the input to your designed system. This was implemented in MATLAB.

L = 40;
n = -L:L; % index n from [-40,-39,....,-1,0,1,...,39,40];
h = (1 - (-1).^n)./(pi*n); %impulse response of Hilbert Transform
h(41) = 0; %Corresponds to the 0/0 term (for 41st term, 0, in n vector above)

xhat = conv(h,xbar); %resultant from Hilbert Transform H(w);

plot(abs(xhat))
user3085701
  • 1
0

Not a true answer to your question but maybe a way of making you sleep better. I believe that you won't be able to be much faster than Matlab in the particular case of what is basically ffts on a matrix. That is where Matlab excels!

Matlab FFTs are computed using FFTW, the de-facto fastest FFT algorithm written in C which seem to be also parallelized by Matlab. On top of that, quoting from http://www.mathworks.com/help/matlab/ref/fftw.html:

For FFT dimensions that are powers of 2, between 214 and 222, MATLAB software uses special preloaded information in its internal database to optimize the FFT computation.

So don't feel bad if your code is slightly slower...

Johannes Rebling
  • 11
  • 1
  • 2