\input vanilla.sty
\input seteps
\font \bigbf=cmb10 scaled\magstep2

\line{}
\vskip -2 cm
\seteps{-1.5 cm}{4.0 cm}{4.0 cm}{c:/stood/481/class/logo481.eps}
\vskip -3 cm


\line{\hskip 3 em\hfill \bigbf LAB 6 \hfill}
\vskip 1 em
\line{\hskip 3 em\hfill\bigbf FOURIER TRANSFORMS II \hfill}
\line{\hskip 3 em\hfill\bigbf \hfill}

\vskip 2 cm


\vskip 1 em\item{} Get LABFTII.mat from PAWS and load LABFTII.mat. This will create the data files (described below) in your workspace that you will need to do the lab. In each of the following exercises, work out the lowest frequency, the highest frequency (Nyquist), the increments to the frequency scale, and the uncertainty  on the frequency,  before doing the transform. Do the same for periods. In each example, experiment with the drift filter and the tapering.

\vskip 1 em\item{1)} \it GRAV \rm is synthetic gravity tide record for Ottawa, hourly samples, in microgal. Fourier transform and tell me what you see.
  
\vskip 1 em\item{2)} \it STCKHLM \rm contains the annual mean sea-level observations (in cm)
at Stockholm from 1774-1984. The first column is the year, the second is
the sea-level. There is a nearly linear decrease in sea-level, of nearly
half a cm/yr, ie a drift. In this case, the drift is actually more
interesting than the residual. You should know the cause.  Fourier transform the sea
level data anyway and see if there is anything interesting in it. There
are small 18.6 yr and 9.3 yr tides, as well a Chandler `pole tide', which
at 14 month period is aliased.  What frequency is it 
aliased to?  This is a very noisy spectrum, and the noise may be caused
by a variety of things, the smallest fluctuations, including the aliased
Chandler are lost in the noise.  

\vskip 1 em\item{3)} \it CO2 \rm contains monthly observations of $CO_2$ at Mauna Loa. Column 1 is year, column 2 is CO2 in ppm and column3 is column 2 minus an annual term. This is a famous data set that is often interpreted as a warning sign for global warming.  Some invalid observations have been replaced by interpolation. There are the invalid data to deal with,  a drift, and obviously an annual signal, which you can see just by looking at the data. Fourier transform the raw data and interpret?

\vskip 1 em\item{4)} \it ILS \rm contains some  latitude observations of the International Latitude Service from 1900 to 1979, at intervals of 1/20 of a year (each row is a year, with twenty equispaced samples in that year). The values are the change in latitude along the Greenwich meridian, in thousandths of a sec of arc. QFT wants to read a single column of data so  do \it ILS=reshape(ILS',76*20). \rm Compute the spectrum of about 10 years of data, and describe
it. Compute the spectrum of the whole data set. What is different, and why? 


\vskip 1 em\item{5)} \it KINGSTN \rm and \it BURLNTN.DAT \rm are hourly lake level data from Kingston, at the
east end of Lake Ontario, and Burlington near the west end. Fourier
transform these and interpret. In the raw data you should see one large
signal that only lasts for a short time. Estimate the frequency. In the
transform, you should see this signal and another, perhaps slightly larger
in amplitude. Why is the latter similar or slightly larger	than the
former in the spectrum, and almost unobservable in the raw record?

\vskip 1 em\item{6)} \it sunspot \rm is a monthly record of the sunspot index from 1749 to 2003. Fourier transform these data and interpret. 


\vskip 1 em\item{7)} \it belany \rm is a record from the Belany Islands earthquake  of March 25 1998, recorded by the superconducting gravimeter at Cantley, Quebec. The first column is the gravity channel output in volts, the second column is the pressure channel output in volts and the third channel is the gravity earth tide at Cantley. The sample interval is 1 minute. Observed gravity, corrected for tides and atmospheric pressure is $g=-78.5\mu$ $gal/V \times col1 - col3 +0.3\mu gal/mb \times 60.08 mb/V \times col2$. Convert all three columns to gravity and look at them separately, as well as gravity after each correction is applied. The final residual contains the surface waves and free oscillations generated by the Belany Island event.  One trick to beat down the noise and make the small signals more noticeable would be to divide
the record up into (say) four segments of equal length, Fourier transform
each segment and compute the power (that is ${real}^2 + {imag}^2$). Finally,  add (stack) the separate power transforms. The reason for computing the power is that the phase in each segment is reckoned from the start of the segment, so even with a pure sinusoid the phase in each segment will differ. By computing the power you make the phase irrelevant, so that when you add the transforms together any \it frequency \rm that is common to all segments will stack. If you tried to stack the phase, you would find that the phase of your signal is random from segment to segment, just like the noise, and will not stack.   This averaging (it is exactly the
same as stacking or averaging in the time domain to suppress noise)  will tend to reduce the amplitude of incoherent noise, perhaps revealing something else (no guarantees).
 

\vskip 1 em\item{8)} Extract a long road Bouguer profile from the complete field school data set. Fourier transform your gravity profile (you may need
to interpolate to get evenly spaced samples). The spectrum should
peak at low frequencies and tail off to high frequencies. What level does
the amplitude spectrum approach as you get towards the Nyquist?  What does
this mean?

\vfill
\eject

\bye 


\heading
ANSWERS TO LAB II
\endheading

\vskip 1 em\item{1)} You should get diurnal, semidiurnal, terdiurnal, quartdiurnal, monthly, and 
semimonthly signals. The monthly is just barely short enough period to be
handled properly. There are annual and semiannual signals as well, which, over
this short record constitute a small drift.
  
\vskip 1 em\item{2)} The linear decrease in sea level is probably mostly post-glacial rebound. This area is still rebounding from the weight of the Fenno-Scandian ice sheet. To work out the aliased frequency suppose we found two sinusoids with frequencies
$\omega_a$  and $\omega_h$ where $\omega_a$ is the aliased version of some high frequency $\omega_h$, that satisfied
$$
e^{i(\omega_a \Delta T)} =e^{i(\pm\omega_h  +{2n\pi\over \Delta T})\times
\Delta T}
$$
where n and m are integers,
so that we are sampling at integer multiples of $\Delta T$, and at
each sample the sinusoids are identical. The
exponentials give the change in phase of each wave. Therefore,
$$
(\omega_a \Delta T) =(\pm\omega_h\times  +{2n\pi\over \Delta T})\times
\Delta T
$$
\item{} or, with circular instead of radial frequencies 
$$
(f_a \Delta T) =(\pm f_h\times  +{n\over \Delta T})\times
\Delta T
$$
\item{} or
$$
f_a  =\pm f_h  +{n\over \Delta T}
$$
\item{} or
$$
f_a  =-f_h  +{2n\over f_N}
$$
\vskip 1 em\item{} where $f_n$ is the Nyquist frequency.	The choice of the $+$ sign would
lead to aliased frequencies outside of our Fourier bandwidth, so we drop
this choice. The file alias.m in geol481 calculate the aliased frequency
of a high frequency signal.

\vskip 1 em\item{} In other words the aliased frequencies of $\omega_h$, are
$\omega_a$, as given by this formula, because they result in the same
samples of the two sinusoids.
\vskip 1 em\item{} Rearranging
$$
f_a=f_N-(f_h-f_N) \quad\text{ with n=1}
$$
\vskip 1 em\item{} So $f_a$ is obtained by folding the frequency axis at $f_N$ and $f_h$
folds over to $f_a$. If $f_h$ is high enough that folding the frequency
axis at $f_N$ aliases $F_h$ to negative frequencies, the folding the axis
again at 0 will bring the alias back to the positive side of the axis.
The twice folded frequency may still be above
the Nyquist, so folding again about the Nyquist, (and again if necessary
about 0) and so on will eventually bring the aliased frequency between 0
and the Nyquist, where it may cause a problem. 

\vskip 1 em\item{}  If the Chandler frequency is .845 cy/yr, and we are sampling at 1 year, then
the aliased frequency is 
$$
\sigma _a = .5 cy/yr -(.845\, cy/yr -.5 cy/yr)= \pm 0.155cy/yr
$$
\vskip 1 em\item{} There is an m-file \it alias.m \rm that will calculate aliases.
You will probably not see a peak at the aliased Chandler frequency unless
you do the stacking suggested in the assignment.
  
\vskip 1 em\item{3)} There is of course an annual signal, and the transform reveals a semi-annual
signal as well, which is not that easy to see in the data. There does not appear
to be much else, other than the seasonal fluctuations, and the secular increase.
There are atmospheric phenomena in this period range, including El-Ni\~no, but
they do not seem to involve  $CO_2$.

\vskip 1 em\item{4)} This motion is dominated by the 14 month Chandler wobble, which is the only
peak you see in the 10 year data set, but when you use all the data, the single
peak splits into a 14 month (Chandler) term, and a yearly term. These are too
close in period for the transform to separate with only ten years of data. Why?
We have about 200 samples in ten years, so the frequency increments are
$$
\Delta f = {1\over (2\times 1/20\, (yr)\times 100})= .4\, cy/yr 
$$
\vskip 1 em\item{} The Chandler frequency is .845 cy/yr and the annual 1, so they are closer 
together (.155 cy/yr) than the spacing on the frequency scale! If we increase the
number of samples to 1520, or 76 years at  $20^{th}$ yr sampling rate, the
frequency increments are .05 cy/yr, so we can resolve  two separate peaks. At
this sampling rate  (20/yr) we would need 30 or more years to resolve the two
peaks.  You may also see a semi-annual peak, caused (like the annual) by the
atmosphere moving about, and perhaps a 24 year period, which some people say is
caused by the ocean, others by the fluid core, and others by the observers
themselves!

\vskip 1 em\item{5)} A seiche is a 'bathtub' oscillation of a lake, bay, or harbour. The
period depends on the length and depth of the body of water. They are most
often caused by wind storms that pile a lot of water up at one end of the
lake. If the storm abates faster than the water can flow back, the lake
will go into a few cycles of sloshing back and forth. For an east
west seiche in lake Ontario the period is about four or five hours, with a
maximum amplitude of about half a meter, and three cycles is typical.  

\vskip 1 em\item{6)} The amplitude spectrum reveals an 11 year cycle, as well as a 100 year
cycle and a 5 or 6 year cycle. The 11 year cycle is the famous sunspot
cycle, but the 100 year and 6 year cycles are less well known, and much
less prominent in the spectrum. 

\vskip 1 em\item{7)} The amplitude spectrum should tail off to the noise level in the data
by the time you get to the Nyquist. What did you estimate from the
repeats, and what do you estimate here. They are probably not the same,
because the repeats measure only one component of 'noise' that due to
operator, or poor drift correction. The 'noise' in the spectrum estimates the 
actual short scale gravity variations that may be present as well, and are
'noise' for your survey. 
\bye
