[ncl-talk] The question about the complex demodulation
Dennis Shea
shea at ucar.edu
Sun Dec 31 17:51:04 MST 2017
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
> <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_200
>> 7/pdfHMANlJLqCE.pdf
>> Kessler: https://faculty.washington.edu/kessler/generals/complex-demo
>> dulation.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>
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20171231/acd3a523/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: unwrap_phase.ncl
Type: application/octet-stream
Size: 1251 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20171231/acd3a523/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: unwrap_phase_matlab.ncl
Type: application/octet-stream
Size: 2122 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20171231/acd3a523/attachment-0001.obj>
More information about the ncl-talk
mailing list