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

# Classic FM Demodulation, using Matlab

Started by ●January 22, 2015

Reply by ●January 22, 20152015-01-22

>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 itfrom>first principles. The method that makes most sense to me is the equationI>see in Rick Lyon's excellent book "Understanding Digital SignalProcessing">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 forany>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

Reply by ●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.comI don't speak MATLAB but a standard method for avoiding jumps is to use atan2(x,y) instead of atan(y/x).

Reply by ●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 ●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 ●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.comYou 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 ●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 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 ●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 ●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.comPerhaps 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.