Guest Post by Willis Eschenbach
While investigating the question of cycles in climate datasets (Part 1, Part 2), I invented a method I called “sinusoidal periodicity”. What I did was to fit a sine wave of various periods to the data, and record the amplitude of the best fit. I figured it had been invented before, so I asked people what I was doing and what its name was. I also asked if there was a faster way to do it, as my method does a lot of optimization (fitting) and thus is slow. An alert reader, Michael Gordon, pointed out that I was doing a type Fast Fourier Transform (FFT) … and provided a link to Nick Stokes’ R code to verify that indeed, my results are identical to the periodogram of his Fast Fourier Transform. So, it turns out that what I’ve invented can be best described as the “Slow Fourier Transform”, since it does exactly what the FFT does, only much slower … which sounds like bad news.
My great thanks to Michael, however, because actually I’m stoked to find out that I’m doing a Fourier transform. First, I greatly enjoy coming up with new ideas on my own and then finding out people have thought of them before me. Some folks might see that as a loss, finding out that someone thought of my invention or innovation before I did. But to me, that just means that my self-education is on the right track, and I’m coming up with valuable stuff. And in this case it also means that my results are a recognized quantity, a periodogram of the data. This is good news because people already understand what it is I’m showing.
Figure 1. Slow Fourier transform periodograms of four long-term surface air temperature datasets. Values are the peak-to-peak amplitude of the best-fit sine wave at each cycle length. The longest period shown in each panel is half the full length of the dataset. Top panel is Armagh Observatory in Ireland. The second panel is the Central England Temperature (CET), which is an average of three stations in central England. Third panel is the Berkeley Earth global temperature dataset. The fourth panel shows the HadCRUT4 global temperature dataset. Note that the units are in degrees C, and represent the peak-to-peak swings in temperature at each given cycle length. Data in color are significant after adjustment for autocorrelation at the 90% level. Significance is calculated after removing the monthly seasonal average variations.
I’m also overjoyed that my method gives identical results to its much speedier cousin, the Fast Fourier transform (FFT), because the Slow Fourier Transform (SFT) has a number of very significant advantages over the FFT. These advantages are particularly important in climate science.
The first big advantage is that the SFT is insensitive to gaps in the data. For example, the Brest tide data goes back to 1807, but there are some missing sections, e.g. from 1836-1846 and 1857-1860. As far as I know, the FFT cannot analyze the full length of the Brest data in one block, but that makes no difference to the SFT. It can utilize all of the data. As you can imagine, in climate science this is a very common issue, so this will allow people to greatly extend the usage of the Fourier transform.
The second big advantage is that the SFT can be used on an irregularly spaced time series. The FFT requires data that is uniformly spaced in time. But there’s a lot of valuable irregularly spaced climate data out there. The slow Fourier transform allows us to calculate the periodogram of the cycles in that irregular data, regardless of the timing of the observations. Even if all you have are observations scattered at various times throughout the year, with entire years missing and some years only having two observations while other years have two hundred observations … no matter. All that affects is the error of the results, it doesn’t prevent the calculation as it does with the FFT.
The third advantage is that the slow Fourier transform is explainable in layman’s terms. If you tell folks that you are transforming data from the time domain to the frequency domain, people’s eyes glaze over. But everyone understands the idea of e.g. a slow six-inch (150 mm) decade-long swing in the sea level, and that is what I am measuring directly and experimentally. Which me leads to …
… the fourth advantage, which is that the results are in the same units as the data. This means that a slow Fourier transform of tidal data gives answers in mm, and an SFT of temperature data (as in Figure 1) gives answers in °C. This allows for an intuitive understanding of the meaning of the results.
The final and largest advantage, however, is that the SFT method allows the calculation of the actual statistical significance of the results for each individual cycle length. The SFT involves fitting a sine wave to some time data. Once the phase and amplitude are optimized (fit) to the best value, we can use a standard least squares linear model to determine the p-value of the relationship between that sine wave and the data. In other words, this is not a theoretical calculation of the significance of the result. It is the actual p-value of the actual sine wave vis-a-vis the actual data at that particular cycle length. As a result, it automatically adjusts for the fact that some of the data may be missing. Note that I have adjusted for autocorrelation using the method of Nychka. In Figure 1 above, results that are significant at the 90% threshold are shown in color. See the note at the end for further discussion regarding significance.
Finally, before moving on, let me emphasize that I doubt if I’m the first person to come up with this method. All I claim is that I came up with it independently. If anyone knows of an earlier reference to the technique, please let me know.
So with that as prologue, let’s take a look at Figure 1, which I repeat here for ease of reference.
There are some interesting things and curious oddities about these results. First, note that we have three spatial scales involved. Armagh is a single station. The CET is a three-station average taken to be representative of the country. And the Berkeley Earth and HadCRUT4 data are global averages. Despite that, however, the cyclical swings in all four cases are on the order of 0.3 to 0.4°C … I’m pretty sure I don’t understand why that might be. Although I must say, it does have a certain pleasing fractal quality to it. It’s curious, however, that the cycles in an individual station should have the same amplitude as cycles as the global average data … but we have to follow the facts wherever they many lead us.
The next thing that I noticed about this graphic was the close correlation between the Armagh and the CET records. While these two areas are physically not all that far apart, they are on different islands, and one is a three-station average. Despite that, they both show peaks at 3, 7.8, 8.2, 11, 13, 14, 21, 24, 28, 34, and 42 years. The valleys between the peaks are also correlated. At about 50 years, however, they begin to diverge. Possibly this is random fluctuations, although the CET dropping to zero at 65 years would seem to rule that out.
I do note, however, that neither the Armagh nor the CET show the reputed 60-year period. In fact, none of the datasets show significant cycles at 60 years … go figure. Two of the four show peaks at 55 years … but both of them have larger peaks, one at 75 and one at 85 years. The other two (Armagh and HadCRUT) show nothing around 60 years.
If anything, this data would argue for something like an 80-year cycle. However … lets not be hasty. There’s more to come.
Here’s the next oddity. As mentioned above, the Armagh and CET periodograms have neatly aligned peaks and valleys over much of their lengths. And the Berkeley Earth periodogram looks at first blush to be quite similar as well. But Figure 2 reveals the oddity:
Figure 2. As in Figure 1 (without significance information). Black lines connect the peaks and valleys of the Berkeley Earth and CET periodograms. As above, the length of each periodogram is half the length of the dataset.
The peaks and valleys of CET and Armagh line up one right above the other. But that’s not true about CET and Berkeley Earth. They fan out. Again, I’m pretty sure I don’t know why. It may be a subtle effect of the Berkeley Earth processing algorithm, I don’t know.
However, despite that, I’m quite impressed by the similarity between the station, local area, and global periodograms. The HadCRUT dataset is clearly the odd man out.
Next, I looked at the differences between the first and second halves of the individual datasets. Figure 3 shows that result for the Armagh dataset. As a well-documented single-station record, presumably this is the cleanest and most internally consistent dataset of the four.
Figure 3. The periodogram of the full Armagh dataset, as well as of the first and second halves of that dataset.
This is a perfect example of why I pay little attention to purported cycles in the climate datasets. In the first half of the Armagh data, which covers a hundred years, there are strong cycles centered on 23 and 38 years, and almost no power at 28 years.
In the second half of the data, both the strong cycles disappear, as does the lack of power at 28 years. They are replaced by a pair of much smaller peaks at 21 and 29 years, with a minimum at 35 years … go figure.
And remember, the 24 and 38 year periods persisted for about four and about three full periods respectively in the 104-year half-datasets … they persisted for 100 years, and then disappeared. How can one say anything about long-term cycles in a system like that?
Of course, having seen that odd result, I had to look at the same analysis for the CET data. Figure 4 shows those periodograms.
Figure 4. The periodogram of the full CET dataset, and the first and second halves of that dataset.
Again, this supports my contention that looking for regular cycles in climate data is a fools errand. Compare the first half of the CET data with the first half of the Armagh data. Both contain significant peaks at 23 and 38 years, with a pronounced v-shaped valley between.
Now look at the second half of each dataset. Each has four very small peaks, at 11, 13, 21, and 27 years, followed by a rising section to the end. The similarity in the cycles of both the full and half datasets from Armagh and the CET, which are two totally independent records, indicates that the cycles which are appearing and disappearing synchronously are real. They are not just random fluctuations in the aether. In that part of the planet, the green and lovely British Isles, in the 19th century there was a strong ~22 year cycle. A hundred years, that’s about five full periods at 22 years per cycle. You’d think after that amount of time you could depend on that … but nooo, in the next hundred years there’s no sign of the pesky 22-year period. It has sunk back into the depths of the fractal ocean without a trace …
One other two-hundred year dataset is shown in Figure 1. Here’s the same analysis using that data, from Berkeley Earth. I have trimmed it to the 1796-2002 common period of the CET and Armagh.
Figure 5. The SFT periodogram of the full Berkeley Earth dataset, and the first and second halves of that dataset.
Dang, would you look at that? That’s nothing but pretty. In the first half of the data, once again we see the same two peaks, this time at 24 and 36 years. And just like the CET, there is no sign of the 24-year peak in the second hundred years. It has vanished, just like in the individual datasets. In Figure 6 I summarize the first and second halves of the three datasets shown in Figs. 3-5, so you can see what I mean about the similarities in the timing of the peaks and valleys:
Figure 6. SFT Periodograms of the first and second halves of the three 208-year datasets. Top row is Berkeley Earth, middle row is CET, and bottom row is Armagh Observatory.
So this is an even further confirmation of both the reality of the ~23-year cycle in the first half of the data … as well as the reality of the total disappearance of the ~23-year cycle in the last half of the data. The similarity of these three datasets is a bit of a shock to me, as they range from an individual station to a global average.
So that’s the story of the SFT, the slow Fourier transform. The conclusion is not hard to draw. Don’t bother trying to capture temperature cycles in the wild, those jokers have been taking lessons from the Cheshire Cat. You can watch a strong cycle go up and down for a hundred years. Then just when you think you’ve caught it and corralled it and identified it, and you have it all caged and fenced about with numbers and causes and explanation, you turn your back for a few seconds, and when you turn round again, it has faded out completely, and some other cycle has taken its place.
Despite that, I do believe that this tool, the slow Fourier transform, should provide me with many hours of entertainment …
My best wishes to all,
w.
As Usual, Gotta Say It: Please, if you disagree with me (and yes, unbelievably, that has actually happened in the past), I ask you to have the courtesy to quote the exact words that you disagree with. It lets us all understand just what you think is wrong.
Statistical Significance: As I stated above, I used a 90% level of significance in coloring the significant data. This was for a simple reason. If I use a 95% significance threshold, almost none of the cycles are statistically significant. However, as the above graphs show, the agreement not only between the three independent datasets but between the individual halves of the datasets is strong evidence that we are dealing with real cycles … well, real disappearing cycles, but when they are present they are undoubtedly real. As a result, I reduced the significance threshold to 90% to indicate at least a relative level of statistical significance. Since I maintained that same threshold throughout, it allows us to make distinctions of relative significance based on a uniform metric.
Alternatively, you could argue for the higher 95% significance threshold, and say that this shows that there are almost no significant cycles in the temperature data … I’m easy with either one.
Data and Code: All the data and code used to do the analysis and make these graphics is in a 1.5 Mb zipped folder called “Slow Fourier Transform“. If you change your R directory to that folder it should all work. The file “sea level cycles.R” is the main file. It contains piles of code for this and the last two posts on tidal cycles. The section on temperature (this post) starts at about line 450. Some code on this planet is user-friendly. This code is user-aggressive. Things are not necessarily in order. It’s not designed to be run top to bottom. Persevere, I’ll answer questions.
Bart says:
May 9, 2014 at 6:20 pm
What you’ve shown is a depth-2 autoregressive model which is right on the edge of instability. You can model and plot it in one line in R like this:
And indeed, the results are indistinguishable from yours.
However, if you modify the parameters even slightly, the series becomes unstable:
What you are showing in your linked file is an oddity that occurs when ar1 is almost 2 and ar2 is almost -1. The procedure is operating right on the edge of instability. If the parameters change even slightly, you get runaway positive feedback. And just below that runaway feedback threshold, you get the kind of cyclical behavior shown in your link.
As a result, it is extremely unlikely that this situation exists in nature. Nature doesn’t run that near to positive feedback. If it did, it would have gone off the charts long ago.
And in fact, if we analyze e.g. the CET as a depth-2 AR model, we get this:
Call: arima(x = CET, order = c(2, 0, 0)) Coefficients: ar1 ar2 intercept 1.2980 -0.5924 9.2213 s.e. 0.0123 0.0123 0.1163Note that ar1 is nowhere near 2, and ar2 is nowhere near -1. So it is far from instability, and therefore it doesn’t give results that look anything like what you show. Instead the results from such a model are somewhat lifelike.
However, pure AR models do poorly at replicating natural datasets. I generally use ARMA models of depth 2, which usually give much more lifelike results than pure AR models of any depth.
As you say, if you are interested, I will be glad to discuss this further with you …
w.
1sky1 says:
May 9, 2014 at 5:28 pm
I sure hope you don’t think this kind of egregious boasting improves your traction here … when a man comes on strong like that, from then on his stock tends to sell at a deep discount. There are people reading this blog who are both smarter and better educated regarding Fourier analysis than either you or I. I doubt that they are “blind to” the “certain basics” you go on to list.
I don’t think I’m blind to them either, but then I was born yesterday, que se yo? Well, except for the autocorrelation question. You say:
Here’s the ACF of the CET data …
Not sure how that is a vital means to distinguish between periodic and random signals. Your comments appreciated.
w.
Willis Eschenbach says:
May 10, 2014 at 12:52 am
“As a result, it is extremely unlikely that this situation exists in nature.”
You couldn’t be farther from the truth. This is exactly how nature works. Ring a bell. Slosh a cup of coffee. Drive over a pothole. Throw a rock in a pond. Twang a guitar string. How long the oscillations last depends on how pure the metal, how viscous the fluid, how old your shock absorbers, how big your rock and deep and wide the pond, how good your strings and stiff your guitar neck.
Nature abounds with such examples. Resonances are ubiquitous. Excite them with a persistent broadband input, and they will vary chaotically in amplitude and phase. But, there is structure. And, that structure can be used for estimation, prediction, and control.
Nature is stable because energy is conserved. But, it does not have to dissipate rapidly. And, the rate of dissipation is what determines how long oscillations which are excited will persist, and how pure they will appear over a given timeline.
I’m not going to argue this with you. I know very well it is true. I deal with it every goddammed day. So, if you want to be an ass about it, go ahead. I’m done.
Bart says:
May 10, 2014 at 2:36 am
It appears you totally missed my point, likely my lack of clarity. My reference to “this situation” meant an AR system such as you described, which is running right on the edge of positive feedback.
Yes, nature has resonances and feedback and damping systems and all the stuff you talk about.
But what nature doesn’t have are systems running on the edge of runaway positive feedback. If it did, such systems would have run away long ago. Since they haven’t, we can assume they don’t exist.
Oh, I see now. You’re right. And you know it very well. So obviously, there is no use in discussing it, you’re done …
Is that what’s called the “scientific method” on your planet? Is that what’s called a “conversation”?
Look, Bart, you set out to school me with some non-representative, non-realistic equation with odd properties that you thought I wouldn’t understand.
Instead of that happening, I pointed out that your equation was a simple two-lag AR model. Not only that, I pointed out that if you change the parameters by a few parts in a thousand, your example goes into runaway feedback. I noted that the only reason it had such variations in amplitude was that it was delicately balanced right on the edge of runaway, and that as such, it doesn’t represent a natural system. I added that such 2-lag AR models don’t look very lifelike in any case.
And now, rather than deal with the fact that you set out to school me and you got schooled instead, now you say you’re “done”. You want to stop the discussion, claim you are 100% right and by golly you know it, so you’re gonna take your ball and go home?
We should all be so lucky.
I could only wish it were true that you are “done”, Bart. But sadly, I’m sure you’ll be back to try your tricks again ..
w.
Bart –
It’s not that the specific resonator you proposed is “wrong”, (although as I tried to gently suggest, above, it is a schoolboy exercise) but that such a high-Q does NOT generally exist in nature. It CAN mathematically exist, and CAN exist digitally (numerically) but not if built of springs, capacitors, biological tissues, etc. A few accidental, small-scale, high-Q resonators (perhaps a cave as a Helmholtz resonator) might be noted.
As Willis suggested correctly, a continuous time resonator of that Q would wander into oscillation. Eighty years ago people had “regenerative” radio receivers that added a positive feedback component. In attempting to “pull in” a weak station, the user set the regeneration (“‘tickler coil”) so high that it not only drifted into a local squeal, but broadcast to the neighborhood!
Willis Eschenbach says:
May 10, 2014 at 3:32 am
“But what nature doesn’t have are systems running on the edge of runaway positive feedback.”
Surely, you jest. What do you think a resonance is? You don’t get any runaway condition because energy is not spontaneously created. But, natural systems can store energy with low losses, and does quite frequently, resulting in sustained oscillations. Drive that resonance with persistent inputs whose bandwidth extends across that resonance, and you get persistent wobbles.
“Look, Bart, you set out to school me with some non-representative, non-realistic equation with odd properties that you thought I wouldn’t understand.”
No, Willis. You’re being paranoid. I set out to explain to you why apparent oscillations seem to come and go in terms you would understand.
Bernie Hutchins says:
May 10, 2014 at 9:37 am
I wasn’t proposing any specific resonator, just attempting to demonstrate a concept. But I assure you, high-Q systems abound in nature.
I’d go into more depth, but why play to a hostile audience? I really don’t know what burr Willis has gotten under his saddle with me, but this I don’t need.
One last thing – you think the Q of the system I gave you is unreasonably high? It had 1% damping, which is higher than we generally assume for worst case structural modes. How about this one?
y(k) = 1.608*y(k-1) – 0.9875*y(k-2) + w(k)
Better? it’s essentially the same system. It has a different gain, hence RMS, but its poles are mapped from the same continuous transfer function poles, just sampled at a slower rate.
Perspective. It matters a lot when you are dealing with a massive system like the Earth with very slow dynamics.
Bart said:
“I wasn’t proposing any specific resonator, just attempting to demonstrate a concept. But I assure you, high-Q systems abound in nature.
I’d go into more depth, but why play to a hostile audience?”
Bart – I already suggested a cave as a possible Helmholtz resonator. But if such high-Q systems “abound” why not mention a few to humor us? I can’t, offhand, imagine any others aside from possible biologically-guided ones. I am not hostile – just honestly inquisitive.
Bart also said:
“How about this one? y(k) = 1.608*y(k-1) – 0.9875*y(k-2) + w(k) Better? it’s essentially the same system. It has a different gain, hence RMS, but its poles are mapped from the same continuous transfer function poles, just sampled at a slower rate.”
Bart – this is just changing the peak from 0.01 times the sampling frequency to 0.1 times the sampling frequency. Isn’t that all? You didn’t say previously that you were “mapping” a continuous function. What is your map? What is the analog prototype? In particular, what is the Q of the prototype? No zeros so it looks like it might be impulse invariant ? But I miss the significance of the change.
Bernie Hutchins says:
May 10, 2014 at 11:19 am
“…why not mention a few to humor us…”
I did. You seem to be laboring under the assumption that I am proposing cesium oscillator type Q-factors. Hardly. The Q factor of the system I gave as an example is about 50. It didn’t even need to be that high. Just high enough that oscillations can be self-sustained for maybe a cycle or two.
How high does that need to be? For a continuous system with characteristic function s^2 + (wn/Q)*s + wn^2, the time constant is tau = 2*Q/wn = Q*Tn/pi, so the ratio of the time constant to the period is tau/Tn = Q/pi. Let tau be, say, 5 cycles. Then, Q = 5*pi = a little over 15.
“…this is just changing the peak from 0.01 times the sampling frequency to 0.1 times the sampling frequency.”
Yes. The sample frequency went from 1 time unit to 10 of them. It is a zero order hold equivalent, z = exp(s*T), where T is the sample period. It’s not precise, because I rounded the coefficients, but close enough. The frequency was wn = 2*pi*0.01 rad/time-unit. The damping ratio was zeta = 0.01, hence Q = 50.
To Willis:
I am well aware that you are having quite a time with people making nit-picky arguments against you on other threads right now, and you are huddled in a defensive crouch. But, my input here is meant to be constructive, not condescending or supercilious, and possibly open up other avenues you might like to explore when you have the time.
And yes, Bernie, the actual zero-order hold equivalent has a zero. I didn’t bother with that – I just wanted to demonstrate a simple recursion which would produce lightly damped oscillations.
For comparison, a cesium fountain has a Q of about 10 to the 10th power. We’re not talking especially, or even moderately, high Q’s here.
… relative to that.
Bart gave these examples:
“This is exactly how nature works. Ring a bell. Slosh a cup of coffee. Drive over a pothole. Throw a rock in a pond. Twang a guitar string. How long the oscillations last depends on how pure the metal, how viscous the fluid, how old your shock absorbers, how big your rock and deep and wide the pond, how good your strings and stiff your guitar neck.”
Bart – Thanks, I saw those but did not suppose they were intended as examples any more than the digital resonator you gave. The bell, the suspension system, and the guitar string are not found in Nature, which is why I said examples with “biologic guidance” (even engineers are biological!) were possible. (Are you confusing Natural with Physical? I thought this was more about climate sort of things.) I am not sure of the rock in the pond – is that really a resonance (?) – or a radial wave phenomenon. Coffee cup – not too sure but the Q is probably 2-3 at most – my experiment using log-decrement. Possibly you should accept my Helmholst-Cave and seek similar if there are any.
Your analog/digital system transformation seems “non-standard” but I still don’t have any idea why you needed to offer the filtered noise example.
Bernie Hutchins says:
May 10, 2014 at 12:33 pm
Bernie, this is silly. For crying out loud, bang on a hollow log instead of a bell. Twang a twig instead of a guitar string. Is the rock in the pond exciting a resonance? Certainly. The motion is the superposition of the excited normal modes of oscillation with characteristic viscous damping.
The Q factor is roughly the number of cycles you will see after application of impulsive forcing before they damp out. Let’s see, I’ve got a coffee cup right here. Tap the side… so. 1, 2, 3… 14, 15…20, 21… I get at least 30 oscillations before I can’t see them anymore.
Nonstandard? Say what?
Why did I offer it? Because this is how nature works. Resonances abound. Excite them persistently, and oscillations will grow and fade over time. Sometimes they will be evident, and sometimes not. But, that does not mean that they do not have a physical cause which can be modeled to project, at least into the near future, how they will evolve over time.
Willis was suggesting that the oscillations he observed were a fluke with no physical significance, because they were not sustained over time. But, the conclusion of no physical significance does not follow from the observation.
Everything driving the climate is subject to these kinds of dynamics. The oceans oscillate in their basins. Kelvin waves propagate through the atmosphere. The tides are a manifestation of particularly high Q oscillations of the Moon about the Earth, and the Earth about the Sun, with the damping displayed in the orbital recessions.
The Sun, itself, displays such behavior. See, for example, here and the three other previous pictures (click on the left arrow at the side of the picture to navigate backward).
WordPress seems to have truncated my comment here. The missing part is reproduced below. I’ll try to answer any substantive questions that may arise next week, but all ad hominem comments will be summarily ignored.
5. The auto-correlation function (acf) of signals usually found in
nature provides a vital means to distinguish between periodic and random
signal components. In the former case, the acf is likewise periodic; in
the latter case, it decays with lag, i.e., is transient. Unlike periodic
signals, random signals are far from being perfectly predictable. It is
the rate of decay of acf, which is inversely proportional to signal
bandwith in the frequency domain, that determines the time-horizon over
which useful predictions can be made by optimal prediction filters. The
cornerstone of random signal analysis is NOT the FFT, but is provided by
the Wiener-Kintchine theorem, which defines the power density spectrum as
the F. integral of the decaying acf. Decimation (in time or in frequency)
algorithms must be employed to make periodograms provided by FFT analysis
approximate the power spectrum.
6. Sinusoidal regression is a least-squares curve-fitting method that is
legitimate ONLY in the case of truly period signals. It can be used very
effectively to determine the amplitude and phase when the period is known A
PRIORI–even in noisy records shorter than a single cycle (as described by
Munk in “Super-resolution of the tides”). The periodicity need not be
commensurable with the data sampling interval. However, even in
principle, it cannot provide a VARIANCE-PRESERVING decomposition of a more
general, random-signal data record, or when the signal period is unknown,
because it provides no set of ORTHOGONAL basis functions. When added up, the
discrete sinusoids thus identified generally do NOT add up to the data
series, as with the FFT.
The “Willisgram” simply sweeps a LS curve-fitting routine through
a comb of periods dictated by the sampling rate under the mistaken
assumption that the maximum amplitude obtained identifies the most
prominent periodicity. It doesn’t even identify the (perhaps
incommensurable) sinusoid which minimizes the residual in the data series.
Nor does the absence of a large amplitude sinusoid in such analysis of
short records prove the absence of irregular, quasi-periodic oscillations
commonly found in nature. The only “advantage” it offers is that of
satisfying naive presuppositions of what geophysical signals are composed.
This merely perpetuates the long tradition in the dismal field of “climate
science” of tinkering up quaint, ad hoc data analysis methods without
credible analytic foundations.
(1) Bart Said:
“Bernie, this is silly. For crying out loud, bang on a hollow log instead of a bell. Twang a twig instead of a guitar string. Is the rock in the pond exciting a resonance? Certainly. The motion is the superposition of the excited normal modes of oscillation with characteristic viscous damping.”
Bart –
Well I think the pebble in the pond is a different problem – as I recall – probably because of the lack of boundary conditions as we would have in, for example, your plucked string. The log and the twig are biologically guided. Still – my Helmholtz cave seems the best example. No matter – I need your “natural” log in a moment.
(2) Bart also said:
“The Q factor is roughly the number of cycles you will see after application of impulsive forcing before they damp out. Let’s see, I’ve got a coffee cup right here. Tap the side… so. 1, 2, 3… 14, 15…20, 21… I get at least 30 oscillations before I can’t see them anymore.”
Bart –
Not much of a way to do science! One method of measuring Q is to measure the “ring time”, tr, (to 1/e or about 37%). Then Q=pi*fo*tr. So if we went to about (1/e)^3 or about 5%, 3*tr, Q would be about the number of oscillations – which is where your rule of thumb apparently comes from. Good luck stopping at 5%. Is 5% “observable” or not. Depends. At what point did you stop? More to the point use instead Q = -pi/Ln(d) where Ln(d) is the “natural log” (see it was handy) of the decrement (the change in amplitude from one cycle to the next). This is what I used with the coffee cup. Mine died very quickly. I don’t see how you could have possibly seen 30 oscillations! But I did find that there are at least three resonant modes you can get – so this is a bad example.
(3) Bart also said:
Nonstandard? Say what?
Bart –
Nonstandard because people don’t use the direct map z=exp(sT) because it doesn’t guarantee anything (preserve anything from the analog response). Mostly Bilinear-z and IIV are used, occasionally forward and/or backward difference, but not the matched z. That was out, or should have been, about 40 years ago. Not that it matters since you were just making things up as an example.
(4) Bart –
Yes this is silly. What we need to get to is the fact that you can’t name and/or give evidence for any CLIMATE parameter that has a record suggesting a Q of 50 – by your own rule of thumb – 50 observations, decaying. Surface temperature, ocean heat, ice cover, tornadoes, sea level? No reason to expect simplicity.
1sky1 –
I read through your addendum to your earlier post. Likely you will agree that even a moderate-length text-posting such as yours is not an adequate substitute for a complete technical presentation with equations and figures, including well-defined test examples. Of course, that would be a week’s work or so.
One simple thing that is not clear to me is what you mean by:
“ The “Willisgram” simply sweeps a LS curve-fitting routine through
a comb of periods dictated by the sampling rate under the mistaken
assumption that the maximum amplitude obtained identifies the most
prominent periodicity. It doesn’t even identify the (perhaps
incommensurable) sinusoid which minimizes the residual in the data series.
Nor does the absence of a large amplitude sinusoid in such analysis of
short records prove the absence of irregular, quasi-periodic oscillations
commonly found in nature.”
The question is: Does it not “identify” the sinusoidal because it is totally wrong in concept, or is it just imprecise; perhaps even being forced to select from widely spaced choices with a “real” answer in between?
Thanks
Bernie Hutchins says:
May 10, 2014 at 5:57 pm
Wow. Sorry, no. No. No. No. No. Not much point in continuing this dialogue.
Bart said, May 11, 2014 at 10:58 am
“Wow. Sorry, no. No. No. No. No. Not much point in continuing this dialogue.”
Bart –
You are of course free to discontinue. It might show a bit more class (and increase your credibility) if you FIRST responded to my comments, in line with whatever ability you have, and THEN asked to stop. Otherwise it could just be a cop-out!
Your linking to Wikipedia probably was an indication of your level of experience that I should have noted, rather than gotten more technical, and I apologize for that.
Sincere best wishes and thanks for your time.
Bernie
Bernie Hutchins:
The figure of merit in “sinusoidal regression,” as in any LS curve-fitting
scheme, is NOT the amplitude of the fitted sinusoid, but the variance of
the residual. The “Willisgram” shows the former for a PRE-SELECTED discrete
set of periods coinciding with INTEGER multiples of the data-sampling
interval, but fails to show the latter. This leaves unresolved the SINGLE
sinusoid amongst all periods, possibly INCOMMENSURABLE with the sampling
interval, that best fits the available data in the LS sense.
But more importantly, even that best-fitting sinusoid is not necessarily a
good–or even reasonable–representation of the continuous-time signal that
gave rise to the data. Only in cases where the signal actually contains
strictly periodic components, such as the annual cycle, can sinusoidal
regression be expected to provide meaningful results. In all other cases,
it leads to badly mistaken expectations, as the prognosticators of solar
sunspots rudely discovered with SC24. That is the consequence of mistaking
a narrow-band random process (with strong structural similarities to that
of shoaling swell) for a periodic signal suitable for analysis by repeated
sinusoidal regression. [BTW, for those interested in distinguishing the
difference between truly periodic and narrow-band processes, I can e-mail
the estimated acfs of the monthly S.F. sea-level and of the yearly Zurich
SSN data to any address.]
It needs to be recognized that, because of frequency entanglement in
non-orthogonal decompositions (particularly of long periods in short
records), a wide range of linearly incremented periods will show large
amplitudes that would NOT be there in much longer records. That is the
irretrievable failure of the conception behind the “Willisgram.”
Undoubtedly, the brief sketch of the landscape of proper signal analysis
methods I provided here is wholly inadequate to the task of elucidating a
complex subject to analytic novices. But that’s what a plethora of
introductory texts written since the 1960s undertake. A popular one is
Oppenheim & Shafer’s “Discrete-time Signal Processing.” It’s dismaying that
basic concepts covered there continue to elude so many who feature
themselves “climate scientists.” And it’s disturbing that the ANNUAL
cycle, which is HIGHLY prominent in both sea-level and temperature data, is
nowhere displayed in any results on this and previous related threads.
1sky1 –
Thanks for the reply.
I understand what you are saying. But I think my most fundamental doubt is that there is enough “Fourier” signal in the climate data that any “detection” means much. I think this is your view too?
But it does seem to me that you can get the frequency pretty accurately using Prony’s method (or similar) combined with some least square manipulations to handle the noise that Prony can’t tolerate directly, and when the system has a dominant (not necessarily high-Q) pole pair. This likely works even when an FFT lies through its teeth! NOT saying that it helps with climate data.
I first learned DSP using a 1972 mimeograph version of O&S! I co-taught with the “Halloween” (black/orange covers) versions, and later versions. I didn’t think the students liked it very much, and I never used any text when teaching by myself. Certainly I don’t think you were suggesting that that text was very useful for what you have been saying? Not to say that some very BASIC misunderstandings are not common in the comments. The basic one hurt. The more subtle ones get all of us from time to time (i.e., every day!)
Bernie Hutchins:
To be sure, detection and proper identification of oscillations in the multi-decadal to centennial range from climate records seldom much longer than a century is a real challenge. That challenge is scarcely helped by applying data analysis methods ill-suited to the intrinsic character of the underlying signal. Sadly, confusion between periodic and aperiodic signals is rampant in “climate science,” whose mind-set is rooted in the data = mean + trend + cycle + i.i.d. noise paradigm of classical statistics, with understanding of stochastic processes imprisoned in simple Markov chains. My purpose here was to draw a sharp distinction between periodic, transient, and far more complex random signals encountered in geophysics and to indicate suitable analysis frameworks for each.
It strikes me that Prony’s method is applicable only to transient (finite energy) signals, but not to continuing (finite power) signals, periodic or random. In the latter case, when there no strictly periodic components, the very idea of using LS curve-fitting methods on two or three wave-forms is foolish prima facie. Even “surfer dudes” recognize that no two waves are alike in height and apparent “period.” Few in “climate science” manifest such wisdom. We are truly at the threshold of reliable, meaningful characterization of long-period oscillations. Nevertheless, unlike curve-fitiing routines, sound signal analysis methods provide the best glimpse available of those signal characteristics.
You’re entirely correct that O&S’ latest opus was not my reference nor my recommendation. It is simply a cheaply available text that covers some of bare basics. At least the authors had learned enough from their sabbatical at WHOI to avoid the lapses that made their earlier “Halloween” opus indeed frightful. I wrote the commnent off the top of my head, deliberately refraining from recommending “Signal Analysis” by Papoulis, which is intended for a more advanced audience. BTW, I began hands-on learning as a summer undergraduate apprentice to one of Wiener’s students, who had his own copy of the “yellow peril,” and in loaded up in grad school on “Statistical Theory of Communication” along with geophysical course-work. Always nice to meet someone who has solid grasp of the subject.
1sky1 –
I myself have a detailed and complete understanding of the earth’s climate system; which curiously has an exact counterpart detailed and complete understanding of my vegetable garden: “Sometimes things grow, and sometimes they don’t”.
In a similar vein, I tend to have to discover things for myself. I am not offering that necessarily as a virtue, but as the easiest and perhaps only way I myself can learn. And I’m not too good at learning something until I really need it. Give me someone’s basic idea (or a notion of my own) and let me play with Matlab writing my own code. “The purpose of computation is insight, not numbers.” – Richard Hamming
Prony’s method works for any decay constant, including no decay. If you take four samples like [0, 0.5 sqrt(3)/2, 1] Prony happily tells you the frequency is 1/12 if you specify that it is 2nd order and there is no noise. The FFT lies, and never heard of a frequency 1/12 anyway. Prony is just as happy with a frequency 1/12.123, for example. And it is linear. A perhaps limited, but a valuable tool. I’m still trying to see what it can do. I learned Prony from Tom Parks who remarked that “people are always re-inventing Prony’s method”. At that time, I noted that the methods was coming up on 200 years old (invented in 1795).
Yes I love all the signal analysis and modeling stuff. But many years ago in a senior physics lab I fit some data to every possible math function I could imagine, and plotted every possible graph. My professor, Herb Mahr, praised my industry in doing so, and ever so gently (bless his heart) told me that it didn’t really MEAN anything.
Nothing new under the sun.
Bernie
Bernie:
There indeed is little new under the sun–except in the minds of the “way
cool” ME-generation that treats what isn’t presented on the Web as de facto
non-existent. This, of course, leaves the field wide open for puerile
claims of “invention,” which you’re mature enough to see through.
It turns out that Pony’s method has been tried out on transient geophysical
signals found in seismograms and sea-level disturbances by tsunamis (see:
http://www.reproducibility.org/RSF/book/tccs/nar/paper.pdf and
http://oos.soest.hawaii.edu/pacioos/about/documents/technicalreports/TechRep_PacIOOS_SL_Forecast_07-17-2012.pdf).
If I read correctly (with only ~1/2hr available in my schedule for Web
browsing), in effect it fits an ARIMA model to the data, reminiscent of EOF
decomposition in “singular” spectrum analysis, or fitting a purely AR model
to the sample acf in Burg’s “maximum entropy” spectral estimation algorithm,
thereby extending it to longer lags than the short data record can
robustly estimate. I.e., it’s an attempt to replace genuine empirical
information with a presumptive model.
What many decades of experience in empirically-based research has taught me
is that all attempts to find free lunch ultimately leave the serious
scientific appetite hungry. When the rubber of preconception meets the
hard road of reality, it’s the rubber that usually breaks. But that is rare
in “climate science,” where data are much too often manufactured, or
presented, to conform to preconception. An example of this is the sample
acf of the exceptionally noisy monthly CET, shown by Willis on a
LOGARITHMIC scale and only for lags of 2yrs and longer–a transparent
attempt to induce the naive conclusion that long-term temperature series
conform to the exponential decay of a “red noise” process. That’s how weeds
from the garden, masterfully promoted on the Web, get sold as roses.