Calculate and interpret the instantaneous frequency

Question

I'm new to the principle of calculating the instantaneous frequency, and came up with a lot of questions on it. You find them all in a bullet-point list at the end of this text. The text might be a little long, excuse me for that, but I really tried to work on that problem on my own.

So I'm interested in the instantaneous frequency $f(t)$ of a real valued signal $x(t)$. The calculation is done with the help of an analytic signal $z(t) = x(t) + j y(t)$, where $y(t)$ is the Hilbert transformation of $x(t)$.

To calculate the instantaneous frequencies from the analytic signal $z(t)$ I followed the paper:

The calculation of instantaneous frequency and instantaneous bandwidth by Arthur E. Barns from 1992. In this paper he introduces multiple methods to calculate the instantaneous frequency. I write down, all formulas he proposed (and I used) in a moment.

For "learning", I played around with a very simple, and two a little more complex signals, in MATLAB, and wanted to get their instantaneous frequencies.

Fs = 1000;                                            % sampling-rate = 1kHz
t = 0:1/Fs:10-1/Fs;                                    % 10s 'Timevector'
chirp_signal = chirp(t,0,1,2);                         % 10s long chirp-signal, signal 1
added_sinusoid = chirp_signal + sin(2*pi*t*10);        % chirp + sin(10Hz), signal 2
modulated_sinusoid = chirp_signal .* sin(2*pi*t*10);   % chirp * sin(10Hz), signal 3

The plots in time domain of those three signals look the following:

The plots of all instantaneous frequencies I got after applying all methods from the paper, are the following:

Instantaneous frequencies of pure chirp signal: Instantaneous frequencies of chirp signal with added sinusoid: Instantaneous frequencies of modulated chirp signal: Please note, that in all three images, the y-axis of plot 3 and 4 are zoomed in, so the amplitudes of those signals are very tiny!

The first possibility to get from the analytic signal to the instantaneous frequency is: \begin{equation} f_{2}(t) = \frac{1}{2\pi} \frac{d}{dt}\theta(t) \end{equation} where $\theta(t)$ is the instantaneous phase. I think this is the most common method used today, at least on MATLAB's webpage it is calculated that way. The code looks the following:

function [instantaneous_frequency] = f2(analytic_signal,Fs)
    factor =  Fs/(2*pi);
    instantaneous_frequency = factor * diff(unwrap(angle(analytic_signal)));
    % Insert leading 0 in return-vector to maintain size
    instantaneous_frequency = [0 instantaneous_frequency];
end

In the paper Barns now suggest (or rather said compiles) four other ways to calculate instantaneous frequencies from the analytic signal. He also mentions the upper formula, but is the opinion that it is impractical due to ambiguities in the phase. I guess, he did not know of the unwrap() method, or to be more precise the math behind it. (I, myself did learn about that method just today, when looking at some other source codes on instantaneous frequencies)

In his paper, the formula has label Number (2), therefore, I gave the f(t) the index 2. All other indexes correspond the same way to their numbers in the paper.

Because of the ambiguities in phase, he rather suggests:

\begin{equation} f_{3}(t) = \frac{1}{2\pi}\cdot\frac{\overbrace{x(t)y'(t)}^\text{a}-\overbrace{x'(t)y(t)}^{\text{b}}}{\underbrace{x(t)^2}_{\text{c}}+\underbrace{y(t)^2}_{\text{d}}} \end{equation} I introduced the symbols "a","b","c" and "d", to make the programming a little easier:

function [instantaneous_frequency] = f3(analytic_signal,Fs,T)
    x = real(analytic_signal);
    y = imag(analytic_signal);
    diff_x = diff(x);
    diff_y = diff(y);
    factor = Fs/(2*pi);
    a = x(2:end).*diff_y;
    b = y(2:end).*diff_x;
    c = x(2:end).^2;
    d = y(2:end).^2;
    instantaneous_frequency = factor * ((a-b)./(c+d));
    % Insert leading 0 in return-vector to maintain size
    instantaneous_frequency = [0 instantaneous_frequency];
end

Then Barner gives three more formulas which he names "instantaneous frequency approximations": \begin{equation} f_{9}(t) = \frac{1}{2\pi T}\cdot\text{arctan}\left[\frac{\overbrace{x(t)y(t+T)}^{\text{a}}-\overbrace{x(t+T)y(t)}^{\text{b}}}{\underbrace{x(t)x(t+T)}_{\text{c}}+\underbrace{y(t)y(t+T)}_{\text{d}}}\right] \end{equation}

function[instantaneous_frequency] = f9(analytic_signal, Fs, T)
    x = real(analytic_signal);
    y = imag(analytic_signal);
    factor = Fs/(2*pi*T);
    a = x(1:end-T).*y(1+T:end);
    b = x(1+T:end).*y(1:end-T);
    c = x(1:end-T).*x(1+T:end);
    d = y(1:end-T).*y(1+T:end);
    instantaneous_frequency = factor.*atan((a-b)./(c+d));
    % Append 0 to return-vector to maintain size
    instantaneous_frequency = [instantaneous_frequency zeros(1,T)];
end

\begin{equation} f_{11}(t) = \frac{1}{4\pi T}\cdot\text{arctan}\left[\frac{\overbrace{x(t-T)y(t+T)}^{\text{a}}-\overbrace{x(t+T)y(t-T)}^{\text{b}}}{\underbrace{x(t-T)x(t+T)}_{\text{c}}+\underbrace{y(t-T)y(t+T)}_{\text{d}}}\right] \end{equation}

function [instantaneous_frequency] = f11(analytic_signal, Fs, T)
    x = real(analytic_signal);
    y = imag(analytic_signal);
    factor = Fs/(4*pi*T);
    a = x(1:end-2*T).*y(1+2*T:end);
    b = x(1+2*T:end).*y(1:end-2*T);
    c = x(1:end-2*T).*x(1+2*T:end);
    d = y(1:end-2*T).*y(1+2*T:end);
    instantaneous_frequency = factor.*atan((a-b)./(c+d));
    % Append and insert 0s to maintain size
    instantaneous_frequency = [zeros(1,T) instantaneous_frequency zeros(1,T)];
end

\begin{equation} f_{14}(t)=\frac{2}{\pi T}\left[\frac{\overbrace{x(t)y(t+T)}^{\text{a}}-\overbrace{x(t+T)y(t)}^{\text{b}}}{\underbrace{(x(t)+x(t+T))^2}_{\text{c}}+\underbrace{(y(t)+y(t+T))^2}_{\text{d}}}\right] \end{equation}

function [instantaneous_frequency] = formula14(analytic_signal, Fs, T);
    x = real(analytic_signal);
    y = imag(analytic_signal);
    factor = 2*Fs/(pi*T);
    a = x(1:end-T).*y(1+T:end);
    b = x(1+T:end).*y(1:end-T);
    c = (x(1:end-T)+x(1+T:end)).^2;
    d = (y(1:end-T)+y(1+T:end)).^2;
    instantaneous_frequency = factor * ((a-b)./(c+d));
    % Append and insert 0s to maintain size
    instantaneous_frequency = [instantaneous_frequency zeros(1,T)];
end

In all 3 approximation Formulas T was set to Fs (T = Fs = 1000 = 1s), as suggested in the paper.

Now my question's are:

• Formulas f2 and f3 return the same result for the pure chirp signal. I think that's good, as they do calculate the same. The three approximation methods do not return the same, not even something which is close to it! Why is that the case? (I hope it is not just an programming-bug...)
• Although they return the same, especially at the end of the plot they start to 'wiggle' a lot. What is the explanation for that? I first thought of something like aliasing, but my sampling frequency is quite high, compared to the signals' frequency, so I think that can be excluded.

• At least f2 and f3 seem to work appropriate on a pure chirp signal, but all methods, including f2 and f3 seem to fail horrible, when it comes to more than one frequency in the signal. In reality having more than one frequency in a signal is rather always the case. So how can one get the (more or less) correct instantaneous frequency?

  • I actually don't even know what to expect, when more than one frequency is present in the signal. The calculation returns one number for a given point in time, so what should it do when, like here, more frequencies are present? Return the average of all frequencies or something like that?

• And my probably most important question is, how is that handled in a real and elaborated software? Let's say I want to know the instantaneous frequency of the modulated signal at 1.75 s, and I chose method f2, than I can be 'lucky' and get a number close to 6[Hz] which is most likely the correct answer, or I pick my results few samples next to it and suddenly I get some wired, way to high, result, because I unfortunately picked a value in the spike. How can this be handled? By postprocessing it with mean or even better a median filter? I think even that might get really difficult especially in regions where many spikes are next to each other.

And one last, not so important question, why is it that most papers I find on instantaneous frequencies are from the area of geography, especially in calculating seismographic events like earthquakes. Barne's paper also takes that as an example. Isn't the instantaneous frequency interesting in many areas?

That's it so far, I'm very thankful for every reply, especially when somebody gives me tips on how to implement it in an real software project ;)

Kind regards, Patrick

Answer 1

Not really an answer but maybe helpful: Personally I found that the concept of instantaneous frequency is only useful for sufficiently narrow band signals.

Consider the simple example of two steady sine waves, say 100Hz and 934Hz. In this case you can certainly define and calculate the instantaneous frequency (in whatever way you want) but what should the result be? What possible insight or properties can the instantaneous frequency have that say anything meaningful about the signal? Applying the concept of instantaneous frequency to signals that have multiple frequencies at the same time, simply doesn't make a lot of sense.

That's why you get decent results for the sweeps but odd curves for the Sweep+sine. It's also the reason why you see the wiggles a the high part of the sweep. The bandwidth of the signal gets too high to assign a single frequency number to it and so the results jumps around.

Answer 2

At least f2 and f3 seem to work appropriate on a pure chirp signal, but all methods, including f2 and f3 seem to fail horrible, when it comes to more than one frequency in the signal. In reality having more than one frequency in a signal is rather always the case. So how can one get the (more or less) correct instantaneous frequency?

as Hilmar suggests, the Hilbert transform (or "Analytic Signal") method does not work on wide-band because there are more than one frequency component. you can do this method only for a single sinusoidal component.

so, with the Analytic Signal approach, what you want to do is make use of this identity:

$$ \arctan u - \arctan v = \arctan \left( \frac{u-v}{1+uv} \right) $$

if $|u-v|$ is small enough, which you can derive Barner's "$f_9$" formula from.

but there must be only one time-varying sinusoid in the Hilbert transform calculation to do that correctly. and you better line up the "in-phase" component with the output of the Hilbert transform (which is delayed with a causal FIR filter). otherwise you'll get crap.

Answer 3

Sorry to provide an answer a year after the fact, but I stumbled across this post while searching for articles on this very topic. Your questions reflect the widespread disagreements and interpretations of "instantaneous frequency" that have plagued the field since its inception. Numerous people will tell you, as some of the answers here, that IF is only applicable to "narrowband" or "mono-component" signals. In fact, that is not true: sometimes the IF obtained by the Hilbert transform is perfectly well-behaved for broadband and/or "multi-component" signals. One quantity that has been proposed that avoids many of these difficulties is the "weighted average instantaneous frequency (WAIF)", which can be measured using a spectrogram.

See Loughlin in J. Acoust. Soc. Am., Jan. 1999. Other good papers on IF and common misconceptions are by Picinbono (IEEE Trans. Sig. Proc., March 1997) and Vakman (IEEE Trans. Sig. Proc., April 1996).

Answer 4

Wow, what a huge question. I'm going to answer the not-so-important question first:

And one last, not so important question, why is it that most papers I find on instantaneous frequencies are from the area of geography, especially in calculating seismographic events like earthquakes. Barne's paper also takes that as an example. Isn't the instantaneous frequency interesting in many areas?

The reason is that the seismographic system "vibroseis" is used in the oil industry to do seismic surveys. The trucks I've linked to vibrate from about 5 Hz to about 90 Hz and can be made to do chirp signals. There is much money in the oil industry, and processing the returns from these signals can be very, very lucrative. Hence, many people have spent many hours analyzing such signals, including looking at instantaneous frequency techniques.

---

As for your more important questions: generally, doing arithmetic differences and calculating arctangents on discrete-time signals is a Bad Thing$^{\rm TM}$. This is because discrete-time frequency estimates need to be calculated using "circular arithmetic" (AKA vector arithmetic).

Check out this paper.

Better approaches tend to use "phase weighted averagers" as implemented here. Or here for a direct link to the matlab.