Classic FM Demodulation, using Matlab

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'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. 	

I found the following sequence to eliminate concerns about phase jumps:

    dt = 1/Fs;  % Fs=sample rate
    ht = hilbert( xr );  % Hilbert transform of real signal Vout
    yd = ht( 2:length(ht) ) .* conj( ht(1:length(ht)-1) );
    freq_hz=angle(yd)/(2*pi*dt);  % units of Hz

although I'm not exactly sure about the theory. 

When I plot freq_hz versus time, I noticed some negative amplitude data
points occurring periodically. Could someone confirm that this algorithm is
indeed supposed to plot the instantaneous frequency versus time of signal
xr? If so, how to explain negative values, or do I have something wrong
here?	 

_____________________________		
Posted through www.DSPRelated.com

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).

>
>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

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-]

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.

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

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-]

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-]

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.