[ncl-talk] The question about the complex demodulation
Nai Suan
shif422 at gmail.com
Sun Dec 17 02:47:44 MST 2017
Dear Dennis Shea,
Thank you for your huge effect.
I agree with you, a minus sign does not affect the amplitude but affect the phase.
If I understand well, you gave me two reasons to explain the difference of the phases in the Bloomfield’s book and the NCL example.
1) you show that the ‘atan2' or ‘atan’ functions may effect the phase plot. I do not agree with you. It can not be resolved the ‘jump’ in the phase plot.
2) you think that the different lowpass filters may be a reason. I do not agree with you too. I have used eight filter methods to carry out the complex demodulation. The results show that the filter method can slightly affect the amplitude and phase, but the filter method can not be used to explain the ‘jump’ in phase plot.
Anyway, there is no an obvious difference in amplitude plot between Bloomfield’s book and NCL example. My question still focus on the differences of phases using the complex demodulation between Bloomfield’s book and NCL example. How to explain the huge and obvious differences? In a word, How to deal with the obvious ‘jump’ in NCL example? The NCL code ‘demod_cmplx_nai’ can not deal with the ‘jump’ too.
You implied that a 'fancy' with the phase plot in Bloomfield’s book? Do you know or guess what technique Bloomfield did to deal with "a non-pretty jump”? I agree that Bloomfield used a technique to show only 4 degrees difference from 358 degrees to 2 degrees, but the technique does not include in the two Fortran subroutines which you sent to me. Do you have any idea to reproduce the phase plot in Bloomfield’s book?
All the best,
Frank
> On Dec 15, 2017, at 9:56 PM, Dennis Shea <shea at ucar.edu> wrote:
>
> 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/atan2.shtml>
> https://www.ncl.ucar.edu/Document/Functions/Built-in/atan.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 <https://www.pmel.noaa.gov/maillists/tmap/ferret_users/fu_2007/pdfHMANlJLqCE.pdf>
> Kessler: https://faculty.washington.edu/kessler/generals/complex-demodulation.pdf <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 <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 <mailto: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 <mailto: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 <mailto: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 <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 <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 <mailto: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 <mailto:ncl-talk at ucar.edu>
> >> List instructions, subscriber options, unsubscribe:
> >> http://mailman.ucar.edu/mailman/listinfo/ncl-talk <http://mailman.ucar.edu/mailman/listinfo/ncl-talk>
> >
>
>
> <demod_cmplx.ncl><demod_bloomfield_1976.f><DEMODU_nist.f>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20171217/dfdf9299/attachment.html>
More information about the ncl-talk
mailing list