Results_uwnd_precip_Emily

**help with part 1:**
I suffered a bit as well. YOu need to make sure that total(power) is the same for all the different kinds of plots. And check that variance is conserved. For IDL users I did the following: pow = abs(fft(x90_tseries))^2 ;;;x90_tseries is just your personal time series PSD = pow/delta_f ;;;;Just take positive half & double it nt = n_elements(x90_tseries) PSD_half = PSD(0:nt/2-1)

;;;CHECK TO MAKE SURE total variance conserve. Notice some funny things here. Just doubling the PSD_half isn't right! ANYWAY, MAYBE THIS IS MORE DISTRACTING THAN HELPFUL, SORRY. total, PSD 2.57788 total, PSD_half*2 2.94272 total, PSD_half[1:*]*2 + psd_half[0] 2.57788 total, PSD_half[1:*]*2 2.21303 total, var_n(x90_tseries) 2.2130355

So, i did the following for psd_half PSD_half[1:*] = PSD_half[1:*]*2

If you remove the mean to begin with, none of these problems arise. total, PSD 2.21304 total, PSD_half*2 2.21303 total, PSD_half[1:*]*2 + psd_half[0] 2.21303 total, PSD_half[1:*]*2 2.21303 total, var_n(x90_tseries) 2.2130355

=The time series I'm using is u-1000 mb wind at 90°E=

=Question 2:=
 * 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).**

** I HAVE NO IDEA IF WHAT I'M DOING IS RIGHT, SO DON'T LOOK TO ME FOR HEL ** P!
1. Power spectrum for u-1000 mb wind



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


 * Here you can see a jump in power at 6 months and 12 months, or the semi-annual and annual harmonics.

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

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


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

6. **Plot the spectrum in your favorite format, after rebinning** f and Pow to coarser frequency bins.
 * I chose to rebin by 4 months.

7. **Plot the spectrum in your favorite format,** **overplotting** the PSD you get //when you **pad the ends of the time series with zeroes**//. No idea if this is correct.


 * I multiplied the padded time series by 10 to get comparable power values. Then I rebinned the padded time series to 120, so the original power spectrum and the padded one could be plotted on the same figure. Black is the padded time series, red is the original time series. Yes, they are similar at high frequencies.



3. Significance testing of peaks: Overplot a red noise spectrum and its 95% significance level (the F test).
1. I plotted the inverse fft of the symmetric power spectrum to plot the autocorrelation. NOT SURE IF THIS IS CORRECT.

**My Code is as follows:**
 * IDL >** autocor = real_part(fft(pow, /inverse))
 * IDL >** months = findgen(240) - 120

ytitle = 'm/s', thick = 4
 * IDL >** plot, months, shift(autocor, 120), title = 'Autocorrelation', xtitle = 'months', $
 * IDL >** hor, 0
 * IDL>** w = where(months eq -1)
 * IDL>** b = shift(autocor, 120)
 * IDL>** r1 = b[w]




 * At a lag of -1 the autocorrelation value = 1.62998.

2. 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. NOT SURE IF THIS IS CORRECT. **My Code is as follows:** ;;#2--create red noise dt = 1 T = -dt/alog(r1) p_red = 2*T[0]/(1 + T[0]^2*f^2)

;;;;Re-scale to have same total power (variance) as above rescaled_p_red = (p_red/total(p_red)*2.58) ;;;total(power) = 2.58


 * rescaled_p_red is what I have over plotted on the below figure.



3. Use the f-test to **overplot** a line indicating the 99% significance level for spectral peaks.


 * The dashed red curve indicates the 99% significance levels for the spectral peaks.
 * 0.02 and 0.05 cycles per month are above the significance line. That's roughly 1 cycle per 20 and 50 months.



4. Cross-spectrum of your 2 variables.

 * 1. Power spectrum of second time series, in my case precipitation: **




 * 2. Compute cross-spectrum:**
 * IDL>** xhat = fft(x90_tserie
 * IDL>** yhat = fft(y90_tseries)
 * IDL>** cross_spec = xhat * conj(yhat)


 * 3. Take real and imaginary parts:**
 * IDL >**Rcross = real_part(cross_spec)
 * IDL >** Icross = imaginary(cross_spec)


 * 4. Plot cross-spectrum**


 * A 20 month time scale shows the highest covariance.

5. **Plot the squared coherency spectrum** (or just "coherence" in lazy language you will often hear). It is (P 2 +Q 2 ) /(xPow) /(yPow). Why is it always 1?? code xhat = fft(x90_tseries) - mean(x90_tseries) yhat = fft(y90_tseries) - mean(x90_tseries) cross_spec = conj(xhat) * (yhat) pow = abs(xhat)^2 powy = abs(yhat)^2
 * My code is as follows: **
 * 2

Rcross = float(cross_spec) Icross = imaginary(cross_spec)
 * 3

coh_spectrum = (rcross^2+icross^2) /(Pow) /(Powy)
 * 5

code



6. **Plot the squared coherency spectrum after rebinning** P, Q, xPow and yPow to coarser frequency bins. ****Here I followed Brian's code to get my plot*********



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



5. 2D FFT

 * 2. Display this 2D fft** as an image or contour plot, as a function of frequency f and zonal wavenumber k.


 * I reversed the power so the stronger power is in the eastward direction (which I think makes more sense). The plot is the alog(power) to show more structure. The mean has also been removed.
 * Non-zoomed in plot:


 * Zoomed in plots:






 * Besides the strong power at very low frequencies, there are lobes of strong power at .018, .034, .055, and .085 cycles per month in both eastward and westward wavenumbers, though the eastward is stronger, except for .085 cycles per month. This frequency corresponds to about 1 cycle per 55 months, 1 cycle per 30 months, 1 cycle per 18 months, and 1 cycle per 12 months, respectively.

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

;;;;;;CODE x_prime = uwnd ; - mean(uwnd) y_prime = precip - mean(precip) xhat = fft(x_prime) x_shift = shift(xhat, 72, 120)

wk = where(abs(k) le 10) wf = where(abs(freq2) le .1) array = intarr(144, 240) array[min(wk):max(wk), min(wf):max(wf)] = 1

filtered_xhat = x_shift * array filtered_xhat = shift(filtered_xhat, -72, -120)

filtered_data = fft(filtered_xhat, /inverse)

!p.multi = [0, 2, 1] contour, filtered_data, lon, time, /fill , nlev = 20, title = 'filtered Uwnd (f<.1, k<10)' , $ xtitle = 'longitude', ytitle = 'time (years)'