Hi all. I've tried to implement an FM demodulation scheme, as suggested by Rick in ch. 13.22 of his 2nd edition (MATLAB code attached at the end of the post). I have based my approach on Rick's eq. 13-117, but I can't get it to work: - When I don't include the denominator i(n)^2+q(n)^2, I get an instantaneous phase estimate that's close to the "correct" one, but with a sign error. - When I include the denominator there is completely destructable ripple in the instantaneous phase estimate DeltaTheta(n) Any help in spotting blunders, errors and flaws is highly appreciated. Rune %%%%%%%%%%%%%%%%%% MATLAB Script %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Sampling parameters: fs=1; % Sampling frequency [Hz] f0=0.05*fs; % Start frequency for linear sweep df=0.0025*f0; % sweep rate df/dt Ns=128; % Length of sweeped signal [samples] % Generate sweeped signal t=reshape([0:Ns-1]/fs,Ns,1); % Time variable in sweep frame d=sin((f0+df*t)*2*pi.*t); % Generate sweep N=1024; % Length of total signal N0=100; % Start index of sweeped frame % x : raw data % tv : time vector for raw data x=zeros(N,1); x(N0:N0+Ns-1)=d; tv=reshape([0:N-1]/fs,N,1); % Design LP filter for band limitation and syncing delays Nf=128; % Length of filter fc=0.4*fs; % Cut-off of filter hr=sinc(reshape([-(Nf-1)/2:(Nf-1)/2],Nf,1)*2*pi*fc).*kaiser(Nf,0.7); xr=filter(hr,1,x); % Real part of data - inphase part % Design Hilbert transformer - Correct? hi=remez(Nf,[0.025,0.975],[1 1],'hilbert'); xi=filter(hi,1,x); % Imaginary part of data - quadrature part %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% xr=xr-xi; % Eh... correct? %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Differentiator, Ref Lyons p 552 Nd=3; diff =[1,0,-1]; % Compute differentials of xr(n) and xi(n) as % indicated by Lyons, eq. (13-117) dxr=filter(diff,1,xr); dxi=filter(diff,1,xi); DeltaTheta=zeros(length(xi),1); idx=find((xr.^2+xi.^2)>10e-12); DeltaTheta(idx)=(xr(idx).*dxi(idx)-xi(idx).*dxr(idx)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% BAD, BAD, BAD.... %DeltaTheta(idx)=DeltaTheta(idx)./(xr(idx).^2+xi(idx).^2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Adjust time scale after having applied filters ttv=reshape([0:length(xi)-1]/fs,length(xi),1); clf subplot(211) % Plot complex time series plot(ttv,xr,'b',ttv,xi,'r') subplot(212) % Plot estimated instantaneous frequency plot(ttv,DeltaTheta/(2*pi)) hold on % Actual instantanous frequency plot(t+N0,f0+df*t,'r') % Instantaneous frequency delayed after filters and differentiators plot(t+N0+Nf/2+Nd-2,f0+df*t,'g') % Bring figure to the foreground figure(gcf)
FM Demodulation problem (Lyons Ch. 13.22)
Started by ●September 19, 2004
Reply by ●September 20, 20042004-09-20
On 19 Sep 2004 08:01:03 -0700, allnor@tele.ntnu.no (Rune Allnor) wrote:>Hi all. > >I've tried to implement an FM demodulation scheme, as suggested >by Rick in ch. 13.22 of his 2nd edition (MATLAB code attached >at the end of the post). > >I have based my approach on Rick's eq. 13-117, but I can't get >it to work: > >- When I don't include the denominator i(n)^2+q(n)^2, I get an > instantaneous phase estimate that's close to the "correct" one, > but with a sign error. >- When I include the denominator there is completely destructable > ripple in the instantaneous phase estimate DeltaTheta(n) > >Any help in spotting blunders, errors and flaws is highly appreciated. > >Rune >Hi Rune, I just saw your post. I think you FM demodulator is working (almost) fine. If you use the following ("xi = -xi") command, before performing the derivative, to change the phase of the imaginary part by 180 degrees, then the final results show the signal's frequency increasing with time (instead of decreasing with time). xi = -xi; Below is your code with a few changes by me. Please try t run my code. My changes to your code are the "indented" lines. As for the ripples in the FM demod output, aren't those ripples merely the transient response of the 128-tap filters? I could be wrong. Ya' know what seems weird to me? Your "xr=xr-xi" command does not look like the right thing to do. However, ... your FM demodulator seems to work properly whether you use that "xr=xr-xi" command or not!! I can't (off the top of my head) explain that. [-Rick-] % % Filename:FM_Demod_RuneAllnor.m %%%%%%%%%%%%%%%%%% MATLAB Script %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear % Added by Lyons % Sampling parameters: fs=1; % Sampling frequency [Hz] f0=0.05*fs; % Start frequency for linear sweep df=0.0025*f0; % sweep rate df/dt Ns=128; % Length of sweeped signal [samples] % Generate sweeped signal t=reshape([0:Ns-1]/fs,Ns,1); % Time variable in sweep frame d=sin((f0+df*t)*2*pi.*t); % Generate sweep N=1024; % Length of total signal N0=100; % Start index of sweeped frame % x : raw data % tv : time vector for raw data x=zeros(N,1); x(N0:N0+Ns-1)=d; tv=reshape([0:N-1]/fs,N,1); % Design LP filter for band limitation and syncing delays Nf=128; % Length of filter fc=0.4*fs; % Cut-off of filter hr=sinc(reshape([-(Nf-1)/2:(Nf-1)/2],Nf,1)*2*pi*fc).*kaiser(Nf,0.7); xr=filter(hr,1,x); % Real part of data - inphase part % Design Hilbert transformer - Correct? hi=remez(Nf,[0.025,0.975],[1 1],'hilbert'); xi=filter(hi,1,x); % Imaginary part of data - quadrature part % Following added by Lyons. % Compare phase of real & imaginary % parts of the analytic signal "x". figure(10), clf plot(125:325,xr(125:325),'-b') hold on plot(125:325,xi(125:325),'-r'), grid on, zoom on hold off title('Analytic "x", blue=real, red=imag.') % If you use the following "xi = -xi" command, % to change the phase of the imaginary part by % 180 degrees, then the final results show the % signal's frequency increasing with time (instead % of decreasing with time). xi = -xi; % Added by Lyons %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Rune, the following "xr=xr-xi" command % does not look like the right thing to do. % However, your FM demodulator seems to % work properly whether you use that % "xr=xr-xi" command or not!! xr=xr-xi; % Eh... correct? %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Differentiator, Ref Lyons p 552 Nd=3; diff =[1,0,-1]; % Compute differentials of xr(n) and xi(n) as % indicated by Lyons, eq. (13-117) dxr=filter(diff,1,xr); dxi=filter(diff,1,xi); % Following added by Lyons. % Compare phase of the derivatives of % the real & imaginary parts of the % analytic signal "x". figure(20), clf plot(125:325,dxr(125:325),'-b'), grid on, zoom on hold on plot(125:325,dxi(125:325),'-r'), grid on, zoom on hold off title('Derivative of Analytic "x", blue=real, red=imag.') DeltaTheta=zeros(length(xi),1); idx=find((xr.^2+xi.^2)>10e-12); DeltaTheta(idx)=(xr(idx).*dxi(idx)-xi(idx).*dxr(idx)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% BAD, BAD, BAD.... %DeltaTheta(idx)=DeltaTheta(idx)./(xr(idx).^2+xi(idx).^2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Adjust time scale after having applied filters ttv=reshape([0:length(xi)-1]/fs,length(xi),1); figure(1) % % Added by Lyons clf subplot(211) % Plot complex time series plot(ttv,xr,'b',ttv,xi,'r') grid on, zoom on, % Added by Lyons % Following added by Lyons text(340,1.2,'Is not the imaginary part (red) out') text(340,0.7,'of phase from what is should be?') text(340,0.2,'Zoom in to see what I mean.') text(340,-0.3,'Should not the Hilbert xform') text(340,-0.8,'of the blue cosine be a sinewave?') subplot(212) % Plot estimated instantaneous frequency plot(ttv,DeltaTheta/(2*pi)) hold on % Actual instantanous frequency plot(t+N0,f0+df*t,'r') % Instantaneous frequency delayed after filters and differentiators plot(t+N0+Nf/2+Nd-2,f0+df*t,'g') grid on, zoom on, % Lyons % Bring figure to the foreground figure(gcf)
Reply by ●September 21, 20042004-09-21
r.lyons@_BOGUS_ieee.org (Rick Lyons) wrote in message news:<414ec836.391447765@news.sf.sbcglobal.net>...> On 19 Sep 2004 08:01:03 -0700, allnor@tele.ntnu.no (Rune Allnor) > wrote: > > >Hi all. > > > >I've tried to implement an FM demodulation scheme, as suggested > >by Rick in ch. 13.22 of his 2nd edition (MATLAB code attached > >at the end of the post). > > > >I have based my approach on Rick's eq. 13-117, but I can't get > >it to work: > > > >- When I don't include the denominator i(n)^2+q(n)^2, I get an > > instantaneous phase estimate that's close to the "correct" one, > > but with a sign error. > >- When I include the denominator there is completely destructable > > ripple in the instantaneous phase estimate DeltaTheta(n) > > > >Any help in spotting blunders, errors and flaws is highly appreciated. > > > >Rune > > > > Hi Rune, I just saw your post. > > I think you FM demodulator is working > (almost) fine. > > If you use the following ("xi = -xi") command, > before performing the derivative, > to change the phase of the imaginary part by > 180 degrees, then the final results show the > signal's frequency increasing with time (instead > of decreasing with time). > > xi = -xi;OK, thanks.> Below is your code with a few changes by me. > Please try t run my code. > My changes to your code are the "indented" lines. > > As for the ripples in the FM demod output, aren't > those ripples merely the transient response of the > 128-tap filters? I could be wrong.I think you could well be right. It doesn't seem to matter, though, since they appear to be end effects and nothing more. If upi uncomment the "BAD" line, you'll see ripple that appears to be correlated with the zero-crossings of the frequency sweep.> Ya' know what seems weird to me? Your "xr=xr-xi" command > does not look like the right thing to do. > However, ... your FM demodulator seems to > work properly whether you use that > "xr=xr-xi" command or not!! I can't (off the top of > my head) explain that.You commented somewhere that the Hilbert transform doesn't seem to be of correct phase. If you comment out the "xr=xr-xi;" line, the phase relation seems to get more like what it should be.> > [-Rick-]Thanks for your help. I�ll have to look into this some more. Rune> % > % Filename:FM_Demod_RuneAllnor.m > > %%%%%%%%%%%%%%%%%% MATLAB Script %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > > > clear % Added by Lyons > > % Sampling parameters: > > fs=1; % Sampling frequency [Hz] > f0=0.05*fs; % Start frequency for linear sweep > df=0.0025*f0; % sweep rate df/dt > Ns=128; % Length of sweeped signal [samples] > > % Generate sweeped signal > t=reshape([0:Ns-1]/fs,Ns,1); % Time variable in sweep frame > d=sin((f0+df*t)*2*pi.*t); % Generate sweep > > N=1024; % Length of total signal > N0=100; % Start index of sweeped frame > > > % x : raw data > % tv : time vector for raw data > x=zeros(N,1); > x(N0:N0+Ns-1)=d; > tv=reshape([0:N-1]/fs,N,1); > > % Design LP filter for band limitation and syncing delays > Nf=128; % Length of filter > fc=0.4*fs; % Cut-off of filter > hr=sinc(reshape([-(Nf-1)/2:(Nf-1)/2],Nf,1)*2*pi*fc).*kaiser(Nf,0.7); > > xr=filter(hr,1,x); % Real part of data - inphase part > > % Design Hilbert transformer - Correct? > hi=remez(Nf,[0.025,0.975],[1 1],'hilbert'); > > xi=filter(hi,1,x); % Imaginary part of data - quadrature part > > % Following added by Lyons. > % Compare phase of real & imaginary > % parts of the analytic signal "x". > figure(10), clf > plot(125:325,xr(125:325),'-b') > hold on > plot(125:325,xi(125:325),'-r'), grid on, zoom on > hold off > title('Analytic "x", blue=real, red=imag.') > > % If you use the following "xi = -xi" command, > % to change the phase of the imaginary part by > % 180 degrees, then the final results show the > % signal's frequency increasing with time (instead > % of decreasing with time). > > xi = -xi; % Added by Lyons > > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > > % Rune, the following "xr=xr-xi" command > % does not look like the right thing to do. > % However, your FM demodulator seems to > % work properly whether you use that > % "xr=xr-xi" command or not!! > > xr=xr-xi; % Eh... correct? > > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > > % Differentiator, Ref Lyons p 552 > Nd=3; > diff =[1,0,-1]; > > % Compute differentials of xr(n) and xi(n) as > % indicated by Lyons, eq. (13-117) > dxr=filter(diff,1,xr); > dxi=filter(diff,1,xi); > > % Following added by Lyons. > % Compare phase of the derivatives of > % the real & imaginary parts of the > % analytic signal "x". > figure(20), clf > plot(125:325,dxr(125:325),'-b'), grid on, zoom on > hold on > plot(125:325,dxi(125:325),'-r'), grid on, zoom on > hold off > title('Derivative of Analytic "x", blue=real, red=imag.') > > > DeltaTheta=zeros(length(xi),1); > idx=find((xr.^2+xi.^2)>10e-12); > DeltaTheta(idx)=(xr(idx).*dxi(idx)-xi(idx).*dxr(idx)); > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > %%%%% BAD, BAD, BAD.... > %DeltaTheta(idx)=DeltaTheta(idx)./(xr(idx).^2+xi(idx).^2); > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > > % Adjust time scale after having applied filters > ttv=reshape([0:length(xi)-1]/fs,length(xi),1); > > figure(1) % % Added by Lyons > clf > subplot(211) > % Plot complex time series > plot(ttv,xr,'b',ttv,xi,'r') > grid on, zoom on, % Added by Lyons > > % Following added by Lyons > text(340,1.2,'Is not the imaginary part (red) out') > text(340,0.7,'of phase from what is should be?') > text(340,0.2,'Zoom in to see what I mean.') > text(340,-0.3,'Should not the Hilbert xform') > text(340,-0.8,'of the blue cosine be a sinewave?') > > subplot(212) > % Plot estimated instantaneous frequency > plot(ttv,DeltaTheta/(2*pi)) > hold on > > % Actual instantanous frequency > plot(t+N0,f0+df*t,'r') > > % Instantaneous frequency delayed after filters and differentiators > plot(t+N0+Nf/2+Nd-2,f0+df*t,'g') > grid on, zoom on, % Lyons > > % Bring figure to the foreground > figure(gcf)
