[ncl-talk] The question about the complex demodulation

Dennis Shea shea at ucar.edu
Mon Jan 1 12:41:38 MST 2018


[1] frqdem = period/tofloat(ntim)
        The demodulation frequency should equal to 1/period.

Yes ... I should have coded what I wrote.   :-(

===
[2] frqcut = 0.010
        The cut-off frequency should equal to 0.5*frqdem

No. There is relation between frqdem and frqcut.
The latter is strictly related to the desired low-pass filter
characteristics.
===
[3]   nwt    = 41

This was used only because Bloomfield (1976) used that number  of weights
for his least squares low pass filter.
NCL's 'demod_cmplx' uses a different low-pass filter (filwgts_lanczos). I
have no idea how they compare in their response functions. Likely, they are
a bit different. That is why getting exact results is problematical.

++++++++++++++++
I updated the spectral analysis and complex-demodulation page:

https://www.ncl.ucar.edu/Applications/spec.shtml

The sunspot demodulation example shows two figures with different cutoff
frequencies. All other parameters are the same.
+++++++++++++++++

re: Moreover, the updated one can not correctly run in NCL 6.4.0.

The 'unwrap_phase' I attached used a new feature of NCL: a 'true'


*elseif    https://www.ncl.ucar.edu/future_release.shtml
<https://www.ncl.ucar.edu/future_release.shtml>*

===
The unwrap_phase_matlab.ncl could be used.

I have attached a version of unwrap_phase that should run with 6.4.0

Good Luck
D


On Mon, Jan 1, 2018 at 12:43 AM, Nai Suan <shif422 at gmail.com> wrote:

> Dear Dennis Shea,
>         Thank you for your updating.
>         Some assignments are still obscure in the updated one, which cause
> that the results are still different with Bloomfield’s book.
>         1. frqdem = period/tofloat(ntim)
>         The demodulation frequency should equal to 1/period.
>         2. frqcut = 0.010
>         The cut-off frequency should equal to 0.5*frqdem.
>         3. nwt    = 41
>         The widow size should equal to 2T-1=21.
>         Moreover, the updated one can not correctly run in NCL 6.4.0.
>         Happy new year.
> Feng
>
>
>
>
> > On Jan 1, 2018, at 1:51 AM, Dennis Shea <shea at ucar.edu> wrote:
> >
> > FYI:
> >
> > I have updated the demodulation pages. I have added 'Complex
> Demodulation' descriptive text section.
> >
> > https://www.ncl.ucar.edu/Applications/spec.shtml
> >
> > The Scripts and Images have been updated to use a new function
> (unwrap_phase; 6.5.0) that 'unwraps' the 'wrapped' phases returned by
> atan2. This replicates the results of Matlab's 'unwrap' function. FYI: I
> have attached the function which will available with the 6.5.0 release.
> >
> > https://www.ncl.ucar.edu/Document/Functions/
> Contributed/unwrap_phase.shtml
> >
> > Happy New Year
> > D
> >
> >
> >
> >
> >
> > On Wed, Dec 27, 2017 at 6:44 AM, Dennis Shea <shea at ucar.edu> wrote:
> >
> >
> > Sent from my iPhone
> >
> > Begin forwarded message:
> >
> >> From: Nai Suan <shif422 at gmail.com>
> >> Date: December 27, 2017 at 8:40:34 AM EST
> >> To: Dennis Shea <shea at ucar.edu>
> >> Subject: Re: [ncl-talk] The question about the complex demodulation
> >>
> >> Dear Dennis Shea,
> >>      Thank you for your information.
> >>      I guess that I found the real answer.
> >>      1) you are right. The calculation of the ‘frqdem’ is wrong in NCL
> example. The right formula is freqdem = 1/11 ; (1/central_period)
> >>      2) The real reason of the difference between NCL and
> Bloomfield(2000) is the use of the function ‘unwrap’. In Bloomfield (2000),
> the phase angles are corrected to produce the smoother phase. I used the
> ‘unwrap’ function in matlab (https://www.mathworks.com/
> help/matlab/ref/unwrap.html).
> >>      Thank you again for your help.
> >>      All the best,
> >> Feng
> >>
> >>> On Dec 21, 2017, at 5:12 AM, Dennis Shea <shea at ucar.edu> wrote:
> >>>
> >>> Hello,
> >>>
> >>> I am going on vacation over the Christmas holidays.
> >>>
> >>> The following is offline from ncl-talk.
> >>> ------------------------------------------------------------
> -----------------
> >>> Maybe you can do some investigation on your own?
> >>>
> >>> Any tool can be used to performdemodulation.
> >>> Other tools offer demodulation as a function or via examples:
> >>>
> >>> Matlab: http://www.mathworks.com/help/signal/ref/demod.html
> >>> R:          https://gist.github.com/mike-lawrence/8565774
> >>> Python: https://currents.soest.hawaii.edu/ocn760_4/_static/complex_
> demod.html
> >>> IDL:       https://github.com/callumenator/idl/blob/master/
> external/Harris/demod.pro
> >>>
> >>> The R example is similar to NCL except it uses a Butterworth low-pass
> filter rather than a Lanczos filter.
> >>>
> >>> Unfortunately, most 'demodulation' are described using electrical
> engineering/communications terminology.
> >>> ========
> >>>
> >>> The steps for complex demodulation about a specific demodulation
> frequency ('frqdem' , 'omega=2*pi*frqdem) follow.
> >>> I have attached all the fortran codes from the book:  demoddriver
> actually calls the appropriate subroutines.
> >>>
> >>> ----
> >>> Bloomfield's program does detrend the sunspot series.
> >>> In NCL, this is left up to the user.   I doubt very much that this has
> any significant effect.
> >>>
> >>> The main computational parts are:
> >>>
> >>> (a) Calculate the real & imaginary compoments.
> >>> Technically, a time series is multiplied by a complex exponential at
> the demodulation frequency. Bloomfield's fortran code:
> >>>
> >>>      do i=1,n
> >>>          arg   = (i-1)*omega
> >>>          d1(i) =  x(i)*cos(arg)*2
> >>>          d2(i) = -x(i)*sin(arg)*2
> >>>       end do
> >>>
> >>> (b) A low-pass filter is applied independently to the real and
> imaginary components from (a).
> >>> Bloomfield uses "least squares lopass filter using convergence factors.
> >>>
> >>>           subroutine lopass(x,n,cutoff,ns)
> >>>
> >>> (c) Calculate the amplitudes and phases using the smoothed values from
> (b). An excerpt from Bloomfield's fortran code:
> >>>
> >>>       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
> >>> ------------------------------------------------------------
> ------------------------------------------------------------
> ------------------
> >>> Parts (a) and (c) are straightforward. NCL uses array syntax rather
> than 'do' loops. SPecifically:
> >>>
> >>>    arg = ispan(0,ny-1,1)*frqdem*2*pi
> >>>    ARG = conform(y, arg, ndim)     ; all dimensions
> >>>                                    ; complex demodulate
> >>>    yr  =  y*cos(ARG)*2             ; 'real' series
> >>>    yi  = -y*sin(ARG)*2             ; 'imaginary'
> >>> ---
> >>>    lanczos
> >>>
> >>> ;---Calculate the amplitudes and phases  (really, no need for where
> function
> >>>    amp = sqrt(yr^2 + yi^2)
> >>>    pha = where(amp.eq.0,  0, atan2(yi, yr))  ; pha = atan2(yi,yr)
> >>>
> >>>  I am sure that lanczos. Also, in my example, I set frqdem=1.0/22  ...
> It should be 1.0/11
> >>>
> >>> Originally, I tried to use 'subroutine lopass'. However, Bloomfield
> read several constants from a file. The values of these constants are not
> available (near as I can tell). The descriptions are:
> >>>
> >>> c     np2 - a power of 2
> >>> c     nom - (np2/2*pi)*omega
> >>> c     npa  - (np2/2*pi) times the pass frequency of the filter
> >>> c     nst   - (np2/2*pi) times the stopfrequency of the filter
> >>>
> >>> These are used to calculate assorted quantities earlier in the program:
> >>>
> >>>   omega = 2*pi*float(nom)/float(np2)
> >>>   cutoff   = pi*float(npa+nst)/float(np2)
> >>>   ns        = np2/(nst-npa)
> >>>
> >>> I had to guess at 'np2'.
> >>>
> >>> Anyway, the calculated values returned by 'lopass' were just not
> realistic.
> >>> As a result, I decided to use NCL's lanczos filter.
> >>>
> >>> https://www.ncl.ucar.edu/Document/Functions/Built-in/
> filwgts_lanczos.shtml
> >>>
> >>> I am sure the filter works fine. What should 'fca' and 'fcb' be set
> to? I used ihp=0 (lowpass),
> >>>
> >>> I have attached 2 tests.
> >>>
> >>> [1]
> >>>
> >>> %> ncl aTEST.ncl
> >>>
> >>> Basically, the same code as on the website. NCL automatically loads
> all the pertinent libraries. Still, I have explicitly included the 6.4.0
> 'demod_cmplx' function.
> >>>       ihp    = 0   ; low pass
> >>>
> >>>       fcb    = totype(-999, typeof(pi))
> >>>       wgt    = filwgts_lanczos (nwt, ihp, frqcut, fcb, nsigma)     ;
> nwt=41
> >>>      ;wgt    = filwgts_lanczos (nwt, ihp, frqdem, fcb, nsigma)
> >>>
> >>>
> >>> [2] bTest.ncl
> >>>
> >>> This uses Bloomfield's fortran code via a shared object. See:
> https://www.ncl.ucar.edu/Document/Tools/WRAPIT.shtml
> >>>
> >>> &>
> >>>
> >>> WRAPIT demod_ncl.f
> >>>
> >>> %> bTEST.ncl
> >>>
> >>>
> >>>
> >>>
> >>>
> >>>
> >>>
> >>> On Wed, Dec 20, 2017 at 12:46 AM, Nai Suan <shif422 at gmail.com> wrote:
> >>> Dear Dennis Shea,
> >>>     Sorry to disturb you again.
> >>>     Do you have any idea to explain the difference of phases in
> Bloomfield’s book and NCL example?
> >>>     Merry Christmas!
> >>>     All the best,
> >>> Frank
> >>>
> >>>
> >>>> On Dec 17, 2017, at 10:47 AM, Nai Suan <shif422 at gmail.com> wrote:
> >>>>
> >>>> 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/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
> >>>>> >
> >>>>>
> >>>>>
> >>>>> <demod_cmplx.ncl><demod_bloomfield_1976.f><DEMODU_nist.f>
> >>>>
> >>>
> >>>
> >>> <sunspot_wolf.1700_1960.txt><demod_ncl.f><aTest.ncl><bTest.
> ncl><demod_cmplx_640.ncl>
> >>
> >
> > <unwrap_phase.ncl><unwrap_phase_matlab.ncl>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180101/3cb2eed1/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: unwrap_phase.ncl_640
Type: application/octet-stream
Size: 1282 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180101/3cb2eed1/attachment.obj>


More information about the ncl-talk mailing list