results_sst_precip_Johnna


 * What is the fundamental (lowest) frequency possible in this time series? (2*pi)/240months or 1 cycle per 20 years
 * What is the bandwidth or spectral resolution (Δf) of the spectrum you will create from it? (2*pi)/240months
 * What is the Nyquist (highest resolvable) frequency in this time series? 1 cycle / (2Δt), where Δt is 1 month, corresponding to 1 cycle / 2 months.
 * Based on the above, make a 1D frequency array f to use as the x axis on your plots in part 2. f = [1:240months/2]/240months

== 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 == **4. Plot the spectrum** as f*power vs. log(period). Area under the curve should still be proportional to total power (variance).
 * 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".
 * 2. Plot the spectrum** as an indefinite integral (cumulative power) vs. period or log period.
 * 3. Plot the spectrum** as f*power vs. log(f). Area under the curve should still be proportional to total power (variance).
 * 5. Plot the spectrum**as log(Pow) vs. log(f).
 * 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[] aka [|Flicker noise] aka[] . If the slope is -2, you have [|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 [|implications of][|finding a power law may be less profound than they appear].
 * 6. Plot the spectrum in your favorite format, after rebinning** f and Pow to coarser frequency bins. (Rebinning commands are in HW3 question 4).


 * 7. 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.
 * 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)
 * Also make a new frequency array corresponding to this longer series.
 * At high frequencies, the two should be almost identical. They are pretty close!


 * 8. 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.).

== 3. Significance testing of peaks: Overplot a red noise spectrum and its 99% significance level (the F test). ==
 * Estimate **your lag-1 autocorrelation value r1 for your field1.**
 * r1=0.9959


 * Use r1 to **create and overplot** the power spectrum of an autoregressive (AR(1)) or "red noise" process with the same r1 and same variance as your series.
 * 1) Use the [|F test] to**overplot**a line indicating the 99% significance level for spectral peaks.

== 4. Cross-spectrum of your 2 variables. (see section 3.1.2 of Hsieh handout). == **2. Compute the cross-spectrum** by complex multiplication: Cross = xhat .* conj(yhat) **4. Plot a cumulative spectrum of the in-phase part P (or R) as in 2.2 above.** Show that it ends up at the covariance of ts1 and ts2, mean( (ts1-mean(ts1)) * (ts2-mean(ts2)) ). What timescales contribute the most to the overall correlation (or covariance) of your 2 time series, according to this plot?
 * 1. Compute the 2 FFTs** of your 2 time series. xhat = fft(ts1); yhat = fft(ts2); **and plot the spectrum of your field 2** in your favorite depiction from Part 2.6 above, if you didn't already.
 * See Above, already plotted.
 * 3. Separate**the cross spectrum into its real and imaginary parts.
 * R and I in Hsieh (handout) section 3.1.2
 * I often see them called P(f) and Q(f) (the "in-phase" and "quadrature" parts). Quadrature means 90 degrees out of phase: sin and cos components.
 * R_Cumulative(end) =mean((ts1-mean(ts1)).*(ts2-mean(ts2))) = 2.2155

**5. Plot the squared coherency spectrum** (or just "coherence" in lazy language you will often hear). It is (P2+Q2) /(xPow2) /(yPow2). Why is it always 1?? (Hsieh 3.33-3.36)


 * 6. Plot the squared coherency spectrum after rebinning**P, Q, xPow and yPow to coarser frequency bins.
 * Reasoning: For physically real phenomena operating in a general, broad frequency band (like ENSO), the variables x and y will have the same phase relationship for all frequencies since they are physically linked, so averaging (rebinning) won't decrease coherency much. For random, non physically linked fluctuations of x and y, the in-phase and quadrature components will both be random (positive at one frequency, negative at the next), so the averaging will weaken the coherency when the phase relationship is random. See Hsieh handout, (3.37).




 * 7. Plot the phase spectrum**arctan(Q/P), and interpret the phase relationship between your variables in a frequency band where there is strong coherence (like ENSO) by showing this phase relationship using time series plots zoomed in to one dominant case of this strong oscillation (like an ENSO event).

<span style="font-size: 1.3em; margin: 0px; padding-bottom: 0px; padding-left: 0px; padding-right: 0px; padding-top: 5px;"> 5. 2D FFT
**1. Compute the 2D FFT (xhat) of your primary field's time-longitude section (fft2 in Matlab,** **fft in IDL).**
 * 2. Display this 2D fft**as an image or contour plot, as a function of frequency f and zonal wavenumber k (the number of wavelengths per 2.5 degree).
 * You have an f array: you just need to make a k array, following the logic in question 1 above in the x direction.
 * Shift the Power array (with fftshift in Matlab, or the shifting keyword in IDL8's FFT routine) to put low frequencies in the middle of the image.


 * Can you interpret this spectrum directly in terms of the size and orientation of stripes seen in the time-longitude section image?
 * Again we see that the spectrum depends more on time than space.

**3. Show a filtered data image:** Multiply the xhat 2-dimensional array by zero (mask it out) in some spectral region. Maybe low pass (exclude wavenumbers and frequencies > 10 times the fundamental). Or try band pass: can you find and zero out the El Nino peak? Anything that interests you. Explore, then teach us what you learned. Transform back to x-t space with ifft in Matlab or fft(xhat, /inverse) and display filtered and unfiltered data.

**4. Extra credit:** Explore filtering a bit more, and figure out better how the 2D spectrum relates to the rebinned variance diagram from HW3, problem 4.t