# [ncl-talk] The question about the complex demodulation

Dennis Shea shea at ucar.edu
Fri Dec 15 13:56:21 MST 2017

```This is a rather long and, perhaps, confusing response. It was written
episodically.

---
re: *"same mistake" *   .... I am not sure to what you are referring. The
use of a minus sign?

The minus sign does not affect the amplitudes because the quantities are
squared:  amp = sqrt(x*x+y*y)
It does affect the phase. Also, use of *atan2 *or *atan *to determine the
phase are other issues to consider.

https://www.ncl.ucar.edu/Document/Functions/Built-in/atan2.shtml
https://www.ncl.ucar.edu/Document/Functions/Built-in/atan.shtml

PMEL:   https://www.pmel.noaa.gov/maillists/tmap/ferret_users/
fu_2007/pdfHMANlJLqCE.pdf
Kessler: https://faculty.washington.edu/kessler/generals/complex-
demodulation.pdf

PMEL references the same source as NCL.
Reference: Bloomfield, P. (1976). Fourier Decomposition of Time Series, An
Introduction. Wiley, New York, 258pp.

===
I have attached the (slightly updated) fortran subroutines from Bloomfield
but will include here also::
*demod_bloom1976 *(page 148)
*polar_bloom1976* (page 150)

---------------------------
subroutine *demod_bloom1976*(x,n,omega,d1,d2)
implicit none
c
c Fourier Analysis of Time Series: An Introduction
c P. Bloomfield (1976: Wiley); *p148*
c Minor code changes to make code more 'modern'
c
integer n                        ! INPUT
real    x(n), omega, d1(n), d2(n)
integer i                        ! LOCAL
real    arg

do i=1,n
arg   = (i-1)*omega

*         d1(i) =  x(i)*cos(arg)*2         d2(i) = -x(i)*sin(arg)*2*
end do

return
end
c ---------------------------
subroutine *polar_bloom1976*(x,y,n)
implicit none
c
c Fourier Analysis of Time Series: An Introduction
c P. Bloomfield (1976: Wiley); *p150*
c Minor code changes to make code more 'modern'
c
integer n                        ! INPUT
real    x(n), y(n)
integer i                        ! LOCAL
real    amp, phi

do i=1,n

*         amp = sqrt(x(i)*x(i)+y(i)*y(i))         phi =
x(i)= amp
y(i)= phi
end do

return
end
c ---------------------------
c Bloomfield CALL SEQUENCE
c
c omega  = (2*pi*nom)/np2
c cutoff = (pi*(npa+nst))/np2
c ns     = np2/(nst-npa)
c call demod_bloom1976 (x,n,omega,d1,d2) ! DEMODULATION
c call lopass_bloom1976(d1,cutoff,ns)    ! SMOOTH COSINE TERM
c call lopass_bloom1976(d2,cutoff,ns)    ! SMOOTH SINE TERM
c lim = n-2*ns
c call polar_bloom1976(d1,d2,lim)        ! AMPLITUDE & PHASES
===============================
*Bloomfield's book uses a different low-pass filter than NCL (Lanczos) so
results are not directly ('exactly') comparable.*
===============================
NCL's  (attached) *demod_cmplx* uses array notation (no 'do' loops):

...
yr  =  y*cos(ARG)*2             ; 'real' series
yi  = *-*y*sin(ARG)*2              ; 'imaginary'
...
wgt    = filwgts_lanczos (nwt, ihp, frqcut, fcb, nsigma)

;---Apply filter weights to raw demodulated series (nwt/2 point 'lost; at
start/end)

yr = wgt_runave_n_Wrap ( yr, wgt, 0, 0)
yi = wgt_runave_n_Wrap ( yi, wgt, 0, 0)

;---Use the low pass filtered values; *atan2 *uses y/x => yi/yr

amp  = sqrt(yr^2 + yi^2)
pha  = where(amp.eq.zero, zero, *atan2*(yi, yr))    ; primary phase

pha2 = conform(pha, zero, -1)                     ; secondary phase
(Bloomfield)
pha2 = where(pha.ge.zero, pha-pi2, pha+pi2)

==================================
You could take NCL's attached '*demod_cmplx*'; copy it to create your own
function (say):

*demod_cmplx_nai*
and change what you want .... for example.... remove the minus sign

yr  =  y*cos(ARG)*2             ; 'real' series
yi  =   y*sin(ARG)*2              ; 'imaginary'

Play with it.

==================================
I did not do anything 'fancy' with the phase plot. Hence, there can be
'jumps' . Think 358 degrees to 2 degrees. 'We' know it is only 4 degrees
different but the plotting does not. Hence, a non-pretty jump in the plot.

===================================
Archived National Institute of Standards fortran software.

*DEMODU.f*  (attached): Does the same thing as Bloomfield/NCL.
The use af AMPL and PHAS names is because this was part of a larger library
that reused array space.

The following

http://www.itl.nist.gov/div898/software/datapac/DEMOD.f

uses

Y1(I)=X(I)*COS(6.2831853*F*AI)
Y2(I)=X(I)*SIN(6.2831853*F*AI)

but then for phases uses *atan *(not *atan2*)

Z(I)=*ATAN*(Y1(I)/Y2(I))     !  atan(yr/yi)*  ... ?*implicitly
handles - sign?

I am done.

Cheers
D

On Wed, Dec 13, 2017 at 8:31 AM, Nai Suan <shif422 at gmail.com> wrote:

> Dear Dennis Shea,
>         I found a same mistake in two documents Kessler: U. Washington:
> Complex Demodulation and PMEL: Complex Demodulation. The first formula
> should be conected by a ‘plus’ sign not ’minus’ sign.
>         All the best,
> Feng
>
>
>
> > On Dec 13, 2017, at 10:02 AM, Nai Suan <shif422 at gmail.com> wrote:
> >
> > Dear Dennis Shea,
> >       If I am right, the phases of the 1700-1960 sunspot numbers has a
> unique solution. This means Bloomfield’s 2000 (Second Edition) is wrong or
> the 1st complex demodulation example in NCL is wrong. I have no idea for
> another option.
> >       You mentioned NCL example was translated from Bloomfield’s
> subroutine. Is there any possible reason that the translation has some
> problems?
> >       Anyway, if you are interested, I can send Bloomfield’s 2000 book
> to you.
> >       All the best,
> > Frank
> >
> >
> >> On Dec 13, 2017, at 12:15 AM, Dennis Shea <shea at ucar.edu> wrote:
> >>
> >> If the low-pass filtering (Lanczos) and smoothing (weighted running
> average) don't explain the differences... I do not know.
> >> I do not have Bloomfield's 1976 book readily available. That used the
> 1700-1960 sunspot numbers. As I recall, NCL's plots were 'identical' to
> those in the 1976 book.
> >> ----
> >> http://www.ncl.ucar.edu/Applications/spec.shtml
> >>
> >> Data Analysis: Spectral analysis, Complex Demodulation
> >>
> >> The 1st complex demodulation example uses the 1700-1960 sunspots
> >> http://www.ncl.ucar.edu/Applications/Scripts/demod_cmplx_1.ncl
> >> ----
> >>
> >> The 'demod_cmplx' function can be seen at:
> >>
> >> %> less  \$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl
> >>
> >> Basically, it is a translation (fortran-77  ===> NCL) of Bloomfield's
> 1976 f77 subroutine.
> >>
> >>
> >>
> >>
> >>
> >>
> >> On Tue, Dec 12, 2017 at 4:03 AM, Nai Suan <shif422 at gmail.com> wrote:
> >> Dear colleague,
> >>         I have a question about how to explain the phase of the complex
> demodulation. I found that the phases of the sunspot in the figures 7.9 and
> 7.11 (see the attachment) in Bloomfield (2000)’s book (chapter 7 complex
> demodulation) were totally different with the ‘demod_cmplx_1_lg.png (see
> the attachment)’. Which is right? I believe the original sunspots in two
> examples are same and the different low-pass filter methods can not be used
> to explain this difference. How to explain the difference?
> >>         All the best,
> >> Frank
> >>
> >>
> >>
> >> >
> >> >
> >> >
> >> >
> >> >
> >>
> >>
> >> _______________________________________________
> >> ncl-talk mailing list
> >> ncl-talk at ucar.edu
> >> List instructions, subscriber options, unsubscribe:
> >> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
> >
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20171215/3e7d9899/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: demod_cmplx.ncl
Type: application/octet-stream
Size: 3537 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20171215/3e7d9899/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: demod_bloomfield_1976.f
Type: application/octet-stream
Size: 1635 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20171215/3e7d9899/attachment-0001.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: DEMODU_nist.f
Type: application/octet-stream
Size: 2112 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20171215/3e7d9899/attachment-0002.obj>
```