Reply by January 24, 20152015-01-24
On Sunday, January 25, 2015 at 3:22:37 PM UTC+13, all4dsp wrote:
> >Perhaps you can tell us why, when 99.9% of all radio receivers don't use a > = > >Hilbert Transform, why you feel you need one here. I assume you are trying > = > >to recover phase. The usual method is to use a PLL on the carrier or IF > sig= > >nal - or a discriminator - quad detectors are normally the best though due > = > >to phase noise of PLLs. Also, if you want to recover phase then FFTs can > be= > > used using complex cepstrum instead of Hilbert Transforms. > > > > I measure the FM modulated carrier in the lab using bench instruments. From > that, I need to determine the modulation amplitude using any valid > technique via post-processing in Matlab. I'd don't prefer Hilbert Transform > over any other technique. However, in reading about how to perform FM > demodulation, the Hilbert Transform is often mentioned. Thus, my interest > to use it. If there are other mathematical techniques to perform FM > demodulation, they'd work too. So far I've summarized the two methods I'm > aware of above. > > _____________________________ > Posted through www.DSPRelated.com
I would have put the Hilbert Transform at the bottom of any list on FM demodulation. Quad detector being the most common. For ordinary FM demodulation is quite simple, just like AM after you defferentiate - that's a discriminator. Most of the methods work by a form of band limited differentiation. The PLL is 3dB better but as I said, making a VCO with low enough phase noise can be a problem. If I had x+jy I would use 1/(1+(x/y)^2)X differential of x/y found using a band limited differentiator of some sort as used in software radio.
Reply by all4dsp January 24, 20152015-01-24
>Perhaps you can tell us why, when 99.9% of all radio receivers don't use a
=
>Hilbert Transform, why you feel you need one here. I assume you are trying
=
>to recover phase. The usual method is to use a PLL on the carrier or IF
sig=
>nal - or a discriminator - quad detectors are normally the best though due
=
>to phase noise of PLLs. Also, if you want to recover phase then FFTs can
be=
> used using complex cepstrum instead of Hilbert Transforms. >
I measure the FM modulated carrier in the lab using bench instruments. From that, I need to determine the modulation amplitude using any valid technique via post-processing in Matlab. I'd don't prefer Hilbert Transform over any other technique. However, in reading about how to perform FM demodulation, the Hilbert Transform is often mentioned. Thus, my interest to use it. If there are other mathematical techniques to perform FM demodulation, they'd work too. So far I've summarized the two methods I'm aware of above. _____________________________ Posted through www.DSPRelated.com
Reply by January 24, 20152015-01-24
On Friday, January 23, 2015 at 9:22:03 AM UTC+13, all4dsp wrote:
> I've got a real voltage versus time waveform (e.g. xr) of a carrier > frequency that is FM modulated. I need to see the modulation waveform to > determine its peak-to-peak amplitude. While I do see a lot of code online > for FM demodulation, I'm hesitant to use it unless I can understand it from > first principles. The method that makes most sense to me is the equation I > see in Rick Lyon's excellent book "Understanding Digital Signal Processing" > 2nd Edition (see equation 9-6 or 13-11): > > F(t) = d (phi(t)) / dt = d ( atan( xi / xr ) ) / dt > > where xi is Hilbert transform of input signal xr. My Matlab code to > implement this is: > > dt = 1/Fs; % Fs=sample rate > ht = hilbert( xr ); % Hilbert transform of real signal Vout > phi = atan( imag(ht) ./ xr ); % take arctan > freq = ( phi(2:end) - phi(1:end-1) ) / dt; % compute instantaneous > freq by derivative of phi > freq_hz = freq/(2*pi); % convert rad/sec to Hz > > However, it seems there's a problem when Matlab computes atan() in that > there are discrete phase jumps when phi crosses the pi/2 boundary to > re-appear at -pi/2. I tried computing phi=unwrap(...) but haven't had any > success removing these jumps. > > As I assume FM demodulation is pretty popular, I was hoping someone could > share with me how to eliminate these phase jumps. Thanks in advance for any > help. > > _____________________________ > Posted through www.DSPRelated.com
Perhaps you can tell us why, when 99.9% of all radio receivers don't use a Hilbert Transform, why you feel you need one here. I assume you are trying to recover phase. The usual method is to use a PLL on the carrier or IF signal - or a discriminator - quad detectors are normally the best though due to phase noise of PLLs. Also, if you want to recover phase then FFTs can be used using complex cepstrum instead of Hilbert Transforms.
Reply by Rick Lyons January 24, 20152015-01-24
On Fri, 23 Jan 2015 19:09:25 -0600, "all4dsp" <51070@dsprelated>
wrote:

>If I obtain FM demodulated results from this code: > > dt = 1/Fs; % Fs=sample rate > ht = hilbert( xr ); > phi = atan2( imag(ht), real(ht) ); % take arctan > freq = ( phi(2:end) - phi(1:end-1) ) / dt; % compute instantaneous >freq by derivative of phi > freq_hz = freq/(2*pi); % convert rad/sec to Hz > >and compare it with results using this code: > > dt = 1/Fs; % Fs=sample rate > ht = hilbert( xr ); > yd = ht( 2:length(ht) ) .* conj( ht(1:length(ht)-1) ); > freq_hz=angle(yd)/(2*pi*dt); % units of Hz > >Their FFT spectrums do not match.
Hi, In your first snippet of code above, if you "unwrap" the phase of the 'phi' sequence before performing the differentiation then your two spectra will better match each other. [-Rick-]
Reply by Rick Lyons January 24, 20152015-01-24
On Fri, 23 Jan 2015 19:09:25 -0600, "all4dsp" <51070@dsprelated>
wrote:

>If I obtain FM demodulated results from this code: > > dt = 1/Fs; % Fs=sample rate > ht = hilbert( xr ); > phi = atan2( imag(ht), real(ht) ); % take arctan > freq = ( phi(2:end) - phi(1:end-1) ) / dt; % compute instantaneous >freq by derivative of phi > freq_hz = freq/(2*pi); % convert rad/sec to Hz > >and compare it with results using this code: > > dt = 1/Fs; % Fs=sample rate > ht = hilbert( xr ); > yd = ht( 2:length(ht) ) .* conj( ht(1:length(ht)-1) ); > freq_hz=angle(yd)/(2*pi*dt); % units of Hz > >Their FFT spectrums do not match. > >Assuming the first code's results are correct (because I can understand the >code from first principles), does anyone know why the 2nd code doesn't >work? I *think* it is supposed to do the same thing (e.g. FM demodulate), >but then again I found it online and can't explain it from first >principles. Can anyone shed some light on the 2nd code's approach and why >it should (or shouldn't) result in the same thing as the first code?
Hi all4dsp, As far as I know that "2nd code" technique originated in a 1995 Masters thesis by Jim Shima. I'll send you a copy of that thesis later today. [-Rick-]
Reply by all4dsp January 23, 20152015-01-23
If I obtain FM demodulated results from this code:

    dt = 1/Fs;  % Fs=sample rate
    ht = hilbert( xr );  
    phi = atan2( imag(ht), real(ht) );  % take arctan
    freq = ( phi(2:end) - phi(1:end-1) ) / dt;  % compute instantaneous
freq by derivative of phi
    freq_hz = freq/(2*pi);  % convert rad/sec to Hz
    
and compare it with results using this code:
    
    dt = 1/Fs;  % Fs=sample rate
    ht = hilbert( xr );  
    yd = ht( 2:length(ht) ) .* conj( ht(1:length(ht)-1) );
    freq_hz=angle(yd)/(2*pi*dt);  % units of Hz

Their FFT spectrums do not match. 

Assuming the first code's results are correct (because I can understand the
code from first principles), does anyone know why the 2nd code doesn't
work? I *think* it is supposed to do the same thing (e.g. FM demodulate),
but then again I found it online and can't explain it from first
principles. Can anyone shed some light on the 2nd code's approach and why
it should (or shouldn't) result in the same thing as the first code?	 

_____________________________		
Posted through www.DSPRelated.com
Reply by January 22, 20152015-01-22
On Friday, January 23, 2015 at 9:22:03 AM UTC+13, all4dsp wrote:
> I've got a real voltage versus time waveform (e.g. xr) of a carrier > frequency that is FM modulated. I need to see the modulation waveform to > determine its peak-to-peak amplitude. While I do see a lot of code online > for FM demodulation, I'm hesitant to use it unless I can understand it from > first principles. The method that makes most sense to me is the equation I > see in Rick Lyon's excellent book "Understanding Digital Signal Processing" > 2nd Edition (see equation 9-6 or 13-11): > > F(t) = d (phi(t)) / dt = d ( atan( xi / xr ) ) / dt > > where xi is Hilbert transform of input signal xr. My Matlab code to > implement this is: > > dt = 1/Fs; % Fs=sample rate > ht = hilbert( xr ); % Hilbert transform of real signal Vout > phi = atan( imag(ht) ./ xr ); % take arctan > freq = ( phi(2:end) - phi(1:end-1) ) / dt; % compute instantaneous > freq by derivative of phi > freq_hz = freq/(2*pi); % convert rad/sec to Hz > > However, it seems there's a problem when Matlab computes atan() in that > there are discrete phase jumps when phi crosses the pi/2 boundary to > re-appear at -pi/2. I tried computing phi=unwrap(...) but haven't had any > success removing these jumps. > > As I assume FM demodulation is pretty popular, I was hoping someone could > share with me how to eliminate these phase jumps. Thanks in advance for any > help. > > _____________________________ > Posted through www.DSPRelated.com
You don't need atan at all. Analytically if you differentiate y/x you get 1/(1+(y/x)^2)X d/t(y/x) - something like that anyway.
Reply by Rick Lyons January 22, 20152015-01-22
On Thu, 22 Jan 2015 16:21:16 -0600, "all4dsp" <51070@dsprelated>
wrote:

>> >>I don't speak MATLAB but a standard method >>for avoiding jumps is to use atan2(x,y) >>instead of atan(y/x).
Hi all4dsp, Dilip's (dvsarwate) suggestion is a good one.
>Correct me if I'm wrong, but I don't think I can substitute x=xi and y=xr >here to compute Hilbert transform.
When Dilip suggests using 'atan2(x,y)' he's not referring to computing a Hilbert transform. He's referring to computing the instananeous angle of a complex-valued data sequence when the desired angle lies anywhere between -pi and +pi radians.
>That is, I "think" the x in atan(y/x) is the x-axis data point (e.g. time >in units of seconds for me), and y is y-axis point (e.g .voltage in units >of volts for me). I thought atan2 was only for conversion from/to polar >coordinates. Let me know if I'm way off.
Whoa, wait a moment. Dilip's 'x' is the real part (in-phase part) of a complex-valued amplitude sample and his 'y' is the imaginary part (quadrature phase part) of that complex-valued amplitude sample. His 'x' is your 'xr' and his 'y' is your 'xi'. If you want to avoid the potential pitfalls of computing arctangents you might consider using my 2nd Edition's Equation (13-116). That particular FM demod method, which sadly I did not invent, works pretty well. If you're able to send me a private e-mail with the files that allow me to generate your 'xr' FM input signal and duplicate your Matlab code, I will run your code and try to help. [-Rick-]
Reply by all4dsp January 22, 20152015-01-22
> >I don't speak MATLAB but a standard method >for avoiding jumps is to use atan2(x,y) >instead of atan(y/x). >
Correct me if I'm wrong, but I don't think I can substitute x=xi and y=xr here to compute Hilbert transform. That is, I "think" the x in atan(y/x) is the x-axis data point (e.g. time in units of seconds for me), and y is y-axis point (e.g .voltage in units of volts for me). I thought atan2 was only for conversion from/to polar coordinates. Let me know if I'm way off. _____________________________ Posted through www.DSPRelated.com
Reply by dvsarwate January 22, 20152015-01-22
On Thursday, January 22, 2015 at 2:22:03 PM UTC-6, all4dsp wrote:
> I've got a real voltage versus time waveform (e.g. xr) of a carrier > frequency that is FM modulated. I need to see the modulation waveform to > determine its peak-to-peak amplitude. While I do see a lot of code online > for FM demodulation, I'm hesitant to use it unless I can understand it from > first principles. The method that makes most sense to me is the equation I > see in Rick Lyon's excellent book "Understanding Digital Signal Processing" > 2nd Edition (see equation 9-6 or 13-11): > > F(t) = d (phi(t)) / dt = d ( atan( xi / xr ) ) / dt > > where xi is Hilbert transform of input signal xr. My Matlab code to > implement this is: > > dt = 1/Fs; % Fs=sample rate > ht = hilbert( xr ); % Hilbert transform of real signal Vout > phi = atan( imag(ht) ./ xr ); % take arctan > freq = ( phi(2:end) - phi(1:end-1) ) / dt; % compute instantaneous > freq by derivative of phi > freq_hz = freq/(2*pi); % convert rad/sec to Hz > > However, it seems there's a problem when Matlab computes atan() in that > there are discrete phase jumps when phi crosses the pi/2 boundary to > re-appear at -pi/2. I tried computing phi=unwrap(...) but haven't had any > success removing these jumps. > > As I assume FM demodulation is pretty popular, I was hoping someone could > share with me how to eliminate these phase jumps. Thanks in advance for any > help. > > _____________________________ > Posted through www.DSPRelated.com
I don't speak MATLAB but a standard method for avoiding jumps is to use atan2(x,y) instead of atan(y/x).