Discrete Hilbert transformation of a spatial signal

Question

first of all, sorry for this basic question in signal processing, I'm pretty new in this field and have been stuck in a situation a little beyond my expertise.

The idea is to derive the instantaneous wavenumber direction as an estimation of the direction normal of a 2D field by using a directional Hilbert transform in one of the directions.

As far as I can recall correctly, all I have to do is to do a 1D Fourier transform of each data vector, apply a $\pi/2$ phase-shift in the wavenumber domain and then apply an inverse Fourier transformation to get the Hilbert transform.

So the Hilbert transform is defined as: $$ H(x(t)) = \frac{1}{\pi} PV \int_{-\infty}^{\infty}\frac{x(\tau)}{\tau-t}d\tau $$ where $PV$ denotes the principal value of the integral and $H$ is the Hilbert transform of a vector $x(t)$. This is supposed to be equal to the convolution $\frac{1}{\pi t} * x(t)$ and

$$ F(H(x(t))) = F(\frac{1}{\pi t}) F(x(t)). $$

So the Fourier transformation and multiplication of the two factors on the right side is equal to the Fourier transformation of the Hilbert transformation of my vector, right?

My question now is how to apply the phase shift on the Hilbert transform? is it simply multiplying the product on the right side by $\pi/2$ in the wavenumber domain?

Thanks a lot for any help! Cheers!

UPDATE:

So I coded the Hilbert transformation in $C$ but unfortunately I do not get a phase shift of $90^\circ$ but somewhat different from that and with a wrong amplitude. My code looks like this:

    void hilbertfun(float *vin, float *vout, int n, int d)
    {
static complex *A;           //working array
int i;
float pi = 4. * atan(1.);
int signfun;
static float *rkx;

void rk(float *ak, int n, float d, int ind);

A  = (complex*) alloc1(n,sizeof(complex));
rkx = (float*) alloc1(n,sizeof(float));

//copy input array to complex array
for(i=0; i<n; i++){
    A[i].r = vin[i];
    A[i].i = 0.;
}

// perform inverse fft
fouriertransform(1,n,A);

// calculate wavenumbers
wavenumbers(rkx,n,d,1);



// check wether wavenumbers are positive or negative >> apply phase shift by using the impulse response of the Hilbert transformer
// + apply phase shift
for(i=0;i<n;i++){
  if (rkx[i] > 0) signfun = 1;
  else if (rkx[i] == 0) signfun = 0;
  else if (rkx[i] < 0) signfun = -1;
    A[i].r = - A[i].i * signfun;
    A[i].i = A[i].r * signfun;
}

// perform forward fft
fouriertransform(-1,n,A);

for(i=0;i<n;i++) vout[i]=A[i].r/n; //write into output and normalise

}

The Fourier transform and the normalisation work perfectly fine without the phase shift (Original (red), Output (green)):

Yet, with the phase shift I get this result (Original (red), Hilbert transformed Output (green)):

What am I missing?

Answer 1

You need the Fourier transform of the impulse response of the Hilbert transformer

$$h(t)=\frac{1}{\pi t}\Longleftrightarrow H(j\omega)=-j\textrm{ sign}(\omega)$$

So all positive frequencies are multiplied by $-j$ and all negative frequencies are multiplied by $j$. This corresponds to a phase shift of $-\pi/2$.

EDIT:

I'm referring to your updated question including the code. I can't see if the computation of the wavenumber is OK, but one problem I definitely see is in applying the phase shift. You change A[i].r and then you re-use it (as if it were the orignal value) to compute A[i].i. Furthermore, I think the signs are wrong. If signfun=1, A[i].r=A[i].i and similarly for the computation of A[i].i. But the most important thing is to use a temporary variably to store the old value of A[i].r for computing A[i].i.