| 1 | FUNCTION spectra, datin, nyears, tukey, running, full, bootwin |
|---|
| 2 | |
|---|
| 3 | ; -------------------------------------------------------------------------- |
|---|
| 4 | ; Program to extract data, compute anomaly, do time series analysis. |
|---|
| 5 | ; Computes detrended data using a linear regression and |
|---|
| 6 | ; produces autocovariances and correlograms. |
|---|
| 7 | ; Then computes spectrum with and without |
|---|
| 8 | ; the Tukey window for various values of m. |
|---|
| 9 | ; |
|---|
| 10 | ; Adapted by Eric Guilyardi from a fortran program by Julia Slingo |
|---|
| 11 | ; (c) CGAM January 2002 |
|---|
| 12 | |
|---|
| 13 | ; Input: |
|---|
| 14 | ; ------ |
|---|
| 15 | ; datin : 1D input data (monthly time serie) |
|---|
| 16 | ; nyears: number of years in data |
|---|
| 17 | ; tukey : tukey window in years (20 as a first guess for 100 years) |
|---|
| 18 | ; running: running window in years for standard deviation |
|---|
| 19 | ; variations (15 first guess) |
|---|
| 20 | ; full: 1 = no anomaly, 0: anomaly wrt SC |
|---|
| 21 | ; bootwin = bootstrap window for error bar calculation (in months) |
|---|
| 22 | |
|---|
| 23 | ; Output: |
|---|
| 24 | ; ------- |
|---|
| 25 | ; spectra: structure { std: standard deviation |
|---|
| 26 | ; mean: mean of serie |
|---|
| 27 | ; sc: mean seasonal cycle |
|---|
| 28 | ; sc_range: range of mean SC |
|---|
| 29 | ; sc_std: seasonal cycle STD |
|---|
| 30 | ; anom: anomaly time serie (non-normalized) |
|---|
| 31 | ; spec: spectrum |
|---|
| 32 | ; spectw: normalized spectrum with Tukey window |
|---|
| 33 | ; time: spectra time array |
|---|
| 34 | ; fmax: frequency of maximum power |
|---|
| 35 | ; run_std: running standard deviation |
|---|
| 36 | ; run_sc: running seasonal cycle |
|---|
| 37 | ; anom_stdsc: anomaly time serie using run_sc |
|---|
| 38 | ; time_std: associated time axis |
|---|
| 39 | ; err_std: error bar on std dev + new mean |
|---|
| 40 | ; } |
|---|
| 41 | ; |
|---|
| 42 | ; usage : field=spectra(datin, nyears, tukey, running, full, bootwin) |
|---|
| 43 | ; then field.std is the standard deviation |
|---|
| 44 | ; field.spectw the normalized spectrum, etc... |
|---|
| 45 | ; |
|---|
| 46 | ;; -------------------------------------------------------------------------- |
|---|
| 47 | ; |
|---|
| 48 | ; 0. Initializations |
|---|
| 49 | ; ------------------ |
|---|
| 50 | |
|---|
| 51 | ntime = nyears*12 ; size of time serie |
|---|
| 52 | lfreq = ntime/2 ; number of wavelenghts |
|---|
| 53 | ntim1=ntime-1 |
|---|
| 54 | |
|---|
| 55 | period = fltarr(lfreq) ; frequency axis |
|---|
| 56 | amean = fltarr(12) ; mean seasonal cycle |
|---|
| 57 | ameansc = fltarr(12) ; mean seasonal cycle (work array) |
|---|
| 58 | datad = fltarr(ntime) ; anomaly data normalized |
|---|
| 59 | anom = fltarr(ntime) ; anomaly data (non-normalized) |
|---|
| 60 | stdsc = fltarr(12) ; Seasonal cycle of STD DEV |
|---|
| 61 | weight = fltarr(ntime+1) |
|---|
| 62 | cx = fltarr(ntim1+1) |
|---|
| 63 | rx = fltarr(ntim1+1) |
|---|
| 64 | powerux = fltarr(lfreq) ; spectra |
|---|
| 65 | powerwx = fltarr(lfreq) ; spectra with Tukey window |
|---|
| 66 | |
|---|
| 67 | run_std = fltarr(ntime-(running*12)+1) ; time serie of STD DEV |
|---|
| 68 | sc_run = fltarr(ntime-(running*12)+1) ; time serie of running SC |
|---|
| 69 | anom_stdsc = fltarr(ntime-(running*12)+1) |
|---|
| 70 | time_std = (findgen(ntime-(running*12)+1)+(running*12)/2)/12. |
|---|
| 71 | |
|---|
| 72 | err_std = fltarr(2) |
|---|
| 73 | |
|---|
| 74 | pi = 3.141592654 |
|---|
| 75 | twopin=2.0*pi/float(ntime) |
|---|
| 76 | pirm=pi/float(ntime) |
|---|
| 77 | |
|---|
| 78 | period = (float(ntime)/(float(findgen(lfreq))+1.))/12. |
|---|
| 79 | |
|---|
| 80 | ; |
|---|
| 81 | ; 1. Compute monthly means and anomalies |
|---|
| 82 | ; -------------------------------------- |
|---|
| 83 | ; |
|---|
| 84 | ; mean of time serie |
|---|
| 85 | ; |
|---|
| 86 | smean = mean(datin) |
|---|
| 87 | ; mean seasonal cycle |
|---|
| 88 | ; |
|---|
| 89 | FOR t = 0, 12-1 DO BEGIN |
|---|
| 90 | amean(t) = mean(datin(long(findgen(ntime/12)*12+t))) |
|---|
| 91 | ENDFOR |
|---|
| 92 | ; |
|---|
| 93 | ; range of mean SC |
|---|
| 94 | ; |
|---|
| 95 | sc_range = max(amean)-min(amean) |
|---|
| 96 | ; |
|---|
| 97 | ; interannual anomaly |
|---|
| 98 | ; |
|---|
| 99 | CASE full OF |
|---|
| 100 | '0': datad = datin - reform(amean#replicate(1, long(ntime/12)), ntime) |
|---|
| 101 | '1': datad = datin |
|---|
| 102 | ENDCASE |
|---|
| 103 | |
|---|
| 104 | anom = datin - reform(amean#replicate(1, long(ntime/12)), ntime) |
|---|
| 105 | |
|---|
| 106 | xmean = mean(datad) |
|---|
| 107 | ; |
|---|
| 108 | ; standard deviation |
|---|
| 109 | ; |
|---|
| 110 | sdev = sqrt((moment(datad))[1]) |
|---|
| 111 | |
|---|
| 112 | ; error bar on std |
|---|
| 113 | |
|---|
| 114 | temperr = stat_error(anom, bootwin) |
|---|
| 115 | ; err= +- sqrt(sum(temperr(i)-sdev)**2)/N) |
|---|
| 116 | err_std(0) = sqrt( total((temperr-sdev)^2)/n_elements(temperr) ) |
|---|
| 117 | ; err_std(0) = sqrt((moment(temperr))[1]) |
|---|
| 118 | err_std(1) = mean(temperr) |
|---|
| 119 | |
|---|
| 120 | ; |
|---|
| 121 | ; standard deviation and SC in running window |
|---|
| 122 | ; |
|---|
| 123 | FOR t = 0, ntime-(running*12) DO BEGIN |
|---|
| 124 | run_std(t) = sqrt((moment(datad[t:t+(running*12)-1]))[1]) |
|---|
| 125 | ; dattmp = datin[t:t+(running*12)-1] |
|---|
| 126 | ; ntsc = running*12 |
|---|
| 127 | ; FOR tsc = 0, 12-1 DO BEGIN |
|---|
| 128 | ; ameansc(tsc) = mean(dattmp(long(findgen(ntsc/12)*12+tsc))) |
|---|
| 129 | ; ENDFOR |
|---|
| 130 | ; sc_run(t) = sqrt((moment(ameansc))[1]) |
|---|
| 131 | ; sc_run(t) = max(ameansc)-min(ameansc) |
|---|
| 132 | ; sc_run(t) = mean(dattmp(long(findgen(running/12)*12+ (t MOD 12)))) |
|---|
| 133 | ENDFOR |
|---|
| 134 | |
|---|
| 135 | ; anomaly time serie wrt time varying SC |
|---|
| 136 | |
|---|
| 137 | ; anom_stdsc = datin - [sc_run(0)#replicate(1, long(running*12/2)), sc_run, sc_run(ntime-(running*12))#replicate(1, long(running*12/2-1))] |
|---|
| 138 | |
|---|
| 139 | ; run_std = smooth(run_std, running*12) |
|---|
| 140 | ; run_std[0:(running*12)/2] = 0. |
|---|
| 141 | ; run_std[ntime-(running*12)/2, ntime-1] = 0. |
|---|
| 142 | ; |
|---|
| 143 | ; standard deviation SC |
|---|
| 144 | ; |
|---|
| 145 | FOR t = 0, 12-1 DO BEGIN |
|---|
| 146 | stdsc(t) = sqrt((moment(datin(long(findgen(ntime/12)*12+t))))[1]) |
|---|
| 147 | ENDFOR |
|---|
| 148 | ; |
|---|
| 149 | ; Normalize data by standard deviation |
|---|
| 150 | ; |
|---|
| 151 | datad = datad-xmean |
|---|
| 152 | datad = datad/sdev |
|---|
| 153 | |
|---|
| 154 | ; |
|---|
| 155 | ; 2. Spectra |
|---|
| 156 | ; ---------- |
|---|
| 157 | ; |
|---|
| 158 | ; Compute autocovariance coefficients |
|---|
| 159 | ; |
|---|
| 160 | FOR k = 0, ntim1 DO BEGIN |
|---|
| 161 | cx(k) = (total(datad(0:ntime-k-1)*(shift(datad, -k))(0:ntime-k-1)))/float(ntime) |
|---|
| 162 | ENDFOR |
|---|
| 163 | ; |
|---|
| 164 | ; Compute autocorrelations for plotting |
|---|
| 165 | ; |
|---|
| 166 | rx = cx/cx(0) |
|---|
| 167 | ; |
|---|
| 168 | ; Compute the spectrum with no windowing |
|---|
| 169 | ; |
|---|
| 170 | FOR i=0,lfreq-2 DO BEGIN |
|---|
| 171 | omega=twopin*float(i+1) |
|---|
| 172 | powerux(i) = (cx(0)+total((2.0*shift(cx, -1)*cos(omega*float(findgen(ntim1)+1)))(0:ntim1-1)))/pi |
|---|
| 173 | ENDFOR |
|---|
| 174 | ; |
|---|
| 175 | ; Compute the spectrum with windowing |
|---|
| 176 | ; |
|---|
| 177 | mm=tukey*12 |
|---|
| 178 | pirm=pi/float(mm) |
|---|
| 179 | ; |
|---|
| 180 | ; compute the weights for the Tukey window |
|---|
| 181 | ; |
|---|
| 182 | weight=0.5*(1.0+cos(pirm*findgen(mm+1))) |
|---|
| 183 | ; |
|---|
| 184 | ; Compute windowed spectrum |
|---|
| 185 | ; |
|---|
| 186 | FOR i=0,lfreq-2 DO BEGIN |
|---|
| 187 | omega=twopin*float(i+1) |
|---|
| 188 | powerwx(i) = (weight(0)*cx(0)+total((2.0*shift(weight, -1)*shift(cx, -1)*cos(omega*float(findgen(mm)+1)))(0:mm-1)))/pi |
|---|
| 189 | ENDFOR |
|---|
| 190 | |
|---|
| 191 | pmax = max(powerwx) |
|---|
| 192 | |
|---|
| 193 | ; |
|---|
| 194 | ; Find max amplitude freq |
|---|
| 195 | ; |
|---|
| 196 | fmax = max(powerwx/pmax) |
|---|
| 197 | index = where (powerwx/pmax EQ fmax) |
|---|
| 198 | fmax = period(index) |
|---|
| 199 | |
|---|
| 200 | ; |
|---|
| 201 | ; 3. Organise output |
|---|
| 202 | ; ------------------ |
|---|
| 203 | ; |
|---|
| 204 | |
|---|
| 205 | spectra = {std:sdev, mean:smean, sc:amean, sc_range:sc_range, sc_std:stdsc, anom:anom, spec:powerux, spectw:powerwx/pmax, time:period, fmax:fmax, run_std:run_std, sc_run:sc_run, anom_stdsc:anom_stdsc, time_std:time_std, err_std:err_std} |
|---|
| 206 | |
|---|
| 207 | return, spectra |
|---|
| 208 | |
|---|
| 209 | END |
|---|
