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.
Thanks Willis. Good explorations.
You said “Looking for regular cycles in climate data is a fools errand.”
But the imprints of the 24 hrs, seasonal and orbital annual cycles will be there. Also the Sun will mark the record with its “around” 11-years cycle. ENSO and PDO will mark it.
Having a cyclic nature is not the same as having a regular cyclic nature.
Nick Stokes says:
The Fourier Transform goes back to 1822. It takes the form:
F(ω) = ∫ f(t) cos(ωt+b) dt
where b is a phase. You can drop b; it’s enough to calculate a version with sin and cos to cover all phases (just use the sum formula for cos).
===
That’s a cosine transformation, the full fourier is
F(ω) = ∫ f(t) exp(iωt) dt ; i or j in E eng world is srqt(-1)
exp(iωt) = i.cos(ωt)+ sin(ωt)
FT ( or FFT ) is a complex spectrum containing both “cos” and “sin” terms which thus contains the phase information .
atan2(Re/Im) ; atan2 being arctan that knows the difference between the +ve and -ve quadrants.
BTW Nick, nice confirmation of the equivalence of SFT and padded FFT. Shows the results are a useful estimation of the spectrum for broken data.
I’m ultimately more comfortable with several short periods (with gaps) rather than one long period that has been smoothed. The short periods are real, the long period is unreal.
Tweaking numbers smacks of ‘climate science’. Inevitably you end up with results you want to see.
Science battles should be over apodization, not sampling of Internet polls.
There is a FFT algorithm called FOURG that allows for any N number of input points as long as N is non-prime. I’ve used N LOG N for so long I forgot. I think FOURG is also known as the Cooley – Tukey algorithm. I have it somewhere in FORTRAN IV so Gavin would be proud.
Greg says: May 4, 2014 at 5:11 am
“That’s a cosine transformation, “
Not with the arbitrary phase term. Try with b=π/2.
That’s my point there. Expand the cos and you get a lin comb of sin and cos transform. The latter are all you need.
Willis re figure 5: “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.”
Again you insert your interpretation in place of what the data shows. It does not show “no sign of” it shows a reduced amplitude at 24…25 years.
freq ave(20 33) = 24.9
modulation = 101.5y.
Now since your data sample for 2nd half is 102, there may be some data processing reason for this apparent modulation, or it maybe physical. In either case it is incorrect to conclude there is “no sign” of it. You are just not reading the signs correctly.
The combined amplitude of those three peaks is about the same at “first half” and “full”.
Doesn’t all this analysis show that Mannipulated temperature series strips out natural cycles? Do we have anything close to raw data to work with anymore?
It may be fun playing with software etc. but that is how Mann found a hockey stick. Please point my nose towards something real, like a temperature or something. Like something in these cycles that relate to real life.
@Greg.
I agree that you can approximate the spectrum using irregularly sampled data by quadrature or padding. The problem is that if you express the FFT in matrix form, its eigenvectors are orthogonal. When irregularly sampled, this is no longer true. The problem arises (in my experience) when you try to use the spectrum for Fourier interpolation. This has occurred quite often in my area where I sampled a “signal” spatially and the sampling depends on the geometry of the signal in question.
I’m still worried because I don’t understand what has been done here. There is reference to “fitting” and “optimisation”, which is not the case in a DFT which is essentially a correlation. If there has been “optimisation” of some objective function, then I suspect that what is being discussed is not a DFT (as is well known, least squares fit of co/sin waves by least sqares does not give the correct result).
Steve from Rockwood says:
“2. FFTs are limited to values of 2^n (16, 32, 64, 128 …).”
This is a commonly held belief but its not true. FFT’s can be computed using any prime numbers. The ultimate performance in this respect would be the MIT fftw code (which is what my company uses, under licence, for very large 3D fft)
And as I said in the previous post, what Willis doing is a Fourier analysis but by fitting each sine wave one at a time in the time domain. Its not new or original and its very inefficient.
I’ll ask a question by first summarizing what various folks have pointed out above.
(which has an average 
) and minimizing the quantity ![\sum\limits_{n}^{}\left[ S_n - \mu_S - r_k \sin(2\pi t_n /T_k + \phi_k)\right] ^2](https://s0.wp.com/latex.php?latex=%5Csum%5Climits_%7Bn%7D%5E%7B%7D%5Cleft%5B+S_n+-+%5Cmu_S+-+r_k+%5Csin%282%5Cpi++t_n+%2FT_k+%2B+%5Cphi_k%29%5Cright%5D+%5E2&bg=ffffff&fg=000&s=0&c=20201002)
for various sine frequencies 
by finding optimal values of 
and 
. You see as among this approach’s advantages that 
has the same units of measure as the data 
.
, with all frequencies 
so chosen that the sinusoids (
) have respective integer numbers of periods over the 
record and are thus orthogonal: to the extent that the discrete Fourier Transform approximates the continuous Fourier Transform, adding them together reduces the error.![\int\limits_{0}^{T_1}\left[ a(t)-\mu_a - r\left[ b(t)-\mu_b\right] \right]^2dt](https://s0.wp.com/latex.php?latex=%5Cint%5Climits_%7B0%7D%5E%7BT_1%7D%5Cleft%5B+a%28t%29-%5Cmu_a+-+r%5Cleft%5B+b%28t%29-%5Cmu_b%5Cright%5D+%5Cright%5D%5E2dt&bg=ffffff&fg=000&s=0&c=20201002)
, one sets to zero the derivative of that quantity with respect to 
and solves for 
, resulting in 
. So correlation is just a short cut to error minimization. And 
‘s units are the same as 
‘s if 
is dimensionless, so, despite the fact that folks are wont to distinguish between the time and frequency domains, I don’t see why you couldn’t assign the Fourier coefficients the input data’s units of measure.
and the sinusoidal component 
?
Although I didn’t have the patience to run down all the R functions you use, my guess is that you’re taking the (detrended) data
In apparent contrast, a DFT (which the FFT is an algorithm for computing efficiently) works by correlation
With the possible exception of the orthogonality constraint and your ignoring the errors that omitted data impose, though, these are actually the same things. To minimize
That was the wind-up. Here’s the pitch: In this scheme, what is the mathematical expression for the confidence interval? Can that be expressed in terms of the data
to Thinking Scientist:
I thought you could use any number, not just a prime one. The series of numbers you referred to are not prime
Willis,
Climate cycles don’t have clockwork periodicity. Even if the primary driver is some sort of astrophysical cycle with a stable periodicity, the Earth’s climate doesn’t react exactly the same way each time. On top of that, there are many convolved cycles with different periods and amplitudes. Atmosphere and oceanic circulations are very stochastic processes. In the HadCRUT series, the nominal 60-yr cycle can range from about 50 to 80 years.
If you are looking for highly regular cycles, you won’t see them. If you are looking for coherent patterns of quasi periodic fluctuations, they are hard to miss.
Have you tried your SFT on the GISP 2 ice core d18O series? It would be interesting to see how that compares to Davis & Bohling.
RS says: May 4, 2014 at 5:59 am
“If there has been “optimisation” of some objective function, then I suspect that what is being discussed is not a DFT (as is well known, least squares fit of co/sin waves by least sqares does not give the correct result).”
He’s fitting each sinusoid to the objective function of minimum sum squares deviation (pointwise), varying amplitude and phase. I showed here how that is equivalent to FT (mod of complex form). It gets more complicated if you don’t have orthogonality, as you say (ie not DFT). But the objective function itself is just a quadratic, so you can differentiate and get it in one shot.
I’ve run Willis’ code with test sinusoids. It works quite well – as well as other FTs. I would say that it is quite robust to missing data and not much affected by non-periodicity over the interval, as long as you stick to periods less than 1/3 total, as he does. It’s slow – about 20 secs on my PC for 200 yrs; FFT is sub-second.
You’d expect all that. FT’s with orthogonality do minimise that sum squares. And above the third harmonic, deviation from orthogonality due to ends isn’t too bad.
sorry Nick
that all sounded like IPCC talk with the if, buts and maybes. I designed an electronic system based on that and it didn’t do shit. I won’t design anything unless I have zero error bars or less than .001%
For all the people who know things about FFT, how well does it work for heteroscedastic data?
It’s easy to see that the errors in the data sets change with time, as seen in this image that I made a while back. [The image has some smoothing and fits that can be ignored. I made it back when I thought Lovejoy was saying that error as a function of time was constant and small, and I wanted to double-check that for myself.]
Regression can get around that by assuming that the error has the same distribution type, even if the parameters of that change over time. It will screw up any error estimates, but the coefficients should be fine. I have no idea what this means for FFT, though, which is why I ask. Thanks!
to Indpndt:
I’m not much of an expert on FFT, just generally familiar with it’s principles. Essentially you take a curve/line and you can break it down to sine waves of various properties. What u decide based on those sine waves is up to you. You should be able to ‘reassemble’ those sine waves to the original curve/line. The points on your original diagram are not much use in this. Use the blue line.
Thanks, David Middleton. Good observation.
“the Earth’s climate doesn’t react exactly the same way each time.”
That, I think, is because our planet responds to the few regular cycles from ever changing points in its unfolding history.
A new winter starts after a new fall, so even if the astronomical conditions are the same, the response will be new.
So, climate is not theoretical-clockwork cyclical, maybe a clock with a loosely fitting mechanism that at times gains, at other times looses. And when El Niño winds it up real tight it shifts to a new regime that only later La Niñas will relax.
alex says:
May 4, 2014 at 1:25 am
Well … no. My SFT gives the same result as a DFT … but by an entirely different method.
Yes, and so does the SFT, only slower. There is no difference in the result. However, there is a large difference in the method.
So you claim. However, when I’ve actually tested my method by removing data and then putting it back it, I find it works quite well. I think it may relate to the fact that I’m NOT transforming the data from the time space to the frequency space and back again as is done by a Fourier transform. I arrive at the same point, but by a different method that as far as I can tell from my testing is NOT subject to the Gibbs effect. So what did you find when you actually tested my method? Oh … right …
Come back when you’ve actually used my method and you know what you are talking about. Right now, you’re just revealing your own poor math education by making claims about a mathematical method that you have never tested in your life.
w.
RS says:
May 4, 2014 at 3:56 am
Thanks, RS. I am indeed optimizing, meaning that I adjust the phase and the amplitude of the sine wave to get the best fit.
And whether I am “calculating Fourier components” or not is of little interest, because I arrive at the same endpoint. Go figure …
w.
Hi Willis
There seem to be two alexs on this site . I’m Alex with a capital A – I’m the nice one. I don’t denigrate, I just ask questions or make statements politely.
Seeing different cycle lengths at different times brings to mind what you see in a negative-feedback process control environment when trying to regulate an output under changing conditions. My experience is with PID tuning for regulating pressures/flows in a natural gas pipeline. It’s very easy to come up with control parameters for any given situation, but it can be very tricky to handle any and all conditions you are likely to encounter, especially as they change. Obviously, you must when safety and reliability are at stake.
In any particular set of conditions, there will always be an oscillation and its cycle length and amplitude will be fairly uniform. Change any condition that affects the system (upstream pressure, downstream pressure, downstream demand, varying set-points), and the system will stabilize with an oscillation with a different amplitude and/or cycle length. In my situation, things changed all the time based on temperature, time-of-day, season, etc.
What brought this to mind was Willis’s earth-thermostat hypothesis. Assuming that there is such a thing, if whatever affects climate never changes then you wouldn’t expect to see the length or amplitude of the resulting cycles change over time. The apparent fact that they did change implies that something, probably several things, did change and, in my mind, likely changes all the time.
In my closed-system environment, it was a bitch to come up with control algorithms that adequately handled the limited range of conditions likely to be encountered. I can’t even imagine what it would take to understand what’s happening with our planet – certainly not today’s global climate models. But in all cases, I would expect to see cycles like the ones Willis found and to see the length and amplitude change over time as well. With so many unpredictable inputs (volcanoes, cosmic rays, whatever), I don’t think you would ever be able to predict these, only see them after the fact.
RS says:
May 4, 2014 at 3:56 am
Thanks, RS. Since I’m not wandering into the frequency domain, I’m not sure why this would apply. Once again, all of you guys are making claims based on what the Fourier transform does. I’m doing something different and arriving at the same endpoint.
And yes, if you remove any data from a complete dataset, you will NOT calculate the true spectrum of the signal … the question is, how resilient is the method to gaps in the data. The standard Fourier transform is not resilient at all. The slow Fourier transform is surprisingly resilient.
Please, folks, test my method before making claims, no matter how theoretically correct you might think that they are. Yes, if there are gaps in the data the results will be different, in part simply because we have less data. But here’s an example of what I’m talking about. This shows the Armagh periodogram, along with a periodogram of the Armagh data with a full 5% of it removed at random, and also a periodogram of the Armagh data with a five-year chuck of it cut out from 1900-1905:
Can you do that with a standard Fourier transform?
You guys are making claims based on the traditional method of creating the periodogram. I’m not using that method.
w.