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.
Dan in Nevada says:
May 4, 2014 at 7:51 am
“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.”
I can imagine in a long pipeline, apart from variations in pressure and temperature induced by operation of transmission equipment, superimposed variations in volume/density of the oil plus changes in the dimensions of the pipe by ambient temperature and atmospheric pressure, changes in the slopes of the pipe, nearby traffic vibrations etc….that one could expect occasional, unpredictable coincidence of oscillation peaks would cause a serious problem.
Earth’s climate with global/celestial/geological/biological factors, all the oscillations of the sea and atmosphere, irregular contributions from volcanoes, earthquakes, hurricanes, tornadoes,…. yeah modelling and control are out of the question at the current state of knowledge.
In Willis’s work, looking at some of the comments of Greg and others, it seems to me the longer “rough” periods have the biggest effect on temperature over time and much of the short term ones are closer to something like noise saddled on these durable ones. Maybe 50 to 80 years is how grainy the significant oscillating signals are and the in between ones should not be seen as having disappeared.
A very good review of basics of Fourier Transform can be found in Numerical Recipes The Art of Scientific Computing.. This is an excellent text for many digital analysis techniques. The subroutines have been written in many different computer languages and are available. Care must be taken when using FFT that the data are treated properly, usually the endpoints of the raw data are tapered and care must be utilized to not alias data by using too crude sampling.
Consider the little “valley” between the 21 year and 24 year peaks in figure 4, top chart. Notice also that a little peak fills the valley just to the right at 28.
Whereas the chart just below it shows a nice peak at 22 and a nice null to the right at 28, also a peak at 8.
This spectrum is typical of frequency (or phase) modulation in radio engineering in a non-linear mixer. Normally radio waves ignore each other, just as ocean waves normally ignore each other.
But where a nonlinearity exists, these waves won’t ignore each other, each modifies the other, producing artifacts that reveal the nature of the nonlinearity and what went into it.
This nonlinearity can exist in the data itself and isn’t a natural phenomenon. If you simply close the gap of missing data, you have introduced a nonlinearity that will mimic phase or frequency modulation, splitting peaks into two (or more) components and creating sidebands that fill nulls.
In nature a big non-linearity is freezing water — a phase change of water but it also forces a phase shift in the temperature signal. As energy continues to dwindle, temperature dwindles until it starts forming ice — Central England is the name of the chart. In this case splitting peaks suggests a possible natural nonlinearity.
Therefore it is likely that a real 22 year phenomenon exists but is being modulated in a phase-shifting (frequency modulating) way to produce energy peaks at 21 and 24. Notice the peak at 8 is also split by what appears to be a proportionally similar amount. If your data is continuous then this is a natural nonlinearity; but if your data is NOT continuous then what you see is an artifact of missing — or ALTERED — data.
By the way, the technique can be used to discover edited photographs, where one photo has been spliced onto another. If the seams are blended their frequency (or period) spectrum will differ substantially at the seam where the two photos were blended; also, if different cameras were used or at different resolutions the DFT of the lateral image frequencies will differ in the pasted portion of the image. This itself can be imaged by assigning a color or brightness to the highest frequency found in any small range of pixels (such as JPEG’s 16×16 blocks). Interestingly, JPEG has already done the DCT (discrete cosine transform) and directly encodes the highest frequency found.
“But where a nonlinearity exists, these waves won’t ignore each other”
That’s good point to remember, FT et al are all built on the assumption of a linear system where a signal can be chopped up at will any group of components can be added to make a reconstruction of the signal. That does not mean the system has to be linear in it’s operation but at least the variable being analysed has to be linearly additive.
This was one of my problems with the Steucker and Timmerman paper last year based on spectral analysis of trade wind speeds. Wind speeds are not linearly additive and thus not susceptible to spectral analysis.
Since the kinetic energy of the wind is proportional to wind speed squared, this may be a better choice.
http://climategrog.wordpress.com/?attachment_id=283
Michael Gordon says: “This spectrum is typical of frequency (or phase) modulation in radio engineering in a non-linear mixer. Normally radio waves ignore each other, just as ocean waves normally ignore each other”.
Can you be more specific about that splitting and how it applies to those peaks. I noted those five peaks fitted rather well to being read as amplitude modulation.
http://wattsupwiththat.com/2014/05/03/the-slow-fourier-transform-sft/#comment-1628239
Now looking at the triplets in “full” panel of the graph and calculating the average _freq_ of the outer peaks:
21 and 41 => 27.7y 86y
24 and 34 => 28y 163y
Udar says: “This problem is well known in DFT processing, and it is dealt with by applying window function, which tapers endpoints of data set to 0 in a specific fashion, that depends on particular window function.
What is interesting about SFT is that it could be made immune to this problem, given proper selection of optimization function and proper constraints on the way it finds its solution.”
===
This may be a significant advantage of SFT for longer periodicities which are one of the main points of interest in climate.
Usually Fourier spectral analysis requires a tapering window fn which is quite destructive and distorting to the longer periods but gives more accurate results of the high frequency end. In electronics and audio this is often not a problem since the periods being examined are orders of magnitude shorter than the available samples. Sadly not so in climate.
I generally use an extended cosine window fn that is flat over 80% and just tapers the last 10% at each end. This usually works quite well once the autocorrelation has been removed with a diff, but it’s still a compromise.
Could you more specific about what you consider to be “proper selection” and “proper constraints”?
re:
21 and 41 => 27.7y 86y
24 and 34 => 28y 163y
What would be interesting is to find the phase of the implied 86 and 163 periodicities and get more accuracy than the crude nearest year values for 21,41 and 24,34. With more precise calculations, it should be clearer whether those longer periods harmonics since they will have large uncertainty as they stand.
There is a possibility a whole series of harmonics there : 21,42,84,168. Now there’s nothing magic about that but it would simply be the Fourier components of non-harmonic circa 168y form. Now in a 200y record there’s no reason to see that as a repetitive cycle on it own, it’s just the long term record.
What is interesting is it’s apparent interaction with the circa 28y periodicity. It seems that this long term variation is modulating the 28 y periodicity
It is interesting that I also found a circa 28y modulation in Pacific trade wind data linked above.
A.M. side-bands: 28.6583 4.4450 -28.3212
One thing that comes out of longer records, like BEST and these individual series is that 18th and 19th century showed much larger variability in temp than 20th c.
Now is the current reduction due to a damping of more CO2 radiative “blanket” effect or just a long term modulation? If it’s the former we should be glad of our atmospheric CO2, if not we can we expect the larger variability to return (irrespective of CO2) and with it more “extreme weather”?
I’ve no idea what the 28y period is about but if it’s in Armagh air temps and Pacfic trade winds , it gets interesting.
Martin Gordon says:
“But it is easy to understand. Not everyone is a mathematician. I am not. I can struggle with calculus but since this problem is to be solved numerically anyway, which by its nature is discrete/finite, it isn’t an exact representation of the integral anyway.”
Fair enough if this and the previous article were written as ” how to understand fourier transforms by simple time domain examples”. But they are not, they are written as though some new method of analysis were invented. And when I pointed out to Willis that this was effectively calculating the ft in the time domain he did not understand.
ThinkingScientist says:
May 5, 2014 at 12:18 am
Hogwash. I knew some unpleasant person would come along with this lame accusation, which is one of the reasons why I wrote:
So no, your claim just reveals that your animosity exceeds your ability to read.
w.
The normal way to do FFT on irregular data is to do sine or cubic interpolation to caret an ‘analogue’ time series which is then ‘resampled’ to get a times series to feed to the FFT.
What this SFT is, is exactly the same as using sine interpolations.
And I am afraid to say a cursory inspection of the other ‘differences’ doesn’t hold water either.
I am not decrying the fact that the results are indeed extremely interesting, but merely pointing out that reinventing a wheel doesn’t make it any the less round, in the end 😉
If I had the time I’d code up an algo to do sine interpolations and FFT with a choice of how the results are displayed for any arbitrary time/amplitude series.
RS:
I appreciate your response, but it explains no better than Mr. Eschenbach has been able to what he means by “significance.” Knowing my limitations in statistics, I’m perfectly willing to accept that everyone but me understands perfectly well what he meant. Still, I’ve seen no evidence that even he does.
As to his method generally, it is in my view no less valid in identifying a frequency component in a signal than a DFT is. There’s nothing holy about the particular selection of sinusoids that DFT uses if all you want to do is find a given sinusoidal component in a signal, and, outside of the sinusoid selection, his approach is exactly the same as the DFT: it essentially uses correlation.
The reason for the DFT’s particular sinusoid selection is that it’s more appropriate for analyzing a linear system’s response to a signal. Since a linear system’s response to a sinusoid is a sinusoid of the same frequency and (often) readily determined phase and amplitude, and since adding the DFT-determined components together reconstitutes the original signal, and since a linear system’s response to a sum of components is the sum of its responses to the individual components, the DFT can be used to compute a linear system’s response. To the extent that Mr. Eschenbach’s approach uses a sinusoid selection different from the DFT’s, it can’t be used in that way–but he never contended that it could.
So, while there’s no doubt plenty of truth among the various disputants’ statements, I personally see nothing further of interest other than what Mr. Eschenbach means by “significance.” And the prospect of seeing that clarified seem increasingly remote.
Joe Born says: May 5, 2014 at 2:53 am
‘I appreciate your response, but it explains no better than Mr. Eschenbach has been able to what he means by “significance.”’
In his recent version that he described to you, using lm(), the R regression reports standard error for the coefficients a and b. I think that might be what is meant. It involves assumptions of independence; the residuals after fitting are assumed to be iid. Their variance is estimated by the sample variance, and the variance of a and b just follow from the fact that they are a linear combination of the data. Since Willis end up with sqrt(a^2+b^2) there is a bit of extra calculation there.
The residuals are probably autocorrelated, not iid, so the confidence interval range will be understated.
The primitive wheel reinvented here is equivalent to arbitrarily limiting wavelet extent to a single value and deliberately ignoring every other possible value from across the entire extent spectrum for every single grain.
Why would someone artificially impose such infinitely vision-narrowing blinders?
3 dark possibilities:
1. ignorance
2. deception
3. ignorance & deception
Nick Stokes:![\sqrt{a^2 + b^2}\cos[t + \arctan(b/a)]](https://s0.wp.com/latex.php?latex=%5Csqrt%7Ba%5E2+%2B+b%5E2%7D%5Ccos%5Bt+%2B+%5Carctan%28b%2Fa%29%5D&bg=ffffff&fg=000&s=0&c=20201002)
, exactly what does the fact that the resultant standard error is a certain multiple of the determined sinusoidal-component magnitude tell us about how “significant” our result is if the signal consists exclusively of, say, five equal-amplitude frequency components?
Thanks for the response. It was helpful, but my ignorance of statistics compels me to take some time, which I probably won’t have any more of today, to explain to myself in simple terms what it really means .
That is, if we eliminate the independence issue by, say, correlating with a single sequence
I’m not arguing here; I’m just thinking on “paper.”
Fourier transform is a fancy name for a Fourier series which was in the A level maths syllabus and I subsequently used them in heat transfer calculations involving time temperature distributioin..
Joe Born says: May 5, 2014 at 4:57 am
I think you’re right. If you fit one component, the residuals contain the effect of the noise plus the other components. But the other components would be included in the estimate of variance as if they were noise.
If you fitted the five components jointly, you’d get a much lower estimate. But in general you don’t have that information.
chemengrls says: Fourier transform is a fancy name for a Fourier series
No Fourier transform is a fancy name for a Fourier transform . It’s a mathematical operation which converts from the time domain representation to frequency domain representation .
Fourier series is a fancy name for a Fourier series, a theoretically infinite series of harmonic terms, with individual coefficients.
Go back a read your A level maths, and pay attention this time.
You might also pick up an old copy of the Dover Books Classic “The Measurement of Power Spectra”, by Blackman & Tukey.
It’s been my companion for over 50 years.
Greg says:
“Could you more specific about what you consider to be “proper selection” and “proper constraints”?”
Those were of course weasel way of saying it’s probably going to work but I am not sure how yet.
The problem is that optimization algorithms could produce unexpected results and should be approached carefully. There is always a danger of them finding wrong solution.
That said, after thinking about it a bit more, I think the way Willis does it will work fines is.
To be sure, I’d have to do some testing on it (and by “I” I mean someone else – as I do not, unfortunately, have time to do it myself 🙁 )
But right now, my thinking is that SFT would not assume periodicity of the dataset and would be completely immune to edge issues add that as one more advantage of this method vs regular FFT or DFTs
Joe Born says:
May 4, 2014 at 3:45 pm
I’m using the lm() function in R. The R linear model function reports the statistics of the fit as a p-value. When I say a 90% threshold, I mean a p-value of less than 0.1 for the linear model fit after correcting for autocorrelation. Sorry for the confusion.
w.
Well, “pressure broadening” and its “pathological” distribution could be killing the Cheshire cats, from time to time.
I’m late to this party and don’t have time to slog through all the other comments. So I’ll apologize if anyone else has said this before.
1: You might want to to try visualizing your data as a “waterfall”. Where you take smaller overlapping lengths of time and watch the peaks move in time.
2: This looks very similar to Modal analysis from mechanical engineering. The only problem with this is the “Unit Under Test” and input forcing function is not constant. So it is much harder to compare different locations with a linear transform function. (matching peaks).
Leo Smith says:
May 5, 2014 at 2:35 am
Thanks, Leo. What do you mean by using an interpolation to “caret” an “analogue” time series? How do you decide the frequency of the sinusoidal interpolation function? How do you decide if the sinusoidal curve should go up or down? What does a sinusoidal interpretation of a square wave that is missing 40% of its data look like? And in the square wave case, shouldn’t we be using linear interpolation? These questions and more …
A google search on either “caret” or “sine interpolation” finds little of assistance. I do find things on sinc interpolation [ sine(x)/x ], however, such as this:
Is sinc interpolation what you were referring to? Or did you actually mean sine interpolation?
Also, that paper says that when the data is “clumped”, the sinc interpolation is unsatisfactory. And unfortunately, climate data is often “clumped”, with long stretches missing entirely.
Let me offer up an example. The data for this example is here as a small CSV file. It’s 180 years of an 18-year squarewave, plus the same square wave with about a third of the data points missing. I’ve taken out both long stretches of years, as well as deleting a number of individual months of data. This is the data.
Top panel is the square wave, middle panel is the square wave with the data chopped out of it, and the bottom panel shows what data was removed, since the scale is too large to see a missing data point in the middle panel.
My problem is that your method involves an “interpolation” of that middle panel … and I’m having a hard time conceptualizing what that “interpolation” might look like.
In any case, here’s what the SFT does with the two datasets:
Upper panel is the slow Fourier transform periodogram of an 18-year square wave with the same number of data points as in the lower panel. I wanted to compare two datasets of the same length.
Lower panel is the periodogram of the chopped-up square wave. Given that a third of the data missing, I’d say that the slow Fourier transform is doing a respectable job …
I’d be curious to see what your sine [sinc?] interpolation would do with that data.
Never claimed it did … I pointed out what I said were advantages over the FFT. You seem to be agreeing that at least some of those advantages are real. You say the SFT is equivalent to using an FFT on a sine interpolation of the data. Perhaps so.
I pointed out that I didn’t think that I’d invented the method. You are agreeing with that also, although it remains to be seen if your proposed equivalence is correct. I always get nervous when someone says “a cursory inspection” … for example, my method certainly has the advantage of being easier to explain.
I’d be interested in your response to what I said in bold was the largest advantage of the SFT, which was that the SFT method provides an accurate measure of the statistical significance of each cycle length, based not on theory but on the actual data involved, including the size and nature of the gaps in the data.
Thanks for that, Leo. Perhaps you could at least run your analysis on the square wave CSV data I linked to above and show your results.
Or perhaps you could point me to a simple explanation of the method of “sine interpolation”. I’m having a hard time believing that my output will be the same as the output of the sine interpolation, inter alia in the matter of the statistics of the results. My method, since it doesn’t interpolate, will have a much lower “N” (number of data points) than any analysis of interpolated data. How does sine interpolation handle that issue?
My method gives me a standard and well understood way to calculate the actual statistical significance of the various individual frequencies. So for example, if the missing data has a pattern, say every 10th data point is missing, then that will affect some frequencies more than others. My method allows the determination of exactly how the pattern of the missing data affects the results.
Does sine interpolation of a “resampled” series do the same? I don’t know.
Also, suppose we have a totally irregular array of data points, which have no common sampling frequency smaller than their least common multiplier. Using my method, I do the analysis on those actual points.
Using your method, as I understand it, not one of those actual points will be sampled. Instead, a series of arbitrary regularly spaced points will be created using sine interpolation, and then the calculated values at those arbitrary points will be used for the FFT.
Is my understanding correct? And if so … well, I’m not seeing how doing an FFT on an entire string of interpolated points will give exactly the same answer as me using the SFT on the actual points. It just seems like the extra step of the sine [sinc?] interpolation would have to increase the error of the results. Might not happen, I’ve been surprised before … but I’m not understanding how it would work.
I know your time is short, but any assistance would be appreciated. And my thanks for further discussion of what it is that the SFT might be doing. You say the SFT is exactly equivalent to a standard Fourier transform of the sine (sinc?) interpolation of the data. I’m very interested to see if that is the case.
Regards,
w.
Nick Stokes says:
May 5, 2014 at 3:55 am
Thanks, Nick. To avoid the extra calculation, I just took the “fitted values” from the lm() function.
Then I used lm() again comparing the fitted values themselves to the data. This gives me the statistics for the final combination of the sine and cosine terms.
I think there is a more accurate way, however, that will help avoid false positives. More on that if it works out.
Since I said twice in the head post that I had adjusted for autocorrelation, I’m not sure what this means.
Your contributions appreciated,
w.
Blue Sky says:
May 4, 2014 at 4:24 pm
Thanks, Blue. From the head post:
I’ve been wrong more than once. You’ll have to point out which time you are referring to.
Regards,
w.