results_sst_slp_Chen

//**Here are my codes for computing and plotting.**//
==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. (Note: T= 240 months, Δt=1 month) == > fmin = 1/T or 2 π /T > Δf = 1/T or 2 π /T > fmax = 1/(2Δt) or  π /Δt > f = [0:T/2-1]/T; % f = N/T, N=**0 **,1,2,...,T/2-1, including the 1st value (i.e. mean(ts1), OR > f = 2 π*[0:T/2-1]/T; % f = N/T, N=**0 **,1,2,...,T/2-1, including the 1st value (i.e. mean(ts1) > > f = [1:T/2]/T; % f = N/T, N=1,2,...,T/2, excluding the 1st value, OR > f = 2 π*[1:T/2]/T; % f = N/T, N=1,2,...,T/2, excluding the 1st value
 * 1) What is the fundamental (lowest) frequency possible in this time series?
 * 1) What is the //bandwidth// or spectral resolution (Δf) of the spectrum you will create from it?
 * 1) What is the Nyquist (highest resolvable) frequency in this time series?
 * 1) Based on the above, make a 1D frequency array f to use as the x axis on your plots in part 2.

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
Note: the area under the curve (i.e. variance) = PSD* Δf //X label should be 'frequency (**cycle/month**)', the same applies to the following figures wherever it is mislabeled as 'frequency (month_(-1))'//
 * 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. In the upper plot, several peaks (i.e., 60-, 30- and 12- month in period) stand out without mingling with each other; the plot becomes noisy when plotting scale is changed to log-log scale.

**2. Plot the spectrum** as an indefinite integral (cumulative power) vs. period or log period. Several jumps (e.g. at 12-, 30-, 60-month) of the curve come into notice. They indicate a relatively large portion of the total variance is contained in those periods.

Note: the area under the curve = variance*Δf
 * 3. Plot the spectrum** as f*power vs. log(f). Area under the curve should still be proportional to total power (variance).

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



This is just a mirror image of the the plot from the foregoing part.

**5. Plot the spectrum** as log(Pow) vs. log(f). The slope of curve in the high frequency part is close to -2, probably due to Brownian noise in the data.
 * 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). By doing so, one can get rid of the high-frequency fluctuations. Only the strong signals remain and show as distinct peaks. However, strong signal at certain frequency can weaken or even disappear.


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


 * Adjust the variance of the padded time series spectrum so that it overlays the unpadded spectrum well.


 * At high frequencies, the two should be almost identical.


 * 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 12 month (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).
> >> [r,lags] = xcorr(ts1,'coeff'); > r1 = r(T/2+1) = 0.9959 > T = -1/ln(r1) = 242.4646 > P_red(f) = 2T/(1 + f 2 T 2 )
 * 1) **Estimate** your lag-1 autocorrelation value r1 for your field1.
 * 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.
 * 1) Use the [|F test] to **overplot** a line indicating the 99% significance level for spectral peaks.

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 + f2T2). 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 highlighted 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).
**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) **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?
 * 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))) = -0.1392. Time scale of 12-month contributes most to the overall covariance.

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


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

Should it be symmetric in x-axis??? The phase between SST and SLP is 58 degree at the frequency 0.0813 (i.e. period about 12.3 months). Although covariance of the two variables is dominant at this period, phase spectrum shows they do not change in/out phase???

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.


 * 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 directly in terms of the size and orientation of stripes seen in the time-longitude section image?
 * The above figure shows that the 2-dimensional power spectrum of SST depends more on time than space. Zonally orientated stripes indicate the existence of four prominent frequencies (0.03, 0.0875, 0.17 and 0.25 month^(-1), or 34.3, 11.4, 5.85 and 1.25 months in period) in which dominant portion of temporal variation of SST are contained. Temporal variation of SST concentrates mostly in low wavenumbers.
 * Can you interpret this spectrum in terms of the rebinned variance diagram from HW3, problem 4?
 * The shape of the distribution of the 2-dimensional power spectrum of SST looks like a standing rhombus elongated along the axis of frequency. It indicates much more variance of SST is contained in high frequencies (small scale structures) than in high wavenumbers (small scale structures).
 * On the other hand, rebinning of the variance (from the previous homework) smears the small structures first as the rebinning scale factor increases.
 * Thus, variance of SST decreases faster in frequency domain than in wavenumber domain as rebinning scale factor increases. FFT analysis shown here coincides with the conclusion drawn from HW3, problem 4.

**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. The title of the right panel 'low band pass' and the one in the following figure 'high band pass' are bad terminologies; 'low band pass' is supposed to mean a set of low frequency is removed, while 'high band pass' a set of high frequency is removed.

As shown in the original SST distribution (left panel), reddish stripes sandwiched by blueish stripes at around 260 ° (i.e. 100°W, east Pacific) indicate the occurrences of El Nino. These reddish stripes correspond to the 82-83, 87-88, 91-92, 92-93 and 97-98 El Nino events. In contrast, a low band pass in frequency ((1/12)month -1 <f<(1/240)month -1, that is, 1-year to 20-year in period) is applied to the FFT of SST. Signals of El Nino events are wipped off and cannot be detected in filtered SST (right panel). The upper figure shows the results from a high band pass (2 months to 6 months in period). Although delicate structures are lost, El Nino signal is preserved.

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