Numerical implementation of Hilbert transform

I have been trying to implement the Hilbert transform numerically but I have been having some trouble. The way I decided to do it is the use of Fourier transforms. The fourier transform of the Hilbert transform of f(x) is -i*sgn(k)*F(k), where F(k) is the Fourier transform of f(x).

I have written my own Fourier and inverse Fourier routines but I don't get nice results. Can anyone suggest anything? The picture (pdf) I have included shows my computational result and the value which it should be.

[John BG]

Have you tried the Signal Processing Toolbox function hilbert.m?

MuPAD also has the symbolic hilbert transform function.

Or perhaps you want to use this one:

http://uk.mathworks.com/matlabcentral/fileexchange/29072-hilbert2?s_tid=srchtitle

This is MATLAB, not C C# C++ or else, you don't have to write your functions when there are already working and robust FFT and hilbert functions available.

[Star Strider]

It would definitely help if you describe what you are doing in your functions, what the arguments are, and whether they are scalars or vectors. Without that, it’s impossible for me to understand them. The easiest way to comment them is to put several tabs to the right of each line and then describe what the line is doing. For example:

function y=ft(x,f,k)
n=length(k);                                            % Frequency Vector ‘k’
y=zeros(1,n);                                           % Output Vector ‘y’
for i=1:n
    v_1=exp(-sqrt(-1)*k(i)*x);                          % Calculate Complex Frequency Vector ‘v1’
    v_2=f.*v_1;                                         % Multiply By Time-Domain Signal ‘f’
    y(i)=-real(trapz(v_2,x));                           % Integrate To Get Fourier Component At ...
end

I don’t know if any of these comments are actually correct because I don’t know what your function is doing or what the variables are. That’s the problem.

[Mat Hunt]

John BG, It would be nice if there was a dedicated fourier transform, but there isn't. I really have no idea what fft and ifft does. They're certainly not what I would call a Fourier transform. For a Fourier transform, you require f(x),x and a value k (called the wavenumber). fft simply takes a function f(x) and as a result I have no idea what the hell it's doing. So as a general rule, I write off such things and write things myself, so I KNOW what I am getting.

I want a numerical way of computing the Hilbert transform. I don't have the signal processing suite, so I can't use that function. The other function you suggested also requires the Hilbert function.

Star Strider, The programs are written for clarity. To understand then you actually have to know what a Fourier transform is mathematically. Once you know this, then the programs are pretty straightforward to understand. However I will comment them here:

function y=ft(x,f,k) %computes a FT for wavevector k
n=length(k); %Checks the length of the wavevector
y=zeros(1,n); %The output is the same length
for i=1:n
    v_1=exp(-sqrt(-1)*k(i)*x); %computes exp(-ikx)
    v_2=f.*v_1; %computes the integrand of FT
    y(i)=-trapz(v_2,x); %Integrates integrand
end

function y=ift(k,f_hat,x) %Inverse FT
n=length(x); %Computes length of vector x
y=zeros(1,n); %Answer has same length 
a=1/(2*pi); %Scaling factor of 1/2π
for i=1:n
v_1=exp(sqrt(-1)*k*x(i)); %Computes exp(ikx)
v_2=f_hat.*v_1; %Computes integrand f̂(k)exp(ikx)
y(i)=-a*trapz(v_2,k); %Integrates
end

function y=hilbert(x,f_x,z)

%This takes numerical input (x,f(x)) and evaluates the Hilbert transform at
%x=z;
k=-4:0.001:4; %Selects range of Fourier variable k
f_x_hat=ft(x,f_x,k); %Computes Fourier transform of f(x)
dum=-sqrt(-1)*sign(k).*f_x_hat; %Computes -isgn(k)f̂(k)
y=ift(k,dum,z); %Computes inverse Fourier transform to get the Hilbert transform

[Star Strider]

I understand about Fourier transforms, particularly for signal processing. It is not obvious what you are using them for, and what your variables are.

You can also substitute 1i for ‘sqrt(-1)’ eliminating the need for the evaluation and making for more efficient code.

You obviously know infinitely more about MATLAB coding of the Fourier and related integral transforms than we do, so I will leave you to it.

[Mat Hunt]

I am using them to compute the Hilbert transform. It's a well known result that the Fourier transform of the Hilbert transform is -i*sgn(k)*Fourier transform of the original function.

If you have a look at the pdf files, you will see my result, and what the answer should be. I am confused about how I get there as there is a lot of oscillation in my answer where there shouldn't be.

[Mat Hunt]

I've tried the li thing, it just throws up errors.

Can anyone suggest a solution to the problem.

[Walter Roberson]

"For a Fourier transform, you require f(x),x and a value k (called the wavenumber)"

Have you told Jean-Baptiste that?

[Mat Hunt]

Any actual answers?

[Walter Roberson]

Mat Hunt comments to me

I assume that you're unaware how much the Fourier transform is used in the area of surface wave modelling where the k variable is called different things.

[Walter Roberson]

Mat, you are free to call anything you want a "Fourier transform", and to refuse to acknowledge anything other people call a "Fourier transform" as being what you would call a Fourier transform. But if you do so, then we have no common ground and you are going to need to solve the problem yourself.

The transform invented by Jean-Baptiste Joseph Fourier does not require a "k". There are alternative formulations that involve "k" and those formulations make it easier to approach various computations, but they are not the only possibilities.

fft() and ifft() do not implement the Fourier Transform: they implement the Fast Fourier Transform, and the Inverse Fast Fourier Transform, which are implementations of Discrete Fourier Transforms. The Symbolic Toolbox implements the (continuous) Fourier Transform http://www.mathworks.com/help/symbolic/fourier.html, and does so without explicit reference to the variables you indicate as being required. When I look at your code, I am pretty much certain that you are not interested in the continuous Fourier Transform and are instead interested in a Discrete Fourier Transform. However, I have only been reading code for 40 years so I could be mistaken.

[Mat Hunt]

I would like to implement the Fourier transform as defined as follows:

F(f)(k)=\int_{-\infty}^{\infty}f(x)exp(-ikx)dx

I only have a limited number of points to work with, so I have to do an approximation to the above transform. This will work relatively well if f(x) tends towards zero reasonably rapidly. This is the case I have, So I chose a few points and computed the Fourier transform and the results are actually quite good. The only snag is that you have to know what the limits of your function are.