Guest Post by Willis Eschenbach
I’ve developed a curious kind of algorithm that I’ve humorously called “slow Fourier analysis”, but which might be better described as a method for direct spectral analysis. The basic concept is quite simple. I fit a sine wave with a certain cycle length, say 19 months, to the dataset of interest, and I note the peak-to-peak amplitude of this best-fit sine wave. I repeat this process for the various cycle lengths of interest, and put that all together as a periodogram. The periodogram shows how large of a sine wave is the best fit for the data for each of the cycle lengths. This method has worked a treat, but yesterday, this odd method of mine handed me a big surprise.
Now, I have been doing these periodograms utilizing only the time periods inherent in the data. So if the data was monthly, I’ve looked at all possible integral lengths of monthly cycles, starting with two month, three month, four month, five month cycles, and so on up to cycles with periods as long as a third the number of months in the data.
And on the other hand, if the data was yearly, I’ve looked at all possible integral lengths of yearly cycles, starting with two year cycles, three year cycles, four years, and so on up to a third of the number of years in the data.
(I mistrust reports of the longer cycles, and any cycle with a period longer than a third of the dataset is not really worth calculating—it will fool you every time.)
So I was going along in blissful ignorance when an alert WUWT reader, Greg Goodman, pointed out that there was no reason to restrict myself to integral periods. He noted that I could subdivide the results and gain greater resolution. I was dismissive of his idea. I said that I thought you could do it, and that he was right and it was an interesting idea, but I said I thought you wouldn’t gain any real resolution by doing that. Hang on, let me get my exact words and the graphic I used in that thread to illustrate my point … here’s what I said:
In any case, I’m still not convinced that the procedure will give any real increase in information. Here’s the difference, for example, between monthly and annual data:
I think that if I increase the resolution on the sunspot periodogram while still using the annual data, it won’t look anything like the real monthly data … I’ll report back on that one.
Well … this is the promised “report back”. As you can see, both the monthly and the annual data show the peak at 11 years. The monthly data shows more detail, as it resolves the individual peaks, at 10, 11, and ~ 11.9 years. I thought that the test I outlined would be a good one. I’d see what happened when I sampled annual data at monthly intervals. So here’s what I did. I took the SIDC monthly sunspot data. For the annual data, I didn’t use the SIDC annual data. Instead, I averaged the SIDC monthly data to give me annual data that I knew was a true annual average of that particular monthly dataset. Figure 2 shows the monthly data in red, and the annually averaged monthly data in black.
Figure 2. Monthly SIDC sunspot data (red), and the same monthly SIDC data annually averaged (black).
So to be clear, the data shown as the black line is the annual average of monthly data shown as the red line
Next, I modified my periodogram-generating computer function to allow for fractional time sampling of the data.
Finally, I calculated three periodograms. One is the periodogram of the annually averaged data, shown in yellow/black below. The second is the periodogram of the underlying monthly data, shown in blue.
And finally, I calculated the third periodogram (shown in cyan) using the annually averaged data, but sampled on a monthly basis … Figures 3 and 4 show those results. I must confess, when I saw them my jaw hit the floor.
Figures 3 and 4. Upper panel shows the periodogram of annual sunspot data (yellow/black), monthly sunspot data (blue), and monthly-sampled annual data (cyan). Lower panel is a closeup of the upper panel, showing the period from seven to twenty-five years.
You can see why my jaw hit the floor. The blue line (periodogram of actual monthly data) is almost identical to the cyan line (periodogram of annual data sampled at monthly intervals). Who knew?
It turns out that contrary to my expectation, the information about the strength of the monthly cycles is NOT lost when the individual monthly values are subsumed into annual averages. Somehow, that information still exists and can be retrieved. When I take the annual data and analyze it at monthly frequencies, the data is still there. Not exactly, to be sure, as you can see the process is not perfect or precise.
But it is astoundingly accurate. I absolutely didn’t expect that.
Now, I suppose some signal analysis guy will pop up and say he knew that all along and abuse me for my colossal ignorance … but I taught myself this stuff, and all I can say is, it sure was a surprise to me.
Truly, I don’t understand this result. For many, many years I’ve thought that when you average data in chunks, like say averaging monthly data into yearly data, that all of the monthly information was gone. Lost. Kaput. Irretrievable.
But that doesn’t seem to be the case at all. It seems that very little of the power spectrum information was lost at all as a result of the annual averaging process.
Naturally, of course, this brings up the next question—is this procedure invertible? That is to say, given a periodogram such as the one above, can I run the process in reverse? Can I start from the annual data, calculate the monthly periodogram from that annual data, and then invert the periodogram to give me back the monthly data? That would be really amazing … but I’m pessimistic.
My guess is no, the periodogram can’t be inverted to me give back the monthly data … but given how poor my last pathetic guess (at the head of this post) was, I’ll certainly give it a try. Any assistance gladly accepted.
Like I said above … once again all I can say is, I’ll report back on that one.
Regards to all,
w.
Acknowledgement: My thanks to WUWT reader Greg Goodman for the suggestion to investigate fractional time periods.
For Clarity: If you disagree with something someone says, please quote the exact words you disagree with. It avoids all kinds of misunderstandings.
Data and Code: I’ve collected the R code, the R functions, and the data into a small (22KB) zipped folder called “Oversampling Folder“. It includes the monthly and annual sunspot data as CSV files. If you change your R workspace to the folder, it should be turnkey.
For the original data:
SIDC Monthly Sunspot Data Note that following the advice of Leif Svalgaard, I have increased all of the pre-1947 SIDC sunspot values by 20%, to correct for the change in counting methods at that time. The change is immaterial for this analysis.
Changes In The R Function: The function to generate the periodogram, named “sineptop” (for sine peak-to-peak), was previously called like this:
sineptop(annual_data, frequency = 1)
or like this:
sineptop(monthly_data,frequency = 12)
The “frequency” variable identifies the data as having that many periods per year.
It can still be called like that. In addition, however, the new sineptop syntax for e.g. monthly sampling an annual dataset looks like this:
sineptop(annual_data, frequency = 1, by = 1/12)
If the “by” variable is not specified it is assumed to be 1, so you don’t need to specify it, and the original syntax still works.
Discover more from Watts Up With That?
Subscribe to get the latest posts sent to your email.
Willis, I think Fred Haynie did something very similar a year or two ago. He fitted sinewaves to the past temp record to show the possible future. Not sure though. Everythings a bit hazy at my age.
Bernie Hutchins says: “Inverting a Moving Average Operation”
http://electronotes.netfirms.com/AN407.pdf
Thanks. I will read that in detail later. I’m currently working on an optical problem which is equivalent to two running means burring an optical spectrogram. I’m currently deconvoluting with sinc squared (the two being nominally equal). Deconvolution tends to blow up the noise and it’s already rather a noisy spectrum. I’ll see whether I can find a way improve it in your paper.
Do you have anything more advanced on that technique?
Willis “But although the annual averaging is a boxcar average, is that convolution as well?”
Very much so. A running mean is one of the most simple convolutions. It is equivalent to moving your hand when taking a photo. It blurs the result. An annual average is exactly the same as a running mean that is them sampled once per year.
Of general interest, the same thing applies to de-blurring and refocusing photos. A good explanation of blurring and de-blurring by convolution is given here.
http://www.fmwconcepts.com/imagemagick/cameradeblur/index.php
BTW one trick I’ve used in the past is convert time series data to 1 pixel high grayscale image and pass it to Imagemagick for “deblurring” or other FFT processing. It’s a useful cross-check on results.
CC Squid says:
May 26, 2014 at 10:05 pm
So you have answered your question on 11 year sun spot cycles? ((2010-1750)/25= 10.4 Does this means that the Schwabian cycles are a fact? I love this blog and get great pleasure from your posts. Congratulations!
====
Ah thank you. That probably explains the poorly defined bump on the side of the central peak I noticed in Willis’ early plots, that I don’t find by more conventional processing. This is what windowing or “taper” functions are for. I should have thought of that straight way.
That was one reason why I was suggesting better plotting resolution would be good, since I could not determine clearly where that peak was.
You can see this bump in figure 4
Greg Locock says:
May 26, 2014 at 7:39 pm
You cannot use this method. The reason is that the resolution of a Fourier analysis is limited to B=1/T where T is the total length of the sample. Therefore in a 1000 year sample your frequency lines are spaced at 1/1000 y-1. The only valid frequency lines are for integer numbers of cycles in the total sample.
I suggest you read up on Fourier.
………………
Hi Greg, it is useless. Willis is not reader. He is writer.
Willis, you can recover *exactly* the monthly data if you deconvolve it. Because you know *exactly* the filter you convolved it with you can exactly decomvolve it and get your original data back. That is, if you keep the density of data. If you throw out 1/12th of the dat then you lose some information.
Willis, I do have a challenge for you, though. You seem to be utterly convinced this method can give you data that the FT can’t (it cannot). Go back and do EXACTLY what you just did with the DFT and I think you’ll get similar results.
Willis,
You said:
[quote]
That makes it sound like you could “un-average” the data, which still seems like magic to me … but I was born yesterday, what do I know? I was thinking maybe my method wasn’t invertible because the decomposition is not orthogonal … any clues on how to do the inversion with a non-orthogonal decomposition?
[/quote]
I would agree with you: you can’t just add up the components after your projection, since the signals (base vectors) you’re projecting on are not orthogonal. However, in practise, this is also not the case for Fourier Transforms. In FT, the components are only orthogonal if you have an infinite time window; so you get spurious frequencies there as well in practical applications.
A mathematically stricter solution, I guess, is a kind of principal component analysis; i.e. you find the strongest sine wave (its amplitude and phase), subtract it, find the second strongest signal of what’s left, subtract it, etc. You stop when there is nothing left but noise. All the components you find through this procedure can be added back together, even though they are not orthogonal.
Frank
***
tamino says:
May 26, 2014 at 7:57 pm
***
I’m also “gobsmacked”. A reasonable & readable reply from tamino.
wsbriggs: “It seems the mouth breathers aren’t as prevalent in the responses these days, or maybe I’m missing something.”
yes it seems that the volume of political yabooing has dropped considerably. They must have got bored or worked it out of their systems.
If you convert the years to months, an 11 year cycle requires 132 sample points, whereas a year of data is 12 points. You would have to examine your graph from 2 months to 24 months to see any loss of resolution. Of course you couldn’t do this because the yearly data does not contain this information, having a cut-off (highest resolvable frequency) of 24 months. That is why your periodogram starts at 24 months (2 years). Any difference would be to the left of your starting axis.
Greg Locock says: . The only valid frequency lines are for integer numbers of cycles in the total sample. I suggest you read up on Fourier.
“Valid” for what ?!
So If you have 100 years of data with a pure 30 year cycle , you would not regard fitting 30 y as “valid”. The only “valid” frequencies in the data would be 50,25,12.5 …..
Merrick says: Willis, you can recover *exactly* the monthly data if you deconvolve it.
Have you ever actually done a deconvolution?
An excellent reply by Tamino, instructive and informative, as well as courteous. Imho he nails it, with extra insight on the applicability of this method.
Dare we hope that this cooperative spirit is an indication of more harmonious interactions in the future?
Anually-averaged data is just a bad low-pass filter. It attenuates high frequencies, but the frequencies you looking are low.
Filters have an effect on the time domain, but also you can consider them in the frequency domain. averaging -> attenuating high freqs
Willis, I don’t see how an ≈ 11 year cycle could have blossomed from the data unless it were there already. Give the data to Mikey…he’ll turn it into a hokey stick.
Now there’s some magic.
cn
Willis: “Thanks, Michael. I’m not sure I understand that. For me, “convolution” means sweeping a filter along the length of a dataset, with the value of the filter applied to each point in sequence. But although the annual averaging is a boxcar average, is that convolution as well? ”
What Michael D was trying to say is multiplication in time domain = convolution in freq domain and conversely multiplication in FD = convolution in TD.
Both multiplication and convolution can be done in any order.
It’s quite a neat result once you know what it all means.
Boxcar filter is a convolution. The annual sampling is not. The sampling is multiplication by a train of impulses (called the Shah function). In effect this convolves the frequency spectrum by the FT of Shah, which it happens, is also a Shah but with frequency spacing. What this does when you convolve the real spectrum is superimpose lots of copies of the real FT of the signal side by side. This why you need to anti-alias filter BEFORE re-sampling. Once you’ve done it, the mix is done and you can’t filter it back out.
The problem with averaging is that the RM part , which is implicit has it’s zero at a period twice as long as the Nyquist freq, required for anti-alias filtering. Not only than but it does not cut off cleanly below that but leaks an INVERTED signal which will also alias.
Now if all you have is random noise you can hope it all blows out but if there is some periodic content between nyquist and the RM or in the leaky inverting lobe, it will alias.
This is why averaging ( eg the ubiquitous monthly and annual means in climate) are not good practice.
Sometimes it’s small and does not matter but when you see a “periodic” in the result you never know.
I like Willis’ preemptive requirement.
“For Clarity: If you disagree with something someone says, please quote the exact words you disagree with. It avoids all kinds of misunderstandings.”
Of course it can avoid more than “misunderstandings”.
Hopefully it will also prevent the purposefully mendacious use of,
“So what you are really saying is…” method of distortion and straw man diversion.
The frequency sampling of the FFT is sample_rate/fft_size.
If you want a finer sampling frequency you can zero-pad the FFT, thereby making it longer.
The DFT (slow fft) can do this, but the calculations are proportional to N^2 vs. N*log(N).
By the way of example, say your data is 372 points long. left justify this data in a buffer that is 1024-points long and do the FFT. Your frequency sampling will be finer by the ratio of 1024/372.
In Matlab this would be something like this:
rbuf=zeros(1024,1);
rbuf(1:372)=mydata;
freq_resp=abs(fft(rbuf));
You can zero-pad out as far as you want. Just take care of the normalization (scaling) factor.
For some FFT’s you will have to correct the amplitude by the zero-padding ratio.
hi Willis,
What Greg Goodman just wrote makes sense to me. FWIW, I’m largely self-taught in filtering, but had lots of math education with grad degrees in physics. Other commentors are thus able to speak with more authority on this.
Imagine doing a full convolution using the annual-average “boxcar” filter: that would mean sweeping the boxcar month by month, taking the average from Jan01-Dec01, then from Feb01-Jan02, then Mar01-Feb03, etc. I think I called that (above) a “sliding boxcar average”. If you sub-sample that convolution data by saving only the points representing Jan-Dec (thus tossing out 11 of the 12 data points for each year) then that is (I think) a standard annual average, which is what you were working with, so that’s why I was thinking of your data as a “convolution” with the boxcar filter. It would more accurate to call it a subsampled boxcar convolution I suppose.
As a thought experiment, imagine “reconstructing” the sliding boxcar average from the subsampled boxcar. The best you could do would be to do a linear interpolation between the Jan01-Dec01 value and the Jan02-Dec02, etc. to estimate the intermediate months (Feb01-Jan02, Mar01-Feb02, etc.) values. The difference between the original sliding boxcar and the reconstructed boxcar represents the information that has been lost by subsampling. It is in the nature of a boxcar average that the difference will typically be about 1/12 of the signal in this case.
If you use a non-subsampled (i.e. “dense”) boxcar average, followed by your Slow Fourier, you probably haven’t lost much information (though the boxcar average is blind to certain frequencies – there are zeros in its frequency response). To reconstruct the original: arrange your Slow Fourier frequencies to be spaced evenly (i.e: 0, 1, 2, 3,…N) [cycles/(N years)] then do a Fourier transform on that set (if you are not working with complex numbers, the inverse Fourier is identical to the forward Fourier). With some re-scaling (maybe just 1/N) you should see the original pattern come back. You may need to create a “cosine” version of Slow Fourier, repeat the above using cosine Slow Fourier, and then add the sine and cosine results.
i don’t buy the “one year of college calculus” excuse – you are obviously a math whiz and you are close to doing “standard” Fourier domain stuff, you’d probably enjoy delving into that some more. And that would connect you to a bunch of theorems that are helpful.
It is quite incorrect (as some are suggesting) to believe that frequencies are restricted to integer multiples, except as a particular analysis tool (such as FFT) might require this, and the “uncertainty principle” (limit of resolution) likewise depends on the analysis tool. “Rules” (and limitations) that apply to one tool may not apply to another tool, which has its own limitations.
Consider the three time samples:
[0 0.01745240643728 0.03489949670250]
These might look like a straight line, but the third one tips down ever so slightly. So you might “suppose” that they are the first few samples of a very low frequency sinewave. (They are in fact the sines of 0, 1, and 2 degrees.) If you try to fit them to [Sin(0) Sin(theta) Sin(2 theta) ] (i.e., solve for theta) you discover there are enough equations, but they are non-linear.
So try the FFT. The FFT lies! It says that there is a DC term and a frequency of 1/3 (red curve here). It never heard of the correct frequency 1/360 (blue curve).
http://electronotes.netfirms.com/Pdemo.jpg
Prony’s method (1795) allows us to solve for the exactly correct frequency in terms of the parameters of a network that could have produced the samples. Any bandlimited frequency is possible. For example, sqrt(2)/360, (or 1/3 for that matter).
Limitations on Prony? You have to guess the order correctly, and it doesn’t like noise, but it has no preferred frequencies. Each analysis method should be considered individually.
As I re-read what I just wrote, I realize that I’ve been assuming that your “best fit” of the sine wave was akin to a convolution of the sine and cosine waves with the signal. But there are many ways to do best fits, so that is probably a bad assumption. Also if you are fitting the phase of the sine wave, then you are fitting both a cosine and a sine wave, so you don’t need a “cosine Slow Fourier” you just need to split the sine wave into cos and sin waves with different amplitudes and zero phase:
sin(wt + phi) = cos(phi)*sin(wt) + sin(phi)*cos(wt)
Good Day Willis,
You seem to be falling into the “average kid trap”. The average 5th grader is 5 feet tall and weighs 110 pounds. So, we put 30 desk of the appropriate size in a classroom and expect it to fit the 30 kids that show up. When it is a bad fit, blame the kids, the parents, the budget, the voters…
You are looking for an 11 year sun/climate correlation, because the average solar cycle is 11 years. Solar cycles come in a variety of lengths and intensities. (And, probably vary in other ways we don’t know about or measure yet, but, I digress) If I was going to look for a correlation between solar cycles and climate I would look at the type of cycle.
I would start by categorizing the cycles into 3 groups: shorter than average, average and longer than average. Do groups of shorter cycles lead to warming? Just eyeballing the data, it looks like it might. This also requires looking at how we measure the length of a cycle (since they overlap). Does it matter how we choose the dates a cycle begins and ends? Does it matter how long the overlap is?
I would also look at the number of sunspots smoothed over 11 years. In case I’m not using the term correctly: sunspots for 2005 = the annual spots 2003 through 2013 divided by 11. Just eyeballing it, sustained low sunspot activity seems to link to cold or cooling periods.
I would also look at where the earth was in relationship to the sun at solar maximum. Does it matter if the solar max happens during summer or winter solstice? (Given that the southern hemisphere is covered by mostly ocean and the northern by mostly land.)
I would also look at the length and intensity of the solar minimums.
PS. Whatever drives global climate is certainly more than just the sun, but when a unified theory of climate is finally established, I do expect it will begin with the sun.
*I’m a civilian and not a scientist, if common sense has led me astray, I would appreciate correction.
this post should have a warning … no science follows but statistics in spades …
“Thanks, Michael. I’m not sure I understand that. For me, “convolution” means sweeping a filter along the length of a dataset, with the value of the filter applied to each point in sequence. But although the annual averaging is a boxcar average, is that convolution as well?”
You actually have to apply two steps – first convolve with the boxcar filter, then multiply by a function that contains a unit impulse (delta function) in the middle of every year. This will sample the running boxcar filter annually.
As noted above, convolution and multiplication are Fourier transforms of one another, so the effect on the frequency spectrum is to first multiply by the transform of the one-year boxcar (a sinc function) and then convolve with the transform of the sampling function (which is another sampling function with inverse frequency). The Sinc function is broad, and only cuts off the highest frequencies (although it will distort lower frequencies because it’s not constant), and the transformed sample function has widely separated peaks, which pastes copies of the entire spectrum at regular intervals in frequency. Both of these only mess up the high frequency end, which since you’re only plotting periods 2 years and above you have cut off.
Both convolutions and multiplications can be inverted *so long as they are not zero*. (You just look at it in whichever domain involves multiplication by a fixed function, and divide by it instead.) The Sinc function crosses the axis (is zero) at regular intervals, and the sampling function is zero along most of its length. Because of the latter, you can only invert if the bandwidth is less than the Nyquist frequency, and because of the former, you can’t reconstruct frequency components with periods that are exact integer fractions of a year, although you can theoretically reconstruct the rest of the signal.
Because your signal is of only finite duration, which is equivalent to taking an infinitely long signal and multiplying by a boxcar function the length of your data, this has the effect of convolving your spectrum with a correspondingly narrow Sinc function, blurring the fine details. This blurring spreads the nulls out somewhat, and means it’s difficult to reconstruct frequencies very close to the nulls, as well. And if there’s any measurement noise or rounding, this gets magnified around the nulls by deconvolution, spoiling your result.