[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