<div dir="ltr"><div><div><div><div>This is a rather long and, perhaps, confusing response. It was written episodically.<br><br>---<br>re: <b>"same mistake" </b>   .... I am not sure to what you are referring. The use of a minus sign?<br><br></div><div>The minus sign does not affect the amplitudes because the quantities are squared:  amp = sqrt(x*x+y*y)<br></div><div>It does affect the phase. Also, use of <b>atan2 </b>or <b>atan </b>to determine the phase are other issues to consider.<br><b><br></b><a href="https://www.ncl.ucar.edu/Document/Functions/Built-in/atan2.shtml">https://www.ncl.ucar.edu/Document/Functions/Built-in/atan2.shtml</a><b><br></b><a href="https://www.ncl.ucar.edu/Document/Functions/Built-in/atan.shtml">https://www.ncl.ucar.edu/Document/Functions/Built-in/atan.shtml</a><b><b><br></b></b></div><div><br></div>PMEL:   <a href="https://www.pmel.noaa.gov/maillists/tmap/ferret_users/fu_2007/pdfHMANlJLqCE.pdf" target="_blank">https://www.pmel.noaa.gov/<wbr>maillists/tmap/ferret_users/<wbr>fu_2007/pdfHMANlJLqCE.pdf</a><br></div>Kessler: <a href="https://faculty.washington.edu/kessler/generals/complex-demodulation.pdf" target="_blank">https://faculty.washington.<wbr>edu/kessler/generals/complex-<wbr>demodulation.pdf</a><br><br></div><div>PMEL references the same source as NCL.<br>Reference: Bloomfield, P. (1976). Fourier Decomposition of Time Series, An Introduction. Wiley, New York, 258pp.</div><div><br>===<br></div>I have attached the (slightly updated) fortran subroutines from Bloomfield  but will include here also::<br>     <b>demod_bloom1976 </b>(page 148)  <br>     <b>polar_bloom1976</b> (page 150)<br><br> ---------------------------<br>      subroutine <b>demod_bloom1976</b>(x,n,omega,d1,<wbr>d2)<br>      implicit none<br>c<br>c Fourier Analysis of Time Series: An Introduction<br>c P. Bloomfield (1976: Wiley); <b>p148</b><br>c Minor code changes to make code more 'modern'<br>c<br>      integer n                        ! INPUT<br>      real    x(n), omega, d1(n), d2(n)<br>      integer i                        ! LOCAL<br>      real    arg<br><br>      do i=1,n<br>         arg   = (i-1)*omega<br><b>         d1(i) =  x(i)*cos(arg)*2<br>         d2(i) = -x(i)*sin(arg)*2</b><br>      end do<br><br>      return<br>      end<br>c ---------------------------<br>      subroutine <b>polar_bloom1976</b>(x,y,n)<br>      implicit none<br>c<br>c Fourier Analysis of Time Series: An Introduction<br>c P. Bloomfield (1976: Wiley); <b>p150</b><br>c Minor code changes to make code more 'modern'<br>c<br>      integer n                        ! INPUT<br>      real    x(n), y(n)<br>      integer i                        ! LOCAL<br>      real    amp, phi<br><br>      do i=1,n<br><b>         amp = sqrt(x(i)*x(i)+y(i)*y(i))<br>         phi = atan2(y(i),x(i))        </b>! radians<br>         x(i)= amp<br>         y(i)= phi<br>      end do<br><br>      return<br>      end<br>c ---------------------------<br>c Bloomfield CALL SEQUENCE <br>c<br>c omega  = (2*pi*nom)/np2<br>c cutoff = (pi*(npa+nst))/np2<br>c ns     = np2/(nst-npa)<br>c call demod_bloom1976 (x,n,omega,d1,d2) ! DEMODULATION<br>c call lopass_bloom1976(d1,cutoff,ns)<wbr>    ! SMOOTH COSINE TERM<br>c call lopass_bloom1976(d2,cutoff,ns)<wbr>    ! SMOOTH SINE TERM<br>c lim = n-2*ns<br>c call polar_bloom1976(d1,d2,lim)    <wbr>    ! AMPLITUDE & PHASES<br>==============================<wbr>=<br></div><div><b>Bloomfield's book uses a different low-pass filter than NCL (Lanczos) so results are not directly ('exactly') comparable.</b><br>==============================<wbr>=<br></div><div>NCL's  (attached) <b>demod_cmplx</b> uses array notation (no 'do' loops):<br><br>...<br>   yr  =  y*cos(ARG)*2             ; 'real' series<br>   yi  = <b>-</b>y*sin(ARG)*2              ; 'imaginary'   <br>...<br>   wgt    = filwgts_lanczos (nwt, ihp, frqcut, fcb, nsigma)  <br><br>;---Apply filter weights to raw demodulated series (nwt/2 point 'lost; at start/end)<br><br>   yr = wgt_runave_n_Wrap ( yr, wgt, 0, 0) <br>   yi = wgt_runave_n_Wrap ( yi, wgt, 0, 0) <br>  <br>;---Use the low pass filtered values; <b>atan2 </b>uses y/x => yi/yr<br>   <br>   amp  = sqrt(yr^2 + yi^2)<br>   pha  = where(amp.eq.zero, zero, <b>atan2</b>(yi, yr))    ; primary phase (radians)<br><br>   pha2 = conform(pha, zero, -1)                     ; secondary phase (Bloomfield)<br>   pha2 = where(pha.ge.zero, pha-pi2, pha+pi2)   <br><br>==============================<wbr>====<br></div><div>You could take NCL's attached '<b>demod_cmplx</b>'; copy it to create your own function (say): <b>demod_cmplx_nai<br><br></b></div><div>and change what you want .... for example.... remove the minus sign<br><br>   yr  =  y*cos(ARG)*2             ; 'real' series<br>   yi  =   y*sin(ARG)*2              ; 'imaginary'   </div><div><b></b></div><div><b><br></b></div><div>Play with it.<br><br>==============================<wbr>====<br></div><div>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.<br><br>==============================<wbr>=====<br></div><div>Archived National Institute of Standards fortran software.   <br></div><div><br></div><div><b>DEMODU.f</b>  (attached): Does the same thing as Bloomfield/NCL. <br>The use af AMPL and PHAS names is because this was part of a larger library that reused array space.<br><br></div><div>The following <br></div><div><br><a href="http://www.itl.nist.gov/div898/software/datapac/DEMOD.f" target="_blank">http://www.itl.nist.gov/<wbr>div898/software/datapac/DEMOD.<wbr>f</a><br><br></div><div>uses <br><br>      Y1(I)=X(I)*COS(6.2831853*F*AI)<br>      Y2(I)=X(I)*SIN(6.2831853*F*AI)<br><br></div><div>but then for phases uses <b>atan </b>(not <b>atan2</b>)<b><br><br></b>      Z(I)=<b>ATAN</b>(Y1(I)/Y2(I))<b>     </b>!  atan(yr/yi)<b>  ... ?</b>implicitly handles - sign?<br><b><br></b></div><div>I am done.<br><br></div><div>Cheers <br></div><div>D<b><br></b></div><div><br></div><div><br><br><br><div><br><div><div><br><br></div></div></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Dec 13, 2017 at 8:31 AM, Nai Suan <span dir="ltr"><<a href="mailto:shif422@gmail.com" target="_blank">shif422@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Dear Dennis Shea,<br>
        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.<br>
        All the best,<br>
Feng<br>
<div class="HOEnZb"><div class="h5"><br>
<br>
<br>
> On Dec 13, 2017, at 10:02 AM, Nai Suan <<a href="mailto:shif422@gmail.com">shif422@gmail.com</a>> wrote:<br>
><br>
> Dear Dennis Shea,<br>
>       Thank you for your reply.<br>
>       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.<br>
>       You mentioned NCL example was translated from Bloomfield’s subroutine. Is there any possible reason that the translation has some problems?<br>
>       Anyway, if you are interested, I can send Bloomfield’s 2000 book to you.<br>
>       All the best,<br>
> Frank<br>
><br>
><br>
>> On Dec 13, 2017, at 12:15 AM, Dennis Shea <<a href="mailto:shea@ucar.edu">shea@ucar.edu</a>> wrote:<br>
>><br>
>> If the low-pass filtering (Lanczos) and smoothing (weighted running average) don't explain the differences... I do not know.<br>
>> 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.<br>
>> ----<br>
>> <a href="http://www.ncl.ucar.edu/Applications/spec.shtml" rel="noreferrer" target="_blank">http://www.ncl.ucar.edu/<wbr>Applications/spec.shtml</a><br>
>><br>
>> Data Analysis: Spectral analysis, Complex Demodulation<br>
>><br>
>> The 1st complex demodulation example uses the 1700-1960 sunspots<br>
>> <a href="http://www.ncl.ucar.edu/Applications/Scripts/demod_cmplx_1.ncl" rel="noreferrer" target="_blank">http://www.ncl.ucar.edu/<wbr>Applications/Scripts/demod_<wbr>cmplx_1.ncl</a><br>
>> ----<br>
>><br>
>> The 'demod_cmplx' function can be seen at:<br>
>><br>
>> %> less  $NCARG_ROOT/lib/ncarg/<wbr>nclscripts/csm/contributed.ncl<br>
>><br>
>> Basically, it is a translation (fortran-77  ===> NCL) of Bloomfield's 1976 f77 subroutine.<br>
>><br>
>><br>
>><br>
>><br>
>><br>
>><br>
>> On Tue, Dec 12, 2017 at 4:03 AM, Nai Suan <<a href="mailto:shif422@gmail.com">shif422@gmail.com</a>> wrote:<br>
>> Dear colleague,<br>
>>         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?<br>
>>         Hope to hear your answer.<br>
>>         All the best,<br>
>> Frank<br>
>><br>
>><br>
>><br>
>> ><br>
>> ><br>
>> ><br>
>> ><br>
>> ><br>
>><br>
>><br>
>> ______________________________<wbr>_________________<br>
>> ncl-talk mailing list<br>
>> <a href="mailto:ncl-talk@ucar.edu">ncl-talk@ucar.edu</a><br>
>> List instructions, subscriber options, unsubscribe:<br>
>> <a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" rel="noreferrer" target="_blank">http://mailman.ucar.edu/<wbr>mailman/listinfo/ncl-talk</a><br>
><br>
<br>
</div></div></blockquote></div><br></div>