Part I: Homemade data:
This exercise uses a Butterworth filter
to look at the relationship between high frequency and low frequency
variations. The program filt.d on the software page might help you get off to a faster start, maybe not.
1.) Data: I have put one 4096 file called x552filt in my public
directory on eos.atmos.washington.edu/ftp/pub/dennis/552_Class/.
2.) Spectral analysis: you will find that the spectrum is red, but there is a peak in the variance around 5 days, by doing spectral analysis using the mspec.m program, suitably adapted to these data sets.
The idea is to separate the high frequency variance around 5 days period from the much larger amount of low frequency variance. You can see from the power spectral plot that there is more variance at lower frequencies. This is a little more obvious if you plot the spectrum linearly, rather than logarithmically.
You will note a gap between the low frequency variance and the
bump at 5 days. This is unusually good fortune, as usually you
will have a lot of variance at all frequencies and it is not obvious
where to cut. Now we want to divide the signal into a low-pass
and high-pass component using a Butterworth filter.
3.) Butterworth filter construction: You want to design a butterworth filter that makes a cut between the high and low frequency variance, somewhere below a frequency of 0.2 cycles per time interval. Construct the filter weights using the command,
>>[b,a]=butter(n,Wf);
where n is the order of the filter(N in the notes) and Wf is the position between zero and the Nyquist frequency where you want to make the cut. Note that Wf is two times the frequency in cycles per time step, since the Nyquist frequency is 0.5 cycles per time step. OK, you choose n and Wf. Note that butter seems to fail if n is too big or Wf too small. Now that we have chosen the coefficients we can calculate the frequency response of our filter as follows,
>>[h,w]=freqz(b,a,128);
>>r=abs(h);
>>rlow=r.*r;
>>f=w/pi/2.0;
>>plot(f,r,'--',f,rlow)
You should now have on screen a plot of the response function for the filter you have designed in the same frequency units as the power spectral plot you made earlier. The dashed line is the result of a single pass over the data, the solid line is the response function if you pass over the time series twice in the forward and backward directions to remove phase shifts associated with one pass of the filter. We will be using the double pass method, so the solid line is what you will get. If you are happy with this response function, then we are ready to continue. Make a plot of the response function of your chosen filter to include in your report. Anotate it with the parameters you used:
>>title('Butterworth response: N=n, cut at Wf')
where you have to type in the numbers n and Wf that you used.
>>print
4) High pass filter the data set: Now it is trivial to filter the original time series and construct the low-pass and high-pass filtered data sets.
>>xlow=filtfilt(b,a,x);
>>xhi=x-xlow;
>>t=0:1:4095;
>>plot(t(100:500),xlow(100:500),'--',t(100:500),xhi(100:500))
Now you have a plot of the low-passed(dashed) and the high-passed data stream as a function of time step between time steps 100 and 500. The low-passed is dashed and the high is solid. I am beginning to see some relationship between the high-passed and low-passed data, but it might be better if there were more separation between the frequencies. Let's do an even more low-passed cut with a smaller Wf than before and compare this with the high-passed data.
>>[b,a]=butter(n,Wfsmaller);
>>xlower=filtfilt(b,a,x)
>>plot(t(100:500),xlower(100:500),'--',t(100:500),xhi(100:500))
Look at other parts of the time series to see if the relationship
holds up. Can you imagine a compositing procedure that would
allow you to test your developing hypothesis statistically?
5.) Band-pass filter: Construct a band pass filter to pass only
the variance peak around 5 days and plot its response function.
6.) What have you learned about the 5-day signal and its relation
to the low frequency variance? Could you have learned this from
Fourier spectral analysis?
Part II: Real Data
A. 3-Million-Year Sediment Record
In the same ftp directory you will find a file called raymo.dat.
This is a more than three million year record of oxygen isotopes
from sediment cores drilled in the ocean bottom. The oxygen isotope
measures the abundance of oxygen-18 in seawater, which is related
to the amount of water tied up in ice sheets on land. So it is
a proxy measure of ice ages. I have interpolated the data to
a regular time interval of 3,500 years. This corresponds roughly
to the original spacing of the data I got from Raymo, M. E., 1992:
Global Climate Change: A three million year perspective. Start
of a Glacial, Proceedings of the Mallorca NATO ARW, NATO ASI Series
I, Vol. 3. 3, G. Kukla and E. Went, Eds., Springer-Verlag,
207-233, if you are really interested in understanding this data
set. Here we are just doing data analysis. We have an apriori
expectation of peaks at 100,000 years and 40,000 years. There
are about 850 samples in the original data, and a few more in
my interpolated set. Note that time decreases into the past, with
the first values for the present and the oldest values at the end of
the time series. This is the nutty way paleo people organize time
series data.
1.) Load the data, plot the data.
2) Calculate the power spectrum with a view toward separating
the 100,000 and 40,000 year periodicities. Can you show separate
significant peaks at these two periods using a priori statistics?
3) Use the Butterworth filter to band pass the 40,000 year variability.
Plot the band-passed data. Try to make the cut near the minimum
between the 100,000-year peak and the 40,000-year peak. What
do you see?
4) Use the Butterworth filter to band pass, or low pass the 100,000
year variability. Plot the band-passed data. What do you notice?
What does filtering tell you that spectral analysis alone does
not?
B. Central England Temperature Record
I put a file called cetml1659on.txt into the ftp directory. The first column is the year, then 12 months from January to December, then the annual average.
1. Try to find a spectral peak.
2. Try to filter to show the 65-year period, or any other period you think is interesting.
3. Experiment with different filtrations and present one that you think is the most reasonable.
As usual, make a report.
End