instantaneous frequency

Dear Sir/Madam,
I have some problems with calculating the instantaneous frequency of a
signal which may be so frequent and familiar, but I couldn't solve them.
The signal is a special ground motion acceleration, El Centro. For the
purpose of controlling a structure the frequency in real time is needed.
Using Hilbert transformation results in two problems. The first one is over
shooting and the second one is that negative values are calculated as
instanteneous frequency for some time steps!! I've used this method:

record=load('el_centro.txt');
t=record(:,1);
sig=record(:,2);
H=hilbert(sig);
phi=unwrap(phase(H));
freq=diff(phi)./diff(t);

I would be very grateful if you tell me why the negative frequencies are
calculated and also a better way to calculate the instantaneous frequency
without over shooting. 
Using spectrogram is just a suitable way to plot a surface, however, I
don't know how to use it to calculate what I want.

Sincerely,

On May 16, 7:40&#4294967295;am, "rmgh" <rmgh_65@n_o_s_p_a_m.yahoo.com> wrote:
> Dear Sir/Madam,
> I have some problems with calculating the instantaneous frequency of a
> signal which may be so frequent and familiar, but I couldn't solve them.
> The signal is a special ground motion acceleration, El Centro. For the
> purpose of controlling a structure the frequency in real time is needed.
> Using Hilbert transformation results in two problems. The first one is over
> shooting and the second one is that negative values are calculated as
> instanteneous frequency for some time steps!! I've used this method:
>
> record=load('el_centro.txt');
> t=record(:,1);
> sig=record(:,2);
> H=hilbert(sig);
> phi=unwrap(phase(H));
> freq=diff(phi)./diff(t);
>
> I would be very grateful if you tell me why the negative frequencies are
> calculated and also a better way to calculate the instantaneous frequency
> without over shooting.
> Using spectrogram is just a suitable way to plot a surface, however, I
> don't know how to use it to calculate what I want.

Before you can measure it, you need to define it. Knowing that
frequency is the time rate of change of phase is small help. How do
you measure instantaneous phase? Phase relative to what? Nail these
down, and you fill likely have your answer.

Jerry
--
Engineering is the art of making what you want from things you can get.

On May 16, 7:40&#4294967295;am, "rmgh" <rmgh_65@n_o_s_p_a_m.yahoo.com> wrote:
> Dear Sir/Madam,
> I have some problems with calculating the instantaneous frequency of a
> signal which may be so frequent and familiar, but I couldn't solve them.
> The signal is a special ground motion acceleration, El Centro. For the
> purpose of controlling a structure the frequency in real time is needed.
> Using Hilbert transformation results in two problems. The first one is over
> shooting and the second one is that negative values are calculated as
> instanteneous frequency for some time steps!! I've used this method:
>
> record=load('el_centro.txt');
> t=record(:,1);
> sig=record(:,2);
> H=hilbert(sig);
> phi=unwrap(phase(H));
> freq=diff(phi)./diff(t);
>
> I would be very grateful if you tell me why the negative frequencies are
> calculated and also a better way to calculate the instantaneous frequency
> without over shooting.
> Using spectrogram is just a suitable way to plot a surface, however, I
> don't know how to use it to calculate what I want.
>
> Sincerely,

The method you have used is extremely sensitive to noise. You could
instead try:
freq = angle(H(1:end-1).*conj(H(2:end)));

You'll need to divide by the appropriate scale factor. The above may
be slightly better but is also sensitive to noise.

A positive frequency is one where the instantaneous angle increases
with time. A negative frequency indicates the angle is decreasing with
time. Note - the angle I talking about is the parameter in the complex
exponential i.e. exp(j* theta).

The instantaneous frequency is not well defined in the case when
multiple sinusoids are present, or even in the case of a single
sinusoid plus noise. I would suggest you search this forum. There have
been multiple discussions in the past about instantaneous frequency.

Cheers,
David

On May 17, 8:36&#4294967295;am, Dave <dspg...@netscape.net> wrote:
> On May 16, 7:40&#4294967295;am, "rmgh" <rmgh_65@n_o_s_p_a_m.yahoo.com> wrote:
>
>
>
> > Dear Sir/Madam,
> > I have some problems with calculating the instantaneous frequency of a
> > signal which may be so frequent and familiar, but I couldn't solve them.
> > The signal is a special ground motion acceleration, El Centro. For the
> > purpose of controlling a structure the frequency in real time is needed.
> > Using Hilbert transformation results in two problems. The first one is over
> > shooting and the second one is that negative values are calculated as
> > instanteneous frequency for some time steps!! I've used this method:
>
> > record=load('el_centro.txt');
> > t=record(:,1);
> > sig=record(:,2);
> > H=hilbert(sig);
> > phi=unwrap(phase(H));
> > freq=diff(phi)./diff(t);
>
> > I would be very grateful if you tell me why the negative frequencies are
> > calculated and also a better way to calculate the instantaneous frequency
> > without over shooting.
>
...

> The method you have used is extremely sensitive to noise. You could
> instead try:
> freq = angle(H(1:end-1).*conj(H(2:end)));
>
> You'll need to divide by the appropriate scale factor.

it's the same as the diff(t) that the OP is using which i presume is
constant (if uniformly sampled).

> The above may be slightly better but is also sensitive to noise.

i think it should come out the same except, maybe, at places where
unwrap() doesn't unwrap it the same way.  but if the abs value of the
phase change (between adjacent samples of H()) after unwrap is less
than pi, i would think that the two methods would come out the same,
noisy or not.

but i heartily agree with Dave that, to calculate instantaneous
frequency (say for FM demodulation), calculating the angle between
adjacent samples of the analytic signal (which is really what MATLAB's
hilbert() calculates) is better than unwrapping and calculating the
difference between adjacent angles.

you get a nice consistent formula (without any "if" or branch
instructions) that comes up with the discrete difference of phase
(which can be filtered to get something closer to the actual
derivative of phase, with a little delay).  this filter can be teamed
up with whatever de-emphasis (undoing the pre-emphasis done at the FM
transmitter).

r b-j