2
$\begingroup$

I have a .wav file of human speech, sampled at 48kHz, with useful signal up to 5kHz, and almost no noise. If I do Hilbert transform in MatLab, I get the Analytic signal:

$$ s_\mathrm{a}[n] = s_\mathrm{real}[n] + j \cdot s_\mathrm{hi}[n], $$

where $s_\mathrm{hi}[n]$ is a Hilbert transformation of $s_\mathrm{real}[n]$ from the .wav file. Now for each sample I can calculate the phase:

$$ \phi[n] = \operatorname{Arctan}\left( \frac{s_\mathrm{hi}[n]}{s_\mathrm{real}[n]} \right) $$

and take time derivative from it. That will be the speed of phase changing. I don't want to call it "instantaneous frequency", that's too controversial.

Is any smarter way of doing this in MatLab? Probably this:

Fs/(2*pi)*diff(unwrap(angle(hilbert(s_real))));

I did this in Adobe Audition plugin, in Delphi. Way too many spikes, every time when the envelope $\sqrt{s_\mathrm{real}^2[n] + s_\mathrm{hi}^2[n]} $ was around zero.

$\endgroup$
  • 1
    $\begingroup$ You calculate something and then you are not happy with the results. In order to help you we need to know what you actually want to do and why you did expect that the calculation to come up with something meaningful. In general, the "thing" you calculate is useful for narrow band signals, but for broad band signals like speech, it has little usable information or properties. Why do you think otherwise ? $\endgroup$ – Hilmar Jul 9 '18 at 12:47
  • $\begingroup$ there is no controversy in calling the rate of change of phase of a sinusoid (w.r.t. time), the "instantaneous frequency". that is what it is. but if the input is not a simple sinusoid, but an added collection of sinusoids, there are many instantaneous frequencies, not just one. perhaps that is what you are finding to be controversial? $\endgroup$ – robert bristow-johnson Jul 9 '18 at 17:11
  • $\begingroup$ Thank you for the comments. If the "thing" works well for SSB radio signal, LSB or USB (I'm ham radio operator), because it is narrow band (3kHz on 3.5MHz) - why it wouldn't work when you transfer radio signal to the baseband? It should, I just calculated it wrong. $\endgroup$ – K-man Jul 10 '18 at 7:51
  • $\begingroup$ Before asking the question, I did some homework: term "instantaneous frequency" has quite a history of even academic discussions. I try to stay low, just to find the right algorithm to get the phase change rate of human speech. $\endgroup$ – K-man Jul 10 '18 at 7:59
3
$\begingroup$

I will not address the question whether or not it is meaningful to compute the instantaneous frequency of a speech signal. Instead I will show you a better method for computing the instantaneous frequency from a given analytic signal. This method avoids the phase unwrapping problem by directly computing the instantaneous frequency from the real and imaginary parts of the (discrete-time) analytic signal

$$s[n]=s_R[n]+js_I[n]\tag{1}$$

where

$$ s_I[n] = \mathscr{H} \big\{ s_R[n] \big\} $$

with $\mathscr{H}\{\cdot\}$ denoting the (discrete) Hilbert transform.

First, note that the formula for computing the phase given in your question only works if the real part of the analytic signal is greater than zero. In general you have to use the atan2 function implemented in most programming languages. In Matlab you can of course use angle, as shown in your example code.

Second, the phase becomes irrelevant and undefined if the magnitude of the respective complex number is (close to) zero. So it's pointless trying to compute the phase angle of a complex number with a negligible magnitude. This is true in general, so it doesn't matter how smart your method for computing the instantaneous frequency is.

Let's define the instantaneous frequency of a discrete-time analytic signal as the forward difference

$$f[n]=\frac{\phi[n+1]-\phi[n]}{2\pi T}\tag{2}$$

where $\phi[n]$ is the phase angle of the analytic signal $s[n]$, and $T$ is the sampling period. Note that other definitions such as the backward difference or central difference are also possible and useful.

We will not use $(2)$ directly for computing $f[n]$ because it requires phase unwrapping, which is especially problematic in the presence of noise. With

$$s[n]= \Big|s[n]\Big|e^{j\phi[n]}\tag{3}$$

we get

$$s[n+1]s^*[n]=\Big|s[n+1]\Big|e^{j\phi[n+1]}\cdot \Big|s[n]\Big|e^{-j\phi[n]}=\Big|s[n+1]\Big| \Big|s[n]\Big|e^{j(\phi[n+1]-\phi[n])}\tag{4}$$

(where $^*$ denotes complex conjugation), and, consequently,

$$f[n]=\frac{1}{2\pi T}\arg\Big\{s[n+1]s^*[n]\Big\}\tag{5}$$

The argument should be computed using the atan2 function. The result of the argument computation is only defined for $s[n]\neq 0$ and $s[n+1]\neq 0$.

Note that using backward or central differences in the definition of the instantaneous frequency would have resulted in the following expressions:

$$\begin{align}f_b[n]&=\frac{1}{2\pi T}\arg\Big\{s[n]s^*[n-1]\Big\}\\f_c[n]&=\frac{1}{4\pi T}\arg\Big\{s[n+1]s^*[n-1]\Big\}=\frac12\big(f[n]+f_b[n]\big)\end{align}$$

Expressing Eq. $(5)$ in terms of the real and imaginary parts of the analytic signal gives:

$\begin{align} f[n] &=\frac{1}{2\pi T}\arg\Big\{\big(s_R[n+1]+js_I[n+1]\big)\big(s_R[n]-js_I[n]\big) \Big\}\\ \\ &=\frac{1}{2\pi T}\arg\Big\{s_R[n]s_R[n+1]+s_I[n]s_I[n+1]+j\big(s_R[n]s_I[n+1]-s_R[n+1]s_I[n]\big)\Big\}\\ \end{align}\tag{6}$

More details about instantaneous frequency (and instantaneous bandwidth) can be found in the paper "The calculation of instantaneous frequency and instantaneous bandwidth" by A.E. Barnes.

Also take a look at this related answer.

$\endgroup$
  • $\begingroup$ looks like Matt did it for me (thanks, saves me effort). $\endgroup$ – robert bristow-johnson Jul 9 '18 at 20:31
  • $\begingroup$ Matt, the last equation (7) can be expressed with an $\arctan(\cdot)$ in it as long as we know the change in angle is no more than 90 degrees. $\endgroup$ – robert bristow-johnson Jul 9 '18 at 20:43
  • $\begingroup$ the last makes use of $$ \operatorname{Arg}(x + iy) = \operatorname{atan2}(y,\, x) = \begin{cases} \arctan\left(\frac{y}{x}\right) &\text{if } x > 0, \\ \frac{\pi}{2} - \arctan\left(\frac{x}{y}\right) &\text{if } y > 0, \\ -\frac{\pi}{2} - \arctan\left(\frac{x}{y}\right) &\text{if } y < 0, \\ \arctan\left(\frac{y}{x}\right) \pm \pi &\text{if } x < 0, \\ \text{undefined} &\text{if } x = 0 \text{ and } y = 0 \end{cases} $$ $\endgroup$ – robert bristow-johnson Jul 9 '18 at 23:08
  • $\begingroup$ Thanks Robert for correcting a mistake and making some useful changes. I just removed your final arctan formula, because I think it was wrong. In general, the real and imaginary parts can have both signs, so there simply is no single arctan formula that works in all cases. $\endgroup$ – Matt L. Jul 10 '18 at 7:55
  • 1
    $\begingroup$ @robertbristow-johnson: And the definition of $f[n]$ is really just a matter of taste, none of the two is better I think. I've used the one used in the paper I referred to. $\endgroup$ – Matt L. Jul 10 '18 at 8:12
1
$\begingroup$

instead of calculating angle (from the real and imaginary parts), unwrapping the wrapped angle, and computing a difference, the better way is to compute the difference from the real and imaginary parts directly.

this (computing instantaneous frequency from the hilbert transform) is very similar to the same phase unwrapping problem in computing group delay (from the phase response of an LTI filter). Check out here how to do that. No unwrapping is needed.

$\endgroup$
  • $\begingroup$ i'll come back here and make this answer more explicit, or find a post that had already done that. but i will do that later. $\endgroup$ – robert bristow-johnson Jul 9 '18 at 17:15
  • $\begingroup$ Hilmar asked me what I don't like in my calculations of the "thing". Spikes, as I said, which happen at zero envelope. Matt: can you elaborate more on this please: "the phase becomes irrelevant and undefined if the magnitude of the respective complex number is (close to) zero. So it's pointless trying to compute the phase angle of a complex number with a negligible magnitude." Why undefined? It probably just doesn't make any physical sense. In radio amplifiers the "noise gate" is implemented, in math I guess we should make it zero. $\endgroup$ – K-man Jul 10 '18 at 8:10
  • $\begingroup$ I tried not to deal with the phase angle directly, but used time derivative of arctan, like Barnes Equation (3). Result was the same, with spikes exceeding speech bandwidth. I will try to use this approach above. Still, I think signal around zero should be handled (cut off) manually. Not the elegant solution... $\endgroup$ – K-man Jul 10 '18 at 8:25
  • $\begingroup$ @K-man: The phase of a complex number is undefined if its magnitude is zero. Look at Eq. (3) in my answer, if $|s[n]|=0$ you can assign any value to $\phi[n]$. $\endgroup$ – Matt L. Jul 10 '18 at 16:39
0
$\begingroup$

Phase Change Rate example

I was not correct, suggesting to zeroize the rate. What I meant is on the diagram, where the original speech signal is shown, its Hilbert transform, envelope, wrapped phase, its change rate and the corresponding LSB signal.

Note obvious facts: peaks and valleys of Hilbert signal correspond to the max min change rate of original signal, and vice versa. Envelope of analytic signal is the same as the SSB envelope. Phase is decreasing, its change rate is negative, because it's Lower Side Band. Those spikes, which I was talking about, are happening at the moments when the envelope is close to zero. They even go to positive frequency area, which doesn't make sense at all. I suggest to cut them off, assigning some average red values from the adjacent samples. Same as Matt is saying that any value can be assigned at those moments.

It's easy to do in Paint, but I am looking for some math algorithm to do this automatically. That math above still requires manual exclusion of close-to-zero signals.

PS please ignore the vertical artefacts, that's calculation defect. The phase change rate is reasonably smooth.

$\endgroup$
  • $\begingroup$ Why don't you just threshold the envelope, and assign the previous valid phase rate change value to the current one if the envelope is smaller than the threshold? $\endgroup$ – Matt L. Jul 10 '18 at 20:02
  • $\begingroup$ @Matt: I will do threshold. Interestingly enough, those spikes are pretty wide, when they start - envelope is not zero yet. Initially I thought that it's not enough precision around zero, when making all this math. So I used even 768kHz sampling rate. $\endgroup$ – K-man Jul 11 '18 at 0:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.