<div dir="ltr"><div><div><div><div><div><div><div>[1] re: <br>"In my experience, the cut off frequency of the low-pass filter (frqcut =
0.5*frqdem) will get a wrong amplitude and a wrong phase"<br><br></div><div>Compare with <b>Bloomfield (2000), figure 7.11 which uses a frqcut=0.04968</b><br><br><a href="https://iujfk.files.wordpress.com/2013/04/bloomfield-2000-fourier-analysis-of-time-series-an-introduction-2ed.pdf" target="_blank">https://iujfk.files.wordpress.<wbr>com/2013/04/bloomfield-2000-<wbr>fourier-analysis-of-time-<wbr>series-an-introduction-2ed.pdf</a><br>------<br>with the results of NCL's left figure (Amplitude and Unwrapped-Phase) which are derived via: <b> frqcut = 0.5*frqdem</b>=0.5/11=1/22=<b>0.04545</b> <br><br><a href="http://www.ncl.ucar.edu/Applications/Images/demod_cmplx_1_lg.png" target="_blank">http://www.ncl.ucar.edu/<wbr>Applications/Images/demod_<wbr>cmplx_1_lg.png</a><br><br></div><div><br></div>[a] amplitudes visually match. Bloomfield and NCL use different filters and <i>slightly</i> different low-pass cutoff frequencies so 'exact' match is unlikely. <br></div>[b] unwrapped phases share <b>almost the identical structure</b>. The magnitudes are very different. Why? I do no know. Perhaps, Bloomfield/S+ uses a different algirithm. As noted NCL uses <b>unwrap_phase.<br><br><a href="http://www.ncl.ucar.edu/Document/Functions/Contributed/unwrap_phase.shtml" target="_blank">http://www.ncl.ucar.edu/<wbr>Document/Functions/<wbr>Contributed/unwrap_phase.shtml</a><br></b></div>This matches Matlab's unwrap. <br><br>======<br><br></div>[2] re: "<span class="gmail-m_-4434870793298626144m_-6556111732610160925Apple-tab-span" style="white-space:pre-wrap"></span>Why do you think that frqcut and frqdem are not related?<br><br></div>(a) frqdem: Often, this is suggested by some physical principal or theory. This is used to create the time series of <i><u>unsmoothed</u> </i>real and imaginary components.<br></div><div><br></div>(b) frqcut: <i><u>User</u></i> must decide on what low-pass variability is desirable for research objective(s).<br><br></div>This is why I included:<br><p>
Further, as suggested by the
<a href="http://www.uni-muenster.de/ZIV.BennoSueselbeck/s-html/helpfiles/demod.html"><b>S-Plus demod</b></a> documentation:
</p><pre> <b>To better understand the results of complex demodulation
several lowpass filters should be tried:</b> the smaller the pass band,
the less instantaneous in time but more specific in frequency is the result</pre><div><div><div><div>======<br></div><div>[3] re: <span>(i) "What do you think the dominant frequency of original signal is not equal to the demodulation frequency (frqdem)?</span><span><span> (ii) </span><span>In my opinion, the precondition for complex demodulation is that
the dominate frequency of original signal must be equal to
the demodulation frequency,"</span></span><span>"</span><span>In my opinion, the precondition for complex demodulation is that
the dominate frequency of original signal must be equal to
the demodulation frequency,"<br><br></span></div><div><span>Well, in principle, any <b>frqdem</b> (0<frqdem<0.5) could be specified. However, if the specified frequency is not physically important (your terminology: "dominant"), why would you want to study it?<br></span></div><div><span><br>======<br></span></div><div class="gmail_extra"><span>[4] re: </span><span>t<b>he nwt is suggested as nwt= 2*T-1</b>. <br><br>That just refers to the ?triangle? filter they choose to use.<br><br>=====<br>[5] "</span><span><span>If the estimations of the frqcut and the nwt can be arbitrarily
changed, how to know whether the demodulated amplitude and phase are
right with the suitable values of freqcut and nwt.</span><b>" That is your job as a scientist to decide.<br></b><br>+<br>++++++++++<br></span></div><div class="gmail_extra"><span>As I noted on the 'demodulation page':<br></span><br><span>In fact, complex demodulation (like spectral analysis), could be viewed to be as much 'art' as 'science'.
Different filters, <em>frqcut</em> and <em>nwt must be chosen.<br><br>===========<br></em></span><br><span><em><span>+++++++++<b><br></b></span></em></span><div class="gmail_extra"><span>I think it might be good for you to talk with a statistician ... which I am not!<br><br></span></div><div class="gmail_extra"><span>Good Luck<br></span></div><div class="gmail_extra"><span>D<br></span></div></div><div class="gmail_extra"><span><b></b></span></div><div class="gmail_extra"><span><b><br><br></b></span><div class="gmail_quote">On Mon, Jan 1, 2018 at 1:43 PM, 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:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div style="overflow-wrap: break-word;">Dear Dennis Shea,<div><span class="gmail-m_-4434870793298626144m_-6556111732610160925Apple-tab-span" style="white-space:pre-wrap"> </span>Thank you for updating it according to my comments. Until now, the results of NCL is coincident with the example in Bloomfield (2000)’s book.</div><div><span class="gmail-m_-4434870793298626144m_-6556111732610160925Apple-tab-span" style="white-space:pre-wrap"> </span>I still have a question about the parameter estimations of the frqcut and nwt.</div><div><span class="gmail-m_-4434870793298626144m_-6556111732610160925Apple-tab-span" style="white-space:pre-wrap"> </span>Why do you think that frqcut and frqdem are not related? As I understand, the frqcut must be less than the frqdem. In my experience, the cut off frequency of the low-pass filter (frqcut = 0.5*frqdem) will get a wrong amplitude and a wrong phase, and a smaller cut-off frequency is necessary (frqcut = 0.25*frqdem), when the original signal has very strong other cycles, which are different with the demodulation frequency.</div><span><span class="gmail-m_-4434870793298626144m_-6556111732610160925Apple-tab-span" style="white-space:pre-wrap"> </span>What do you think the dominant frequency of original signal is not equal to the demodulation frequency (frqdem)? Whether is this condition still suitable for the complex demodulation? In my opinion, the precondition for complex demodulation is that the dominate frequency of original signal must be equal to the demodulation frequency, because the precondition for complex demodulation supposes that is the original data contains a perturbed periodic component, whose amplitude and phase are slowly changing ones. For example, it is not suitable for the complex demodulation of sunspot number using the demodulation frequency (1/11), if the sunspot number has very stronger 30-years cycle than 11-year cycle.<br></span><span><span class="gmail-m_-4434870793298626144m_-6556111732610160925Apple-tab-span" style="white-space:pre-wrap"> </span>Moreover, the nwt is suggested as nwt= 2*T-1 in two files 'complex-demodulation.pdf</span><span>’ and ‘pdfHMANlJLqCE.pdf’. However, as we known the nwt depends on the type of the filter that you select. If the estimations of the frqcut and the nwt can be arbitrarily changed, how to know whether the demodulated amplitude and phase are right with the suitable values of freqcut and nwt.</span><div><span><span class="gmail-m_-4434870793298626144m_-6556111732610160925Apple-tab-span" style="white-space:pre-wrap"> </span>All the best,<br></span><div>Feng</div><div><div><div><br><blockquote type="cite"><div><div class="gmail-m_-4434870793298626144h5"><div>On Jan 1, 2018, at 8:41 PM, Dennis Shea <<a href="mailto:shea@ucar.edu" target="_blank">shea@ucar.edu</a>> wrote:</div><br class="gmail-m_-4434870793298626144m_-6556111732610160925Apple-interchange-newline"></div></div><div><div><div class="gmail-m_-4434870793298626144h5"><div dir="ltr"><div><div><div>[1] frqdem = period/tofloat(ntim)<br>
The demodulation frequency should equal to 1/period.<br><br></div>Yes ... I should have coded what I wrote. :-(<br><br>===<br>[2] frqcut = 0.010<br>
The cut-off frequency should equal to 0.5*frqdem<br><br></div>No. There is relation between frqdem and frqcut. <br></div>The latter is strictly related to the desired low-pass filter characteristics.<br><div>===<br>[3] nwt = 41<br><br></div><div>This was used only because Bloomfield (1976) used that number of weights for his least squares low pass filter.<br></div><div>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.<br><br>++++++++++++++++<br></div><div>I updated the spectral analysis and complex-demodulation page:<br><br><a href="https://www.ncl.ucar.edu/Applications/spec.shtml" target="_blank">https://www.ncl.ucar.edu/Appli<wbr>cations/spec.shtml</a><br><br></div><div>The sunspot demodulation example shows two figures with different cutoff frequencies. All other parameters are the same.<br></div><div>+++++++++++++++++<br><br></div><div>re: Moreover, the updated one can not correctly run in NCL 6.4.0.<br><br></div><div>The 'unwrap_phase' I attached used a new feature of NCL: a 'true' <b>elseif<br><br> <a href="https://www.ncl.ucar.edu/future_release.shtml" target="_blank">https://www.ncl.ucar.edu/futur<wbr>e_release.shtml</a><br></b></div><div><br></div><div>===<br></div><div>The unwrap_phase_matlab.ncl could be used.<br><br></div><div>I have attached a version of unwrap_phase that should run with 6.4.0<br><br></div><div>Good Luck<br></div><div>D<br></div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Mon, Jan 1, 2018 at 12:43 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:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Dear Dennis Shea,<br>
Thank you for your updating.<br>
Some assignments are still obscure in the updated one, which cause that the results are still different with Bloomfield’s book.<br>
1. frqdem = period/tofloat(ntim)<br>
The demodulation frequency should equal to 1/period.<br>
2. frqcut = 0.010<br>
The cut-off frequency should equal to 0.5*frqdem.<br>
3. nwt = 41<br>
The widow size should equal to 2T-1=21.<br>
Moreover, the updated one can not correctly run in NCL 6.4.0.<br>
Happy new year.<br>
Feng<br>
<div><div class="gmail-m_-4434870793298626144m_-6556111732610160925h5"><br>
<br>
<br>
<br>
> On Jan 1, 2018, at 1:51 AM, Dennis Shea <<a href="mailto:shea@ucar.edu" target="_blank">shea@ucar.edu</a>> wrote:<br>
><br>
> FYI:<br>
><br>
> I have updated the demodulation pages. I have added 'Complex Demodulation' descriptive text section.<br>
><br>
> <a href="https://www.ncl.ucar.edu/Applications/spec.shtml" rel="noreferrer" target="_blank">https://www.ncl.ucar.edu/Appli<wbr>cations/spec.shtml</a><br>
><br>
> 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.<br>
><br>
> <a href="https://www.ncl.ucar.edu/Document/Functions/Contributed/unwrap_phase.shtml" rel="noreferrer" target="_blank">https://www.ncl.ucar.edu/Docum<wbr>ent/Functions/Contributed/unwr<wbr>ap_phase.shtml</a><br>
><br>
> Happy New Year<br>
> D<br>
><br>
><br>
><br>
><br>
><br>
> On Wed, Dec 27, 2017 at 6:44 AM, Dennis Shea <<a href="mailto:shea@ucar.edu" target="_blank">shea@ucar.edu</a>> wrote:<br>
><br>
><br>
> Sent from my iPhone<br>
><br>
> Begin forwarded message:<br>
><br>
>> From: Nai Suan <<a href="mailto:shif422@gmail.com" target="_blank">shif422@gmail.com</a>><br>
>> Date: December 27, 2017 at 8:40:34 AM EST<br>
>> To: Dennis Shea <<a href="mailto:shea@ucar.edu" target="_blank">shea@ucar.edu</a>><br>
>> Subject: Re: [ncl-talk] The question about the complex demodulation<br>
>><br>
>> Dear Dennis Shea,<br>
>> Thank you for your information.<br>
>> I guess that I found the real answer.<br>
>> 1) you are right. The calculation of the ‘frqdem’ is wrong in NCL example. The right formula is freqdem = 1/11 ; (1/central_period)<br>
>> 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 (<a href="https://www.mathworks.com/help/matlab/ref/unwrap.html" rel="noreferrer" target="_blank">https://www.mathworks.com/hel<wbr>p/matlab/ref/unwrap.html</a>).<br>
>> Thank you again for your help.<br>
>> All the best,<br>
>> Feng<br>
>><br>
>>> On Dec 21, 2017, at 5:12 AM, Dennis Shea <<a href="mailto:shea@ucar.edu" target="_blank">shea@ucar.edu</a>> wrote:<br>
>>><br>
>>> Hello,<br>
>>><br>
>>> I am going on vacation over the Christmas holidays.<br>
>>><br>
>>> The following is offline from ncl-talk.<br>
>>> ------------------------------<wbr>------------------------------<wbr>-----------------<br>
>>> Maybe you can do some investigation on your own?<br>
>>><br>
>>> Any tool can be used to performdemodulation.<br>
>>> Other tools offer demodulation as a function or via examples:<br>
>>><br>
>>> Matlab: <a href="http://www.mathworks.com/help/signal/ref/demod.html" rel="noreferrer" target="_blank">http://www.mathworks.com/help/<wbr>signal/ref/demod.html</a><br>
>>> R: <a href="https://gist.github.com/mike-lawrence/8565774" rel="noreferrer" target="_blank">https://gist.github.com/mike-l<wbr>awrence/8565774</a><br>
>>> Python: <a href="https://currents.soest.hawaii.edu/ocn760_4/_static/complex_demod.html" rel="noreferrer" target="_blank">https://currents.soest.hawaii.<wbr>edu/ocn760_4/_static/complex_d<wbr>emod.html</a><br>
>>> IDL: <a href="https://github.com/callumenator/idl/blob/master/external/Harris/demod.pro" rel="noreferrer" target="_blank">https://github.com/callumenat<wbr>or/idl/blob/master/external/Ha<wbr>rris/demod.pro</a><br>
>>><br>
>>> The R example is similar to NCL except it uses a Butterworth low-pass filter rather than a Lanczos filter.<br>
>>><br>
>>> Unfortunately, most 'demodulation' are described using electrical engineering/communications terminology.<br>
>>> ========<br>
>>><br>
>>> The steps for complex demodulation about a specific demodulation frequency ('frqdem' , 'omega=2*pi*frqdem) follow.<br>
>>> I have attached all the fortran codes from the book: demoddriver actually calls the appropriate subroutines.<br>
>>><br>
>>> ----<br>
>>> Bloomfield's program does detrend the sunspot series.<br>
>>> In NCL, this is left up to the user. I doubt very much that this has any significant effect.<br>
>>><br>
>>> The main computational parts are:<br>
>>><br>
>>> (a) Calculate the real & imaginary compoments.<br>
>>> Technically, a time series is multiplied by a complex exponential at the demodulation frequency. Bloomfield's fortran code:<br>
>>><br>
>>> do i=1,n<br>
>>> arg = (i-1)*omega<br>
>>> d1(i) = x(i)*cos(arg)*2<br>
>>> d2(i) = -x(i)*sin(arg)*2<br>
>>> end do<br>
>>><br>
>>> (b) A low-pass filter is applied independently to the real and imaginary components from (a).<br>
>>> Bloomfield uses "least squares lopass filter using convergence factors.<br>
>>><br>
>>> subroutine lopass(x,n,cutoff,ns)<br>
>>><br>
>>> (c) Calculate the amplitudes and phases using the smoothed values from (b). An excerpt from Bloomfield's fortran code:<br>
>>><br>
>>> do i=1,n<br>
>>> amp = sqrt(x(i)*x(i)+y(i)*y(i))<br>
>>> phi = atan2(y(i),x(i)) ! radians<br>
>>> x(i)= amp<br>
>>> y(i)= phi<br>
>>> end do<br>
>>> ------------------------------<wbr>------------------------------<wbr>------------------------------<wbr>------------------------------<wbr>------------------<br>
>>> Parts (a) and (c) are straightforward. NCL uses array syntax rather than 'do' loops. SPecifically:<br>
>>><br>
>>> arg = ispan(0,ny-1,1)*frqdem*2*pi<br>
>>> ARG = conform(y, arg, ndim) ; all dimensions<br>
>>> ; complex demodulate<br>
>>> yr = y*cos(ARG)*2 ; 'real' series<br>
>>> yi = -y*sin(ARG)*2 ; 'imaginary'<br>
>>> ---<br>
>>> lanczos<br>
>>><br>
>>> ;---Calculate the amplitudes and phases (really, no need for where function<br>
>>> amp = sqrt(yr^2 + yi^2)<br>
>>> pha = where(amp.eq.0, 0, atan2(yi, yr)) ; pha = atan2(yi,yr)<br>
>>><br>
>>> I am sure that lanczos. Also, in my example, I set frqdem=1.0/22 ... It should be 1.0/11<br>
>>><br>
>>> 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:<br>
>>><br>
>>> c np2 - a power of 2<br>
>>> c nom - (np2/2*pi)*omega<br>
>>> c npa - (np2/2*pi) times the pass frequency of the filter<br>
>>> c nst - (np2/2*pi) times the stopfrequency of the filter<br>
>>><br>
>>> These are used to calculate assorted quantities earlier in the program:<br>
>>><br>
>>> omega = 2*pi*float(nom)/float(np2)<br>
>>> cutoff = pi*float(npa+nst)/float(np2)<br>
>>> ns = np2/(nst-npa)<br>
>>><br>
>>> I had to guess at 'np2'.<br>
>>><br>
>>> Anyway, the calculated values returned by 'lopass' were just not realistic.<br>
>>> As a result, I decided to use NCL's lanczos filter.<br>
>>><br>
>>> <a href="https://www.ncl.ucar.edu/Document/Functions/Built-in/filwgts_lanczos.shtml" rel="noreferrer" target="_blank">https://www.ncl.ucar.edu/Docum<wbr>ent/Functions/Built-in/filwgts<wbr>_lanczos.shtml</a><br>
>>><br>
>>> I am sure the filter works fine. What should 'fca' and 'fcb' be set to? I used ihp=0 (lowpass),<br>
>>><br>
>>> I have attached 2 tests.<br>
>>><br>
>>> [1]<br>
>>><br>
>>> %> ncl aTEST.ncl<br>
>>><br>
>>> 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.<br>
>>> ihp = 0 ; low pass<br>
>>><br>
>>> fcb = totype(-999, typeof(pi))<br>
>>> wgt = filwgts_lanczos (nwt, ihp, frqcut, fcb, nsigma) ; nwt=41<br>
>>> ;wgt = filwgts_lanczos (nwt, ihp, frqdem, fcb, nsigma)<br>
>>><br>
>>><br>
>>> [2] bTest.ncl<br>
>>><br>
>>> This uses Bloomfield's fortran code via a shared object. See: <a href="https://www.ncl.ucar.edu/Document/Tools/WRAPIT.shtml" rel="noreferrer" target="_blank">https://www.ncl.ucar.edu/Docum<wbr>ent/Tools/WRAPIT.shtml</a><br>
>>><br>
>>> &><br>
>>><br>
>>> WRAPIT demod_ncl.f<br>
>>><br>
>>> %> bTEST.ncl<br>
>>><br>
>>><br>
>>><br>
>>><br>
>>><br>
>>><br>
>>><br>
>>> On Wed, Dec 20, 2017 at 12:46 AM, Nai Suan <<a href="mailto:shif422@gmail.com" target="_blank">shif422@gmail.com</a>> wrote:<br>
>>> Dear Dennis Shea,<br>
>>> Sorry to disturb you again.<br>
>>> Do you have any idea to explain the difference of phases in Bloomfield’s book and NCL example?<br>
>>> Merry Christmas!<br>
>>> All the best,<br>
>>> Frank<br>
>>><br>
>>><br>
>>>> On Dec 17, 2017, at 10:47 AM, Nai Suan <<a href="mailto:shif422@gmail.com" target="_blank">shif422@gmail.com</a>> wrote:<br>
>>>><br>
>>>> Dear Dennis Shea,<br>
>>>> Thank you for your huge effect.<br>
>>>> I agree with you, a minus sign does not affect the amplitude but affect the phase.<br>
>>>><br>
>>>> 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.<br>
>>>> 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.<br>
>>>> 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.<br>
>>>><br>
>>>> 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.<br>
>>>><br>
>>>> 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?<br>
>>>> All the best,<br>
>>>> Frank<br>
>>>><br>
>>>><br>
>>>>> On Dec 15, 2017, at 9:56 PM, Dennis Shea <<a href="mailto:shea@ucar.edu" target="_blank">shea@ucar.edu</a>> wrote:<br>
>>>>><br>
>>>>> This is a rather long and, perhaps, confusing response. It was written episodically.<br>
>>>>><br>
>>>>> ---<br>
>>>>> re: "same mistake" .... I am not sure to what you are referring. The use of a minus sign?<br>
>>>>><br>
>>>>> The minus sign does not affect the amplitudes because the quantities are squared: amp = sqrt(x*x+y*y)<br>
>>>>> It does affect the phase. Also, use of atan2 or atan to determine the phase are other issues to consider.<br>
>>>>><br>
>>>>> <a href="https://www.ncl.ucar.edu/Document/Functions/Built-in/atan2.shtml" rel="noreferrer" target="_blank">https://www.ncl.ucar.edu/Docum<wbr>ent/Functions/Built-in/atan2.s<wbr>html</a><br>
>>>>> <a href="https://www.ncl.ucar.edu/Document/Functions/Built-in/atan.shtml" rel="noreferrer" target="_blank">https://www.ncl.ucar.edu/Docum<wbr>ent/Functions/Built-in/atan.sh<wbr>tml</a><br>
>>>>><br>
>>>>> PMEL: <a href="https://www.pmel.noaa.gov/maillists/tmap/ferret_users/fu_2007/pdfHMANlJLqCE.pdf" rel="noreferrer" target="_blank">https://www.pmel.noaa.gov/mai<wbr>llists/tmap/ferret_users/fu_20<wbr>07/pdfHMANlJLqCE.pdf</a><br>
>>>>> Kessler: <a href="https://faculty.washington.edu/kessler/generals/complex-demodulation.pdf" rel="noreferrer" target="_blank">https://faculty.washington.edu<wbr>/kessler/generals/complex-demo<wbr>dulation.pdf</a><br>
>>>>><br>
>>>>> PMEL references the same source as NCL.<br>
>>>>> Reference: Bloomfield, P. (1976). Fourier Decomposition of Time Series, An Introduction. Wiley, New York, 258pp.<br>
>>>>><br>
>>>>> ===<br>
>>>>> I have attached the (slightly updated) fortran subroutines from Bloomfield but will include here also::<br>
>>>>> demod_bloom1976 (page 148)<br>
>>>>> polar_bloom1976 (page 150)<br>
>>>>><br>
>>>>> ---------------------------<br>
>>>>> subroutine demod_bloom1976(x,n,omega,d1,d<wbr>2)<br>
>>>>> implicit none<br>
>>>>> c<br>
>>>>> c Fourier Analysis of Time Series: An Introduction<br>
>>>>> c P. Bloomfield (1976: Wiley); p148<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>
>>>>> d1(i) = x(i)*cos(arg)*2<br>
>>>>> d2(i) = -x(i)*sin(arg)*2<br>
>>>>> end do<br>
>>>>><br>
>>>>> return<br>
>>>>> end<br>
>>>>> c ---------------------------<br>
>>>>> subroutine polar_bloom1976(x,y,n)<br>
>>>>> implicit none<br>
>>>>> c<br>
>>>>> c Fourier Analysis of Time Series: An Introduction<br>
>>>>> c P. Bloomfield (1976: Wiley); p150<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>
>>>>> amp = sqrt(x(i)*x(i)+y(i)*y(i))<br>
>>>>> phi = atan2(y(i),x(i)) ! 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) ! AMPLITUDE & PHASES<br>
>>>>> ==============================<wbr>=<br>
>>>>> Bloomfield's book uses a different low-pass filter than NCL (Lanczos) so results are not directly ('exactly') comparable.<br>
>>>>> ==============================<wbr>=<br>
>>>>> NCL's (attached) demod_cmplx uses array notation (no 'do' loops):<br>
>>>>><br>
>>>>> ...<br>
>>>>> yr = y*cos(ARG)*2 ; 'real' series<br>
>>>>> yi = -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; atan2 uses y/x => yi/yr<br>
>>>>><br>
>>>>> amp = sqrt(yr^2 + yi^2)<br>
>>>>> pha = where(amp.eq.zero, zero, atan2(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>
>>>>> You could take NCL's attached 'demod_cmplx'; copy it to create your own function (say): demod_cmplx_nai<br>
>>>>><br>
>>>>> 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'<br>
>>>>><br>
>>>>> Play with it.<br>
>>>>><br>
>>>>> ==============================<wbr>====<br>
>>>>> 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>
>>>>> Archived National Institute of Standards fortran software.<br>
>>>>><br>
>>>>> DEMODU.f (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>
>>>>> The following<br>
>>>>><br>
>>>>> <a href="http://www.itl.nist.gov/div898/software/datapac/DEMOD.f" rel="noreferrer" target="_blank">http://www.itl.nist.gov/div898<wbr>/software/datapac/DEMOD.f</a><br>
>>>>><br>
>>>>> uses<br>
>>>>><br>
>>>>> Y1(I)=X(I)*COS(6.2831853*F*AI<wbr>)<br>
>>>>> Y2(I)=X(I)*SIN(6.2831853*F*AI<wbr>)<br>
>>>>><br>
>>>>> but then for phases uses atan (not atan2)<br>
>>>>><br>
>>>>> Z(I)=ATAN(Y1(I)/Y2(I)) ! atan(yr/yi) ... ?implicitly handles - sign?<br>
>>>>><br>
>>>>> I am done.<br>
>>>>><br>
>>>>> Cheers<br>
>>>>> D<br>
>>>>><br>
>>>>><br>
>>>>><br>
>>>>><br>
>>>>><br>
>>>>><br>
>>>>><br>
>>>>><br>
>>>>> On Wed, Dec 13, 2017 at 8:31 AM, Nai Suan <<a href="mailto:shif422@gmail.com" target="_blank">shif422@gmail.com</a>> wrote:<br>
>>>>> 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>
>>>>><br>
>>>>><br>
>>>>><br>
>>>>> > On Dec 13, 2017, at 10:02 AM, Nai Suan <<a href="mailto:shif422@gmail.com" target="_blank">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" target="_blank">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/Applic<wbr>ations/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/Applic<wbr>ations/Scripts/demod_cmplx_1.n<wbr>cl</a><br>
>>>>> >> ----<br>
>>>>> >><br>
>>>>> >> The 'demod_cmplx' function can be seen at:<br>
>>>>> >><br>
>>>>> >> %> less $NCARG_ROOT/lib/ncarg/nclscrip<wbr>ts/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" target="_blank">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" target="_blank">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/mailma<wbr>n/listinfo/ncl-talk</a><br>
>>>>> ><br>
>>>>><br>
>>>>><br>
>>>>> <demod_cmplx.ncl><demod_bloomf<wbr>ield_1976.f><DEMODU_nist.f><br>
>>>><br>
>>><br>
>>><br>
>>> <sunspot_wolf.1700_1960.txt><d<wbr>emod_ncl.f><aTest.ncl><bTest.n<wbr>cl><demod_cmplx_640.ncl><br>
>><br>
><br>
</div></div>> <unwrap_phase.ncl><unwrap_phas<wbr>e_matlab.ncl><br>
<br>
</blockquote></div><br></div>
</div></div><span id="gmail-m_-4434870793298626144m_-6556111732610160925cid:5AE2D1FB-EDA3-4D6A-94F5-31897809433A"><unwrap_phase.ncl_640></span></div></blockquote></div><br></div></div></div></div></blockquote></div><br></div></div></div></div></div>