Hilbert Transform in C

[cristianoych]

Hello guys,

Anyone know about the Hilbert Transform or Complex envelope theory? Currently, I am doing a project regarding this function to process the Real time input signal frequency. The function has to be written in C. I had done some research and found that Hilbert Transform can be implemented by FFT => Phase Shift => IFFT method. Does anyone know how to write this method in C ? I am very appreciate from you helps. Thanks.

[jemidiah]

I'm sorry, I don't have enough of a background in signal processing to be of much use here. I'd be surprised if anyone frequenting this forum does either, though I suppose it's possible. Have you looked for existing code online? Of course you could expand your search from just C code to C++ or other languages if you're willing to translate.

There seems to be quite a bit of material on the Hilbert Transform, like this, this, and this. These are random links that looked useful in a brief search. I apologize if they're not helpful, though the first one in particular seems to be highly relevant.

Perhaps you'd be best served finding a signal processing forum or asking an expert (say, a professor you're even vaguely acquainted with) directly.

[zaza]

Aren't there loads of methods on the web for implementing Fast Fourier Transforms and their inverse? And a phase shift is just a phase shift.....

I don't know much about the Hilbert Transform, but if it's as you suggest then I would have thought there would be bags of code out there. Have you looked in Numerical Recipes?

[jemidiah]

I had done some research and found that Hilbert Transform can be implemented by FFT => Phase Shift => IFFT method.

I should mention equation (3) from my first link gives the Hilbert transform (in a technically special case the OP likely doesn't care a whit about) as

H{f} (y) = F^-1 {i sgn(k) (F{f} (k))} (y)

where F and F^-1 are the Fourier and inverse Fourier transforms. The middle step doesn't appear to quite be a phase shift--it's more futzing with the negative frequency components somehow.

[cristianoych]

Hello,

Hello jemidiah, the FFT => Phase Shift => IFFT method should be correct, because I study from a reference book. By the way, if you would like to know more about the Hilbert transform, you can take a look on this website: http://www.complextoreal.com/tcomplex.htm

To zaza, may i know what is Numerical Recipes? Sorry i not so good in C.

By the way, thanks to jemidiah and zaza for helps.

Below are some codes that i found:

First

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N 1000


struct complex {float real,imag;};             
void fft(complex*x,int n);
void ifft(complex *x,int n);
void hilbert(float *x,int n,float *y)      
{
     struct compx *z;
     int i,k;
     float temp;
     z=(struct compx *)malloc(sizeof(struct compx)*n);
     for(i=0;i<n;i++)
     {
      z[i].real=x[i];
      z[i].imag=0;
      }
     fft(z,n);
     k=n/2;
     z[0].real=0;
     z[0].imag=0;
     z[k].real=0;
     z[k].imag=0;
     for(i=1;i<k;i++)                      
     {
      temp=z[i].real;
      z[i].real=-z[i].imag;
      z[i].imag=temp;
      } 
     for(i=k+1;i<n;i++)                    
     {
      temp=-z[i].real;
      z[i].real=z[i].imag;
      z[i].imag=temp;
      }
     ifft(z,n);
     for(i=0;i<n;i++)
     x[i]=z[i].real;
     free(z);
     }
      
void fft(){
int i=0,j=0,k=0,l=0;
complex up,down,product;
change();
for(i=0;i< log(size_x)/log(2) ;i++){   
   l=1<<i;
   for(j=0;j<size_x;j+= 2*l ){             
    for(k=0;k<l;k++){       
      mul(x[j+k+l],W[size_x*k/2/l],&product);
      add(x[j+k],product,&up);
      sub(x[j+k],product,&down);
      x[j+k]=up;
      x[j+k+l]=down;
    }
   }
}
}

void ifft(){
int i=0,j=0,k=0,l=size_x;
complex up,down;
for(i=0;i< (int)( log(size_x)/log(2) );i++){   
   l/=2;
   for(j=0;j<size_x;j+= 2*l ){            
    for(k=0;k<l;k++){        
     add(x[j+k],x[j+k+l],&up);
     up.real/=2;up.img/=2;
     sub(x[j+k],x[j+k+l],&down);
     down.real/=2;down.img/=2;
     divi(down,W[size_x*k/2/l],&down);
     x[j+k]=up;
     x[j+k+l]=down;
    }
   }
}
change();
}

Second

/* Perform Hilbert transform and calculate instantaneous magnitude and frequency.
*/

#include <stdio.h>
#include <math.h>

#define LMAX 200000
#define LFILT 128	/* impulse response filter length; must be even!!! */


main(argc, argv)
int argc;
char *argv[];
{

    int i, npt, lfilt;
    static double time[LMAX+1], x[LMAX+1];
    static double xh[LMAX+1], ampl[LMAX+1], phase[LMAX+1], omega[LMAX+1];
    double hilb[LFILT+1], pi, pi2, xt, xht;

    pi = M_PI;
    pi2 = 2*pi;
    lfilt = LFILT;

    if (argc != 1) {
        usage(argv[0]);
        exit(1);
    }
 
    for (i=1; i<=lfilt; i++)
        hilb[i]=1/((i-lfilt/2)-0.5)/pi;

    for (i=1; i<=LMAX && scanf("%lf %lf", &time[i], &x[i])==2; i++) {
        xh[i] = 0.0;
        ampl[i] = 0.0;
    }
    npt=i-1;

    if (npt < lfilt) {
        fprintf(stderr, "%s : insufficient data\n", argv[0]);
	exit(2);
    }

    /* Hilbert transform */
      convol(x, xh, hilb, npt, lfilt);

    /* Ampl and phase */
    for (i=lfilt/2+1; i<=npt-lfilt/2; i++) {
        xt = x[i];
        xht = xh[i];
        ampl[i] = sqrt(xt*xt+xht*xht);
        phase[i] = atan2(xht ,xt);
        if (phase[i] < phase[i-1])
            omega[i] = phase[i]-phase[i-1]+pi2;
        else
            omega[i] = phase[i]-phase[i-1];
    }

    for (i=lfilt/2+2; i<=npt-lfilt/2; i++) 
        printf("%g %.6f %.6f\n", time[i], ampl[i], omega[i]/pi2);

}


convol(source, target, filt, npt, lfilt)
int npt, lfilt;
double *source, *target, *filt;
{
    int i, l;
    double yt;

    for (l=1; l<=npt-lfilt+1; l++) {
        yt = 0.0;
        for (i=1; i<=lfilt; i++) 
            yt = yt+source[l+i-1]*filt[lfilt+1-i];
        target[l] = yt;
    }

    /* shifting lfilt/1+1/2 points */
    for (i=1; i<=npt-lfilt; i++) {
        target[i] = 0.5*(target[i]+target[i+1]);
    }
    for (i=npt-lfilt; i>=1; i--)
        target[i+lfilt/2]=target[i];

    /* writing zeros */
    for (i=1; i<=lfilt/2; i++) {
        target[i] = 0.0;
        target[npt+1-i] = 0.0;
    }
}


usage(prog)
char * prog;
{
    fprintf(stderr, "Usage : %s\n\n", prog);
    fprintf(stderr, " Reads 2 columns of input (time and x) and outputs\n");
    fprintf(stderr, " time and Hilbert transform ampltudes and frequencies.\n");
}

Are those codes perform the hilbert function? Or anyone can provides a new code ? Thanks.

[jemidiah]

The first code uses the formula I quoted in post #4. The second might be using the convolution definition labeled (2) in your link. I'm not qualified to say if they're good implementations or how well they work in edge cases.

Fundamentally, the formula I quoted isn't very hard to apply. You get FFT and IFFT routines, apply the FFT to your signal, multiply the coefficients generated by sgn(k)*i, and apply the IFFT. FFT and IFFT routines should be plentiful, and your first code includes what I hope is a good implementation of each.


It's unclear to me what precisely you meant by "FFT => Phase Shift => IFFT method". If you mean you're phase shifting the Fourier coefficients, then I agree (you're just multiplying by a unit complex number, which is by definition a phase shift). At first I took you to mean phase shifting the time-domain signal, or phase shifting the entire frequency domain signal uniformly, either of which just didn't make sense.


Also, Numerical Recipes is a popular book which gives implementations of numerous numerical problems. Chapter 12 seems to be devoted to the FFT.