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
Instantaneous frequency with Hilbert transforms!!
Started by ●July 17, 2008
Reply by ●July 17, 20082008-07-17
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 youI 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
Reply by ●July 17, 20082008-07-17
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 youI'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
Reply by ●July 17, 20082008-07-17
On Jul 17, 8:53�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!! > > � �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, > Daviddeltou, 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
Reply by ●July 18, 20082008-07-18
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, > > ClayAnother 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
Reply by ●July 18, 20082008-07-18
On Jul 18, 10:06�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!! > > > � �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- 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
