Results_precip_u1000_David

=**FFT homework, using your primary and secondary fields.** =

==**1. Load in your data again, and extract 240-month time series ts1(t) and ts2(t) at some longitude ( the central Pacific, longitude index 80 = 200E = 160W will have a lot of ENSO signal, or pick your own). Take these series from array x and y, not anomx and anomy -- we'll keep the annual cycle in here to help make sure the frequency axis is right (there should be a distinctive peak there). The mean of your time series will appear in the real (cosine) zero frequency bin if you don't remove it.** ==
 * 1) 
 * 2) What is the fundamental (lowest) frequency possible in this time series? 1 cycle per 240 months, f = 1/T
 * 3) What is the //bandwidth// or spectral resolution (Δf) of the spectrum you will create from it? Same as above, 1 cycle per 240 months; Δf = 1/T
 * 4) What is the Nyquist (highest resolvable) frequency in this time series? fmax is 120 cycles per 240 months = 1 cycle per 2 months; fmax = 1/2T
 * 5) Based on the above, make a 1D frequency array f to use as the x axis on your plots in part 2. f = [1:T/2]/T

==**2. Make the 1D power spectrum of your first field time series:** == ==**power per frequency bin Pow = abs(fft(ts1)).^2, Pow = abs(fft(ts1))^2 ** == (Angela) Since I suffered A LOT to figure this out... here is the code to make your power per frequency bin work with the correct units (MATLAB):

X=fft(ts1); %fft of ts1 X(1)=[]; %the first value is the sum of all the other values, so you need to remove it **! the first value equal to the mean of ts1** fftx=X(1:120); %to split your fft in half (don't have to do this) mx=abs(fftx); %absolute value mx=mx/length(ts1); %to get the correct scaling of the fft for power calculation mx=mx.^2; %power mx=mx.*2; %only because I split my fft in half (don't do if did not split it above


 * 1) **Plot the spectrum** as a power spectral density PSD = Δ(variance)/Δf = Pow//Δf vs. frequency f. Label the axes with the right values and units. Area under the curve should be proportional to total power (total variance). Since it's variance in discrete bins, you should ideally use a bar plot or plotting symbols, not just a line plot connecting the "points".
 * You may want to center it on 0 frequency (by shifting the array) to show a symmetric spectrum with positive and negative frequencies.
 * Or you may prefer to just half the spectrum (PSD vs. the absolute value of f -- remember to double the positive frequency part of Pow so area = variance).
 * You may also choose to rebin f and PSD to coarser spectral bands, if the plot is too noisy.


 * 1) **Plot the spectrum** as an indefinite integral (cumulative power) vs. period or log period.


 * 1) **Plot the spectrum** as f*power vs. log(f). Area under the curve should still be proportional to total power (variance).


 * 1) **Plot the spectrum** as f*power vs. log(period). Area under the curve should still be proportional to total power (variance).


 * 1) <span style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0.5em; padding-bottom: 0px; padding-left: 3em; padding-right: 0px; padding-top: 0px;">**Plot the spectrum** as log(Pow) vs. log(f).
 * <span style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0.5em; padding-bottom: 0px; padding-left: 3em; padding-right: 0px; padding-top: 0px;">Why this way? Area under the curve is no longer meaningful. The reason to plot a spectrum this way is to see if it looks like a straight line. If the slope is -1, you have Pink Noise <span style="background-attachment: initial; background-clip: initial; background-color: initial; background-origin: initial; background-position: 100% 50%; background-repeat: no-repeat no-repeat; padding-right: 10px;">[] aka <span style="background-attachment: initial; background-clip: initial; background-color: initial; background-origin: initial; background-position: 100% 50%; background-repeat: no-repeat no-repeat; padding-right: 10px;">[|Flicker noise] aka <span style="background-attachment: initial; background-clip: initial; background-color: initial; background-origin: initial; background-position: 100% 50%; background-repeat: no-repeat no-repeat; padding-right: 10px;">[] . If the slope is -2, you have <span style="background-attachment: initial; background-clip: initial; background-color: initial; background-origin: initial; background-position: 100% 50%; background-repeat: no-repeat no-repeat; padding-right: 10px;">[|Brownian noise (hear an acoustic sample here!)] . Slopes of -3 or -5/3 are predicted for KE (velocity variance) by 3D and 2D turbulence theory (based solely on scaling arguments). No matter what the slope, a straight line implies a power law, although the <span style="background-attachment: initial; background-clip: initial; background-color: initial; background-origin: initial; background-position: 100% 50%; background-repeat: no-repeat no-repeat; padding-right: 10px;">[|implications of finding a power law may be less profound than they appear].

>>>
 * 1) <span style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0.5em; padding-bottom: 0px; padding-left: 3em; padding-right: 0px; padding-top: 0px;">**Plot the spectrum in your favorite format, after rebinning** f and Pow to coarser frequency bins. (Rebinning commands are in HW3 question 4).
 * 2) <span style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0.5em; padding-bottom: 0px; padding-left: 3em; padding-right: 0px; padding-top: 0px;">**Plot the spectrum in your favorite format,** **overplotting** the PSD you get //when you **pad the ends of the time series with zeroes**//. This will highlight the errors associated with making your time series as if it were periodic.
 * <span style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0.5em; padding-bottom: 0px; padding-left: 3em; padding-right: 0px; padding-top: 0px;">to do this part, just make a new data array tspad = [ts*0, ts-mean(ts), ts*0] (IDL) or [ts*0 ts-mean(ts) ts*0]; (Matlab)
 * <span style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0.5em; padding-bottom: 0px; padding-left: 3em; padding-right: 0px; padding-top: 0px;">also make a new frequency array corresponding to this longer series.
 * <span style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0.5em; padding-bottom: 0px; padding-left: 3em; padding-right: 0px; padding-top: 0px;">Adjust the variance of the padded time series spectrum so that it overlays the unpadded spectrum well.
 * <span style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0.5em; padding-bottom: 0px; padding-left: 3em; padding-right: 0px; padding-top: 0px;">At high frequencies, the two should be almost identical.




 * 1) <span style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0.5em; padding-bottom: 0px; padding-left: 3em; padding-right: 0px; padding-top: 0px;">**Extra** credit on 1D spectra: explore end explain some of the virtues of one of the MANY special built-in functions or packages for spectral analysis (periodogram or other special functions in Matlab, IDL wavelet GUI, spectraworks.com package for Mac, etc etc.).

==<span style="font-size: 1.3em; margin: 0px; padding-bottom: 0px; padding-left: 0px; padding-right: 0px; padding-top: 5px;">**3. Significance testing of peaks: Overplot a red noise spectrum and its 95% significance level (the F test).** ==
 * 1) <span style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0.5em; padding-bottom: 0px; padding-left: 3em; padding-right: 0px; padding-top: 0px;">**Estimate** your lag-1 autocorrelation value r1 for your field1.
 * 2) <span style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0.5em; padding-bottom: 0px; padding-left: 3em; padding-right: 0px; padding-top: 0px;">Use r1 to **create and** power spectrum of an autoregressive (AR(1)) or "red noise" process with the same r1 and same variance as your series.
 * 3) <span style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0.5em; padding-bottom: 0px; padding-left: 3em; padding-right: 0px; padding-top: 0px;">Use the <span style="background-attachment: initial; background-clip: initial; background-color: initial; background-origin: initial; background-position: 100% 50%; background-repeat: no-repeat no-repeat; padding-right: 10px;">[|F test] to **overplot** a line indicating the 99% significance level for spectral peaks.



===<span style="font-size: 1.1em; margin: 0px; padding-bottom: 0px; padding-left: 0px; padding-right: 0px; padding-top: 5px;">**Background info for 3.2 and 3.3: //creating the power spectrum of the AR(1) process, your null hypothesis ("Red Noise"), and its F test//** ===
 * <span style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0.5em; padding-bottom: 0px; padding-left: 3em; padding-right: 0px; padding-top: 0px;">//For continuous AR(1) red noise, the autocorrelation r(lag) decays away exponentially with lag: r = exp(-lag/T)//
 * <span style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0.5em; padding-bottom: 0px; padding-left: 3em; padding-right: 0px; padding-top: 0px;">//The time constant T// //for this exponential decay is thus T = -(//Δ//t)/ln( r1 ) where r1 is your lag-1 autocorrelation.//
 * <span style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0.5em; padding-bottom: 0px; padding-left: 3em; padding-right: 0px; padding-top: 0px;">The power spectrum of red noise is P_red(f) = 2T/(1 + f 2 T 2 ). Compute this from your f array and rescale its total variance to match yours.
 * <span style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0.5em; padding-bottom: 0px; padding-left: 3em; padding-right: 0px; padding-top: 0px;">For the F test 99% significance threshold, the curve you plot is just your red noise null hypothesis curve times a factor given in <span style="background-attachment: initial; background-clip: initial; background-color: initial; background-origin: initial; background-position: 100% 50%; background-repeat: no-repeat no-repeat; padding-right: 10px;">[|this table] ([|Ftest99.vonStorchZwiers.png]) (from appendix G of the <span style="background-attachment: initial; background-clip: initial; background-color: initial; background-origin: initial; background-position: 100% 50%; background-repeat: no-repeat no-repeat; padding-right: 10px;">[|vonStorch and Zwiers ebook] ).
 * <span style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0.5em; padding-bottom: 0px; padding-left: 3em; padding-right: 0px; padding-top: 0px;">Background: The F statistic applies to the ratio of variances between 2 processes. The null hypothesis process (red noise here) is exact and analytic, not estimated from data, so it has "infinite" degrees of freedom. Your data spectrum has 2 degrees of freedom (DOFs) per fundamental frequency interval (bandwidth), so I highlighed the 2 DOF number in the table. If you rebin your spectrum to a coarser bandwidth, then you have 4 or 6 or 8 or 10 or more DOFs per bin. Noise will tend to disappear with this bin-averaging, so the threshold for a peak to achieve statistical significance gets lower (you get to use a smaller factor from the F test table, there is another page in the book for DOFs greater than 10). Real peaks arising from physical processes (like ENSO) will tend to produce variance in many nearby frequencies, so a real spectral peak will not shrink as fast with rebinning (or averaging), and may still exceed the threshold for 99% significance.


 * <span style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0.5em; padding-bottom: 0px; padding-left: 3em; padding-right: 0px; padding-top: 0px;">Much more info and discussion (by the professor I learned it from...): www.atmos.washington.edu/~dennis/552_Notes_6b.pdf