results_uwnd_slp_Angela


 * All plots were chosen with the subset at u-wind 1000 and SLP at longitude 160W (East Pacific).**

=**PART 1**=
 * 1) What is the fundamental (lowest) frequency possible in this time series? **The fundamental (lowest) frequency possible is 1 cycle per T, where T is total time in the series. For this particular data set that equates to 1 cycle per 20 years.**
 * 2) What is the //bandwidth// or spectral resolution (Δf) of the spectrum you will create from it? **The bandwidth will be the same as the fundamental frequency, or 2*pi divided by 240 months. Thus, the longer the record (T) the sharper spectral peaks that can be allowed without blurring two distinct frequencies into one.**
 * 3) What is the Nyquist (highest resolvable) frequency in this time series? **The Nyquist (highest) frequency possible is 1 cycle per 2*(delta t), where delta t is 1 month. For this particular data set that equates to 1 cycle per 2 months.**
 * 4) Based on the above, make a 1D frequency array f to use as the x axis on your plots in part 2. **The frequency array f, needs to encompass all possible frequencies from 2 months to 20 years. f=[low:low:high], where low is (1), the interval is the bandwidth (2), and high is (3). Please note: this is only half the frequency spectrum.**

=PART 2= 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 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 do so above)


 * Plot the spectrum** as a power spectral density PSD = Δ(variance)/Δf = Pow//Δf vs. frequency f.


 * As can be seen in these two plots, there is a large annual cycle signal for SLP, but not one for u-wind. Instead what dominates u-wind 1000 is something on a longer timescale with lots of little peaks on the longer wavelength end (smaller frequencies). The largest peak is around 60 months, or ~5 years, this has to be ENSO related!!**


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


 * This plot is helpful to determine if you calculated your power properly as the first value (longest wavelength), or total sum, should equal your variance. Mine checks out (taking into account slight differences from rounding errors) with a total power sum of around 1.6 for u-wind 1000 and 0.8 for SLP. This helps visualize where the largest changes in power occur. For u-wind, you can see that it is longer periods (months) that have the largest effect on power, suggesting that climate features, such as ENSO are playing a large role. On the other hand, for SLP it is clear that the annual cycle has the largest effect on power.**


 * Plot the spectrum** as f*power vs. log(f).


 * Notice that the plots do not appear to follow the same pattern as the frequency. This is simply due to the fact that if the frequency is less than one (which it is for the longer wavelengths), the log is actually negative. However, if you look closely, you can still find the matching peaks at certain times.**


 * Plot the spectrum** as f*power vs. log(period).


 * Now these are a little easier to read than the log(f) plot above, since there is no switch between negative and positive in the log period scale.**


 * Plot the spectrum** as log(Pow) vs. log(f).


 * Surprise, the plot looks like it's upside down! Well, when using log, anything less than 1 has a negative log, which is all of the power per frequency bin values (for both) and the slower time time scales. This makes it a little hard to read, but as you can see, no prominent straight line occurs. Therefore, these two variables most likely do not follow a power law.**
 * 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].


 * Plot the spectrum in your favorite format, after rebinning** f and Pow to coarser frequency bins. (Rebinning commands are in HW3 question 4).


 * I chose to plot the PSD again after a rebinning, because this plot makes the most sense to me and is the easiest for me to read. I rebinned by a factor of 3, because I wanted to keep most of the peaks and simply filter out some of the smaller variability. It is clear for both that they have a seasonal peak (time=6 months). SLP has a strong annual signal and then some, possibly ENSO related or PDO related variability for longer time scales. For u-wind 1000, the largest variability by far occurs on long time scales (approximately 2-7 years) which shows that it is governed by ENSO, rather than seasonal or annual cycles.**


 * lot 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.


 * I am guessing here. I found the fft from the padded ts1 and rebinned by a factor of 3, so that it could be compared to the original f array. Then I multiplied the end power spectrum by a factor of (var(ts1)/var(ts1padded)). Thus, the area under the curve of both power spectrums have almost identical variance (area under the curve). Who knows if this right (I doubt it), but the results are comparable at high frequencies. Help would be appreciated! (Blue is original, Red is padded).**
 * BRIAN SEZ: LOOKS GOOD TO ME. YES HIGH FREQUENCIES HARDLY CARE ABOUT THE PADDING AS THERE ARE LOTS OF REALIZATIONS DEEP WITHIN THE TIME SERIES. LOW FREQUENCIES ARE MORE DISTORTED BY EDGE EFFECTS AND THUS ARE MOST ALTERED WHEN YOU DO SOME PADDING OF THE ENDS.**

==3. Significance testing of peaks: Overplot a red noise spectrum and its 95% significance level (the F test).==
 * 1) **Estimate** your lag-1 autocorrelation value r1 for your field1. **r1 = 0.7703 (u-wind); r2=0.5602 (SLP) The autocorrelation is plotted below for u-wind and SLP.**


 * 1) 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. **(RED LINE)**
 * 2) Use the [|F test] to **overplot** a line indicating the 99% significance level for spectral peaks. **(GREEN LINE)**


 * At the 99% confidence interval, for u-wind, the most prominent peak at 60 months is not deemed significant (not sure if I believe this!), but the 6 month time spike is significant as well as some of the shorter time scale events (high frequencies). For SLP, both the annual and semi-annual time spikes are significant as well as some higher frequency events. (Thanks to Emily for an idea on how to scale these two lines!)**

===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//===
 * //For continuous AR(1) red noise, the autocorrelation r(lag) decays away exponentially with lag: r = exp(-lag/T)//
 * //The time constant T// //for this exponential decay is thus T = -(//Δ//t)/ln( r1 ) where r1 is your lag-1 autocorrelation.//
 * 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.
 * For the F test 99% significance threshold, the curve you plot is just your red noise null hypothesis curve times a factor given in [|this table] (from appendix G of the [|vonStorch and Zwiers ebook]).
 * 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.
 * Much more info and discussion (by the professor I learned it from...): www.atmos.washington.edu/~dennis/552_Notes_6b.pdf

==4. Cross-spectrum of your 2 variables. (see section 3.1.2 of Hsieh handout). (Thanks Changheng for help on these last two sections!)== > **A majority of the power (covariance) occurs on longer timescales (2-5 years), which would link the variability to ENSO. There is also an annual (mainly from SLP) and semi-annual spike that occurs. To isolate ENSO variability, then you just need to remove the annual/semi-annual values before taking the fft.** >
 * 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.
 * 2) **Compute the cross-spectrum** by complex multiplication: Cross = xhat .* conj(yhat)
 * 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.
 * 1) **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?


 * Ends at -0.1105, same as the mean of the multiplication of the two anomalies or covariance. This just shows the last half of the data.**
 * 1) **Plot the squared coherency spectrum** (or just "coherence" in lazy language you will often hear). It is (P 2 +Q 2 ) /(xPow 2 ) /(yPow 2 ). Why is it always 1?? (Hsieh 3.33-3.36)


 * The coherency is always one, because the wavenumbers will describe the same wave structure in both features (regardless of phase offset). Thus, if it's the same wavenumber, then the value will be 1.**




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


 * After rebinning by a factor of 3 (as done before), the coherency is not longer equal to one, because the waves do not necessarily share the same structure anymore.**




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


 * The phase spectrum shows, unsurprisingly, a lot of phase offset for the entire length of the spectrum. As from question 5 in HW3, to get a strong covariance, a longitude offset had to be taken into account to show the effects of ENSO on the two variables. Thus, a weak covariance exists when both features are on the same longitude, as is the case here, and phase offsets occur between the various wave structures.**

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.
 * 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.
 * You may want to display the log or square root of power, so the lowest frequencies don't dominate so strongly and you see more structure.
 * You should probably zoom in on the lowest frequency region.

> Can you interpret this spectrum in terms of the rebinned variance diagram from HW3, problem 4?
 * 1) Can you interpret this spectrum directly in terms of the size and orientation of stripes seen in the time-longitude section image? >>
 * The power spectrum shows that the variance decreases in time faster than in space (as did the plots from HW3). However, I'm not sure I can relate the stripes to the raw data. I think I need some practice looking at these types of spectrums more, but it does look rather awesome.**


 * 3. Show a filtered data image**: Multiply the xhat 2-dimensional array by zero wherever k>10 or f>10. Transform back to x-t space with ifft in Matlab or fft(xhat, /inverse). Display.


 * The low band pass filter removes El Nino/La Nina related events making the whole field more semi-annual in structure with less emphasis on lower frequency events.**


 * The high band pass filter removes the high frequency variability, but leaves in the ENSO-related signals, resulting in an overall smoother picture. Pretty neat if you ask me.**