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.



Peter Sable said in part May 30, 2014 at 8:07 pm
“… O&S is a terrible textbook. …”
Yea – I think 97% of students would agree. I first learned DSP using a mimeographed edition probably about 1973, and have seen three published editions, and worked with students who were required to have it. ( When I taught a course myself, I never used any textbook.) I have heard it was always (or at least originally) intended to be a summary collection rather than a text, which makes sense for the earliest editions – there weren’t many courses. Most students did not like it at all. Later they were more apt to agree that if you had learned the material once, it was a reference they needed on their shelves. But true – you can’t just send someone right to page 733 (of 870 pages!) and expect anything useful to occur.
Nick Stokes:
I appreciate your clarification regarding the bandpass characteristic of LC
circuits, as seen in its freq. response. But we’re back to the misleading
tuning-fork analogy. To be sure, it will respond to any finite-duration
excitation and “ring” thereafter at its resonant frequency. However, that
is far removed from straightforward heterodyning of a PRE-CHOSEN sinusoid
with a data record–especially when that record LACKS any finite-amplitude
sinusoidal components.
That is the usual case in climate “anomaly” series, wherein diurnal and
annual cycles have been removed. Being the product of chaotic physical
processes governed by hyperbolic PDEs, these are oscillatory–but NOT
periodic. The power density spectrum is continuous, instead of
concentrated in discrete lines. Some number will, of course, be produced
with ANY chosen period, but its essential MEANING is not the same. In
practice, sinusoidal regression remains a curve-fitting method suitable
ONLY for periodic signals with KNOWN periods, such as the tidal components.
Bernie:
My conditional offer was not meant to suggest that you stand destitute,
although if cutbacks continue some co-workers may have to purchase tin
cups. 😉
What I liked was the demo of the Gibbs effect you presented here:
http://electronotes.netfirms.com/FFTInterpSquare.jpg
which should remind all practitioners of the pitfalls of truncation in
finite-length DFTs. The Rayleigh criterion, a resolution concept borrowed
from physics, effectively specifies that sinusoidal orthogonality is
achieved when the first zero of the spectral window centered on a frequency
band coincides with the peak response of the adjoining band. Thus sinusoidal
regression in keeping with this criterion will simply replicate what is
obtained by FFTs.
Peter Sable:
While we all agree that O&S may be far from the ideal teaching text, Chapter
11 does cover the critical differences in analyzing periodic vs. random
signals with many examples and is advertised as suitable in part for
undergrads. Bendat & Piersol may be a better read for some minds.
Also it includes a chpater on cepstrum analysis, useful for sonar and radar work.
My objection to plotting power density curves v. period is the destruction
of any sense of the AREA under the density curve, which in many cases with
short records is more important than the exact location of a peak. On the
other hand, periodograms (suitable primarily for analyzing periodic
records, not random series), should be presented as discrete lines, rather
than curves.
1sky1 says: May 31, 2014 at 2:21 pm
“However, that is far removed from straightforward heterodyning of a PRE-CHOSEN sinusoid with a data record”
No, it’s the same. Look again at that time-varying integral solution, written with general L and C, and integrated by parts to get back to V(t):
L*I(t) = cos(ωt)∫ cos(ωτ) V(τ) dτ – sin(ωt)∫ sin(ωτ) V(τ) dτ
again int over all real, but with V(t) of finite duration. ω=1/sqrt(LC)
Each coefficient is exactly the heterodyning integral – multiply the signal by the prechosen (tuning fork) test sinusoid trig(ωτ).
“again int over all real” – wrong, the upper limit is t. And the lower limit might as well be 0, since V(t)=0 outside a finite range,
“The power density spectrum is continuous, instead of concentrated in discrete lines.”
Which is exactly why seeking more plotted resolution that is provided by DFT is desirable.
There have been several comments here to the effect that there is no spectral information available in the data between the Fourier periods ( frequencies) and that zero padding does not provide any more information.
This is a false conclusion that is based on the failure to realise the different between a spectral analysis and a TRANFORMATION.
FT in all its variants provides an equivalent F.Dom representation of the data that is fully equivalent to the T.dom representation and allows it to be reconstructed by iFT This is not a spectral analysis, it is a mathematically equivalent transformation. You can then perform a spectral analysis by looking at the FT but for medium and long periods it will be a crude approximation. That makes it a poor analysis.
The Fourier frequencies are all harmonics of the length of the data sample. It is patently false to suggest that there is no frequency information in the data for other values. The only true limitation is the Nyquist cut-off.
If you do the FFT of 200y of data of a pure 140y sine the FT will not contain a 140y period. But it will contain N ‘incorrect’ frequency components that enable a reconstruction of the correct periodic form over the sample range.
Now if you zero pad to 280y you will find the FT is a lot simpler 😉 and the amplitude stronger than the largest amplitude in N previous frequencies.
This is how progressively padding, by shifting the fixed FT frequencies, does allow a more precise detection of where peaks lie. This concerns the resolution of the analysis. It will never increase the resolution of the data. Many seem to have similar problems differentiating the two.
Of course adding zeros can not improve the resolution of the data. It can improve the resolution of the analysis.
oops: more resolution _than_ is provided by DFT
Re the universe as numbers: to be able to ascribe number and dimension to any object or phenomenon’s properties, it is only necessary that it be consistent from instance to instance. So discovering numbers is discovering consistency. The consistency is not the thing, per se, but it is essential to identification, otherwise it would not be recognizable in the first place. So calling the universe one of numbers is a kind of tautology: “consistencies exist”.
Of course adding zeros can not improve the resolution of the data. It can improve the resolution of the analysis.
Greg: References?
I find there is a dearth of references on how good low frequency information is with an FFT. I would appreciate some references so I can do further reading.
I tried padding with zeros. Yes, I got a higher resolution plot in the low frequency (high period) section of the graph, but the energy was smeared across just as large a width as before – just with less jagged lines.
To show this I tried two sine waves one year in period apart, 32, and 33 years. Can’t tell the difference.
Further messing around shows that at 16 years, I can discern 2/3 year of a difference in overlapping peaks. At 32 years, it’s about 1/4 that. At 64 year’s it’s about 1/16th the resolution as at 16 years. Not unexpected. (note, was actually about 1/18th, so no numerical guessing, should try to derive an actual function symbolically… but I have a day job).
“I tried padding with zeros. Yes, I got a higher resolution plot in the low frequency (high period) section of the graph, but the energy was smeared across just as large a width as before ”
In fact you got higher plotting frequency resolution _throughout_ the spectrum. This was most obvious in long periods when plotting period, not frequency.
Of course it does not change energy distribution, it just shows it in more detail. That is what I was saying about improving the resolution of the analysis.
“Can’t tell the difference.” They have a frequency separation of 1/32-1/33=0.000947 per year, sampled once per month or whatever you did. That is what is meant by resolution of the data that you can’t improve unless you have more samples in the same period, not just more samples.
Sorry, I don’t have any links to hand and I too have a ‘day job’ so I’m not going to do your googling.
“I find there is a dearth of references on how good low frequency information is with an FFT. I would appreciate some references so I can do further reading.”
There’s a dearth because you’re asking the wrong question. The _frequency_ resolution is the same across the spectrum.
Nick Stokes:
Your analogy with a tuning fork or LC circuit subject to arbitrary TRANSIENT
excitation continues to mislead you. That’s a RESONANT physical system,
which acquires and stores energy from TIME-domain input, concentrating into
periodic output at the resonance frequency.
By contrast, the record-length dependent spectral window is a mathematical
construct, whose NON-RESONANT bandpass frequency response is entirely
different in form and operates in the FREQUENCY domain, simply integrating
the spectral content of NON-transient signals in a variance-preserving
manner. The spectral window never disappears and you cannot get a pure
sinusoid unless the input signal consists of spectral lines. Period!
Peter Sable:
You are entirely correct! Zero-padding adds no signal INFORMATION; it
merely increases the COMPUTATIONAL resolution in the frequency domain,
which amateurs mistake for ANALYTIC improvement. This can be seen by
considering the zero effect that padding has upon the sample acf. The only
way that the resolution-limitation of the spectral window can be overcome is
by acquiring more data. By supposing a pre-set period, sinusoidal
regression produces a number, but its significance is dependent entirely
upon the validity of that supposition.
An elementary demonstration of this can made by adding all the sinusoids
that Willis finds for the SSN data. The result is horribly out of scale,
because of non-othogonality of decomposition. And if the largest-amplitude
“periodicity” alone is compared to the data, the mean-square error is huge,
rendering that periodicity quite useless in signal identification or
prediction. The claim of “advantages” is demonstrably pathetic.
Why is an increase in computational resolution not, at least potentially, an analytical improvement?
Increasing the computational resolution can help to correctly identify a period that does not have the good fortune to fall exactly on one of the Fourier harmonics, as I explained here:
http://wattsupwiththat.com/2014/05/26/well-color-me-gobsmacked/#comment-1651187
” And if the largest-amplitude “periodicity” alone is compared to the data, the mean-square error is huge, rendering that periodicity quite useless ”
Quelle surprise! You take one component out of a spectrum and it has a “huge” residual. To follow that logic the whole process is “quite useless”.
However, if you take the largest component from a swept analysis found by progressively increased padding, that finds where the largest signal actually is in the data, rather that hoping it lies near a Fourier harmonic, you will in general get a smaller residual.
Again to take my example, the largest unpadded F cmpt would have a “huge” residual, the 280y padded analysis would have zero residual. Where the input is not itself a pure harmonic function, the result will be somewhere between these two extremes.
“The _frequency_ resolution is the same across the spectrum.”
Ah, and the derivative shown as as percentage error changes with the frequency, which is what I was intuiting but not explaining well.. Okay, I have enough clues to figure out the equation for that.
thanks, I’m catching up after a 15 year hiatus in signal processing. Between the day job and the grayer cells it takes a bit…