# Instantaneous frequency with Hilbert transforms!!

Started by July 17, 2008
```Sorry for post again if posted before!!
About the  Instantaneous frequency with Hilbert transforms.
I still haven't found a clear  answer. This is the code I have so far
in matlab:

[x,fs] = wavread('xx.wav');  %a real speech signal;
xx = hilbert(x);             %analytic signal xx;
n = length(xx);
pha = angle(xx)/2/pi;
ff   = diff(pha);                       %ff is what i want ,that is to
say , ff is the Instantaneous frequency of the signal????

thank you

```
```On Jul 17, 8:45 am, "deltuo" <del...@yahoo.com.cn> wrote:
> Sorry for post again if posted before!!
>    About the  Instantaneous frequency with Hilbert transforms.
>    I still haven't found a clear  answer. This is the code I have so far
> in matlab:
>
>    [x,fs] = wavread('xx.wav');  %a real speech signal;
>    xx = hilbert(x);             %analytic signal xx;
>    n = length(xx);
>    pha = angle(xx)/2/pi;
>    ff   = diff(pha);                       %ff is what i want ,that is to
> say , ff is the Instantaneous frequency of the signal????
>
>   thank you

I believe what you've done is theoretically correct. But calculating
the instantaneous frequency this way is very sensitive to noise. As an
example you can end up with frequencies well outside of your sampling
frequency bandwidth.

Good luck.
Cheers,
David
```
```On Jul 17, 8:45&#4294967295;am, "deltuo" <del...@yahoo.com.cn> wrote:
> Sorry for post again if posted before!!
> &#4294967295; &#4294967295;About the &#4294967295;Instantaneous frequency with Hilbert transforms.
> &#4294967295; &#4294967295;I still haven't found a clear &#4294967295;answer. This is the code I have so far
> in matlab:
>
> &#4294967295; &#4294967295;[x,fs] = wavread('xx.wav'); &#4294967295;%a real speech signal;
> &#4294967295; &#4294967295;xx = hilbert(x); &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; %analytic signal xx;
> &#4294967295; &#4294967295;n = length(xx);
> &#4294967295; &#4294967295;pha = angle(xx)/2/pi;
> &#4294967295; &#4294967295;ff &#4294967295; = diff(pha); &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; %ff is what i want ,that is to
> say , ff is the Instantaneous frequency of the signal????
>
> &#4294967295; thank you

I'm not a matlab guy, but I do see some potential errors here with
your formulation.

First your analytic signal needs to be comprised of both the original
signal and its Hilbert transform. Thus AS(t) = f(t) + j*H(f(t))

where f(t) is the original real valued signal and H(f(t)) is the
Hilbert xform of the original signal.

Also realize that a difference of arctangents is not the same as the
arctangent of a difference.

Ideally the instantaneous freq is w =   ( I(t)*Q'(t) - Q(t)*I'(t) )/
( I(t)*I(t) + Q(t)*Q(t) )

Here I(t) is the original signal and Q(t) is the 90 degree phase
shifted signal.

A standard thing to remember is to high the filter delays matched.
I.e., don't use your original signal mixed with a delayed Hilbert
version.

If you are highly oversampled then your approximation can work.

But as another poster stated this can be rather sensitive to noise.
Also the Hilbert transform is not really the best approach since you
don't need absolute pohase info with respect to the start of the data,
you only need relative phase info.

What is better is to generate a pair of phase orthogonal filters (for
example at +45 degrees and -45 degrees) and have the passband of the
filters just cover the region of interest and wipe out the rest of the
band. In fact you may even create 4 filters that represent the two
nomimal phases and their two derivatives. So with just 4 FIR filters
you can readily have your instantaneous freq (well with a delay from
the filters). Plus if the filters operate (non-zero) over a narrow
range, their overall length may be rather short.

IHTH,

Clay

```
```On Jul 17, 8:53&#4294967295;am, Dave <dspg...@netscape.net> wrote:
> On Jul 17, 8:45 am, "deltuo" <del...@yahoo.com.cn> wrote:
>
> > Sorry for post again if posted before!!
> > &#4294967295; &#4294967295;About the &#4294967295;Instantaneous frequency with Hilbert transforms.
> > &#4294967295; &#4294967295;I still haven't found a clear &#4294967295;answer. This is the code I have so far
> > in matlab:
>
> > &#4294967295; &#4294967295;[x,fs] = wavread('xx.wav'); &#4294967295;%a real speech signal;
> > &#4294967295; &#4294967295;xx = hilbert(x); &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; %analytic signal xx;
> > &#4294967295; &#4294967295;n = length(xx);
> > &#4294967295; &#4294967295;pha = angle(xx)/2/pi;
> > &#4294967295; &#4294967295;ff &#4294967295; = diff(pha); &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; %ff is what i want ,that is to
> > say , ff is the Instantaneous frequency of the signal????
>
> > &#4294967295; thank you
>
> I believe what you've done is theoretically correct. But calculating
> the instantaneous frequency this way is very sensitive to noise. As an
> example you can end up with frequencies well outside of your sampling
> frequency bandwidth.
>
> Good luck.
> Cheers,
> David

deltou,

Actually, what you are doing is nowhere near correct. I do not have
the definition of the matlab hilbert function in front of me, but I
will assume it computes the hilbert transform (normally I design a
hilbert filter with remez). First the angle is not taken of the
hilbert transform output.  The hilbert transform gets you the
imaginary part from the real input (check and see if the sign of the
output is the way you want it, if not adjust it). You will be
calculating the phase of a complex number formed from the real and
imaginary parts, which must be time aligned before the angle
computation. If hilbert() puts in a delay, you need to compensate for
the delay in the real part when you figure out the phase. Not having
them lined up properly will mess things up. Check the bandwidth of
what the hilbert filter transforms.  If it doesn't match the input you
won't be getting what you expect; first option is to similarly
bandlimit a copy of the input signal to use as the real part. What you
are doing with '2/pi' is not right, you can figure this out. You have
not taken into account possible phase wrapping that will show up in
the diff() output. You can do a simple combined diff and unwrap by
calculating the phase of the ratio of successive complex samples. The
order in the ratio must be right.

Finally, be aware that there can be interaction between the calculated
instantaneous amplitudes and frequencies of signals. I'll suspect you
will find it here.

What are you expecting when you apply this to a speach signal? When
there is no signal, but there is noise, what does the instantaneous
frequency tell you? When there is unvoiced speech, what does the
instantaneous frequency tell you? When there is voiced speech, what
does the instantaneous frequency tell you?

What are you planning on doing with result?

Dirk
```
```On Jul 17, 10:46 am, c...@claysturner.com wrote:
> On Jul 17, 8:45 am, "deltuo" <del...@yahoo.com.cn> wrote:
>
> > Sorry for post again if posted before!!
> >    About the  Instantaneous frequency with Hilbert transforms.
> >    I still haven't found a clear  answer. This is the code I have so far
> > in matlab:
>
> >    [x,fs] = wavread('xx.wav');  %a real speech signal;
> >    xx = hilbert(x);             %analytic signal xx;
> >    n = length(xx);
> >    pha = angle(xx)/2/pi;
> >    ff   = diff(pha);                       %ff is what i want ,that is to
> > say , ff is the Instantaneous frequency of the signal????
>
> >   thank you
>
> I'm not a matlab guy, but I do see some potential errors here with
> your formulation.
>
> First your analytic signal needs to be comprised of both the original
> signal and its Hilbert transform. Thus AS(t) = f(t) + j*H(f(t))
>
> where f(t) is the original real valued signal and H(f(t)) is the
> Hilbert xform of the original signal.
>
> Also realize that a difference of arctangents is not the same as the
> arctangent of a difference.
>
> Ideally the instantaneous freq is w =   ( I(t)*Q'(t) - Q(t)*I'(t) )/
> ( I(t)*I(t) + Q(t)*Q(t) )
>
> Here I(t) is the original signal and Q(t) is the 90 degree phase
> shifted signal.
>
> A standard thing to remember is to high the filter delays matched.
> I.e., don't use your original signal mixed with a delayed Hilbert
> version.
>
> If you are highly oversampled then your approximation can work.
>
> But as another poster stated this can be rather sensitive to noise.
> Also the Hilbert transform is not really the best approach since you
> don't need absolute pohase info with respect to the start of the data,
> you only need relative phase info.
>
> What is better is to generate a pair of phase orthogonal filters (for
> example at +45 degrees and -45 degrees) and have the passband of the
> filters just cover the region of interest and wipe out the rest of the
> band. In fact you may even create 4 filters that represent the two
> nomimal phases and their two derivatives. So with just 4 FIR filters
> you can readily have your instantaneous freq (well with a delay from
> the filters). Plus if the filters operate (non-zero) over a narrow
> range, their overall length may be rather short.
>
> IHTH,
>
> Clay

Another of doing it which I think amounts to the same thing as what
Clay said: multiply the signal by a complex sinusoid to move the
center frequency down to 0 Hz ( that means multiply by exp(-
j2*pi*fc*t) ) and then lowpass filter the resulting complex data
stream. You can also decimate the signal if required.

Doing it this way means you don't have to take into account the group
delay of the Hilbert transform.

Cheers,
Dave
```
```On Jul 18, 10:06&#4294967295;am, Dave <dspg...@netscape.net> wrote:
> On Jul 17, 10:46 am, c...@claysturner.com wrote:
>
>
>
>
>
> > On Jul 17, 8:45 am, "deltuo" <del...@yahoo.com.cn> wrote:
>
> > > Sorry for post again if posted before!!
> > > &#4294967295; &#4294967295;About the &#4294967295;Instantaneous frequency with Hilbert transforms.
> > > &#4294967295; &#4294967295;I still haven't found a clear &#4294967295;answer. This is the code I have so far
> > > in matlab:
>
> > > &#4294967295; &#4294967295;[x,fs] = wavread('xx.wav'); &#4294967295;%a real speech signal;
> > > &#4294967295; &#4294967295;xx = hilbert(x); &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; %analytic signal xx;
> > > &#4294967295; &#4294967295;n = length(xx);
> > > &#4294967295; &#4294967295;pha = angle(xx)/2/pi;
> > > &#4294967295; &#4294967295;ff &#4294967295; = diff(pha); &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; %ff is what i want ,that is to
> > > say , ff is the Instantaneous frequency of the signal????
>
> > > &#4294967295; thank you
>
> > I'm not a matlab guy, but I do see some potential errors here with
> > your formulation.
>
> > First your analytic signal needs to be comprised of both the original
> > signal and its Hilbert transform. Thus AS(t) = f(t) + j*H(f(t))
>
> > where f(t) is the original real valued signal and H(f(t)) is the
> > Hilbert xform of the original signal.
>
> > Also realize that a difference of arctangents is not the same as the
> > arctangent of a difference.
>
> > Ideally the instantaneous freq is w = &#4294967295; ( I(t)*Q'(t) - Q(t)*I'(t) )/
> > ( I(t)*I(t) + Q(t)*Q(t) )
>
> > Here I(t) is the original signal and Q(t) is the 90 degree phase
> > shifted signal.
>
> > A standard thing to remember is to high the filter delays matched.
> > I.e., don't use your original signal mixed with a delayed Hilbert
> > version.
>
> > If you are highly oversampled then your approximation can work.
>
> > But as another poster stated this can be rather sensitive to noise.
> > Also the Hilbert transform is not really the best approach since you
> > don't need absolute pohase info with respect to the start of the data,
> > you only need relative phase info.
>
> > What is better is to generate a pair of phase orthogonal filters (for
> > example at +45 degrees and -45 degrees) and have the passband of the
> > filters just cover the region of interest and wipe out the rest of the
> > band. In fact you may even create 4 filters that represent the two
> > nomimal phases and their two derivatives. So with just 4 FIR filters
> > you can readily have your instantaneous freq (well with a delay from
> > the filters). Plus if the filters operate (non-zero) over a narrow
> > range, their overall length may be rather short.
>
> > IHTH,
>
> > Clay
>
> Another of doing it which I think amounts to the same thing as what
> Clay said: multiply the signal by a complex sinusoid to move the
> center frequency down to 0 Hz ( that means multiply by exp(-
> j2*pi*fc*t) ) and then lowpass filter the resulting complex data
> stream. You can also decimate the signal if required.
>
> Doing it this way means you don't have to take into account the group
> delay of the Hilbert transform.
>
> Cheers,
> Dave- Hide quoted text -
>
> - Show quoted text -

Of course it does require a lot more work than taking into account the
Hilbert filter delay, and the one filter delay is really taken into
account by the other filter delay.

But there are other problems previously mentioned here that this does
take care of, like having the same bandlimiting in both the real and
imaginary parts of the IQ signal.

Dirk
```