[ncl-talk] The question about the complex demodulation
Dennis Shea
shea at ucar.edu
Wed Dec 20 21:14:40 MST 2017
One thing ... I set the demodulation frequency at 1/22 years, it should be
1/11 years.
In a separate offline communication, I sent you some information and test
scripts.
Cheers
D
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-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>
>
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20171220/eea87392/attachment.html>
More information about the ncl-talk
mailing list