[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 =
atan2(y(i),x(i)) *! radians
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
(radians)
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,
> > Thank you for your reply.
> > 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?
> >> Hope to hear your answer.
> >> 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>
More information about the ncl-talk
mailing list