---
tags: Showcase, Notes
geometry: margin=3cm,top=3cm
---
# Power Spectrum of Response Functions
> And the shape of the temperature decay!
[](https://hackmd.io/@engeir/SyCB0N-_i)
[](https://github.com/engeir/hack-md-notes/blob/main/power-spectrum.md)
## Temperature decay
We first do a quick comparison between the decaying phase of the temperature time
series and a couple of analytical functions.
All fits are made with the [**lmfit**](https://lmfit.github.io/lmfit-py/) library.
As one can see, the single exponential does not do a good job, while the double
exponential and the power law both follow the decay of the original time series well.
From the statistics provided by **lmfit**, the double exponential fit comes out on top
concerning the `R-squared` value, with a value $0.006222$ higher than what the power law
fit obtained. The fitting parameters are four for the double exponential, and two for
the power law. One should note, however, that the time axis has also been adjusted to
improve the fit of the power law.
### Exponential(s)
In the plots below, we try both with one single exponential and with two superposed
exponentials.
The single exponential was made using the following fitting parameters:
```toml
[[Model]]
Model(exponential)
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 32
# data points = 7200
# variables = 2
chi-square = 0.69728364
reduced chi-square = 9.6872e-05
Akaike info crit = -66541.2750
Bayesian info crit = -66527.5114
R-squared = 0.90800218
[[Variables]]
amplitude: 0.25565769 +/- 0.00155324 (0.61%) (init = 0.1137494)
decay: 4.77277534 +/- 0.02355057 (0.49%) (init = 8.94272)
[[Correlations]] (unreported correlations are < 0.100)
C(amplitude, decay) = -0.915
```
The double exponential was made using the following fitting parameters:
```toml
[[Model]]
(Model(exponential, prefix='exp1_') + Model(exponential, prefix='exp2_'))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 66
# data points = 7200
# variables = 4
chi-square = 0.12422683
reduced chi-square = 1.7263e-05
Akaike info crit = -78957.8737
Bayesian info crit = -78930.3463
R-squared = 0.98360983
[[Variables]]
exp1_amplitude: 0.03284301 +/- 7.1478e-04 (2.18%) (init = 1)
exp1_decay: 25.2601488 +/- 0.80882239 (3.20%) (init = 80)
exp2_amplitude: 0.49291210 +/- 0.00397228 (0.81%) (init = 2)
exp2_decay: 2.30849841 +/- 0.01661270 (0.72%) (init = 2)
[[Correlations]] (unreported correlations are < 0.100)
C(exp1_amplitude, exp1_decay) = -0.983
C(exp2_amplitude, exp2_decay) = -0.923
C(exp1_amplitude, exp2_decay) = -0.884
C(exp1_decay, exp2_decay) = 0.830
C(exp1_amplitude, exp2_amplitude) = 0.661
C(exp1_decay, exp2_amplitude) = -0.598
```
### Power Law
The power law was made using the following fitting parameters:
```toml
[[Model]]
Model(powerlaw)
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 13
# data points = 7200
# variables = 2
chi-square = 0.17138784
reduced chi-square = 2.3810e-05
Akaike info crit = -76644.7702
Bayesian info crit = -76631.0065
R-squared = 0.97738753
[[Variables]]
amplitude: 0.77438211 +/- 0.00346024 (0.45%) (init = 0.5914298)
exponent: -1.39394516 +/- 0.00268833 (0.19%) (init = -1.269584)
[[Correlations]] (unreported correlations are < 0.100)
C(amplitude, exponent) = -0.964
```
## Power Spectrum of the Response function
We first define how PSD is to be calculated. We do this with the Welch method from the
**scipy** library:
```python
def spectrum(signal) -> tuple[np.ndarray, np.ndarray]:
"""Calculate the one sided spectrum of the signal."""
signal = (signal - signal.mean()) / signal.std()
sample_frequency = 365
frequency, power = ssi.welch(
signal, sample_frequency, nperseg=2**13, return_onesided=False
)
frequency_plus = frequency[frequency > 0]
power_plus = power[frequency > 0]
return np.asarray(frequency_plus[1:]), np.asarray(power_plus[1:])
```
The power spectrum of the response function obtained from the CESM LME data set follows
a power law, which can be seen when plotted on the log-log axis. In this case, we see
the weakness of the [**lmfit**](https://lmfit.github.io/lmfit-py/) library (red line in
the plot below), since it does not properly consider data as being on a log-log scale.
That is, the absolute difference on a linear axis is used, which for the small values in
the tail makes for a bad fit when plotted on log-log.
Let us see if the [**powerlaw**](https://github.com/jeffalstott/powerlaw) library can
cook up something better. This is the yellow fit in the figure below, which does come
closer to the original power spectrum (blue). This library returns an exponent of
`1.824` (the library defines the power law as $p(x)\propto x^{-\alpha}$) while **lmfit**
returns an exponent of `-1.306`.
Also included is a double exponential fit, again using the **lmfit** library, which in
this case does not provide a good fit over the full frequency range.
**lmfit** output from power law fit of the power spectra of the response function.
```toml
[[Model]]
Model(powerlaw)
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 35
# data points = 896
# variables = 2
chi-square = 0.00280794
reduced chi-square = 3.1409e-06
Akaike info crit = -11351.2270
Bayesian info crit = -11341.6311
R-squared = 0.90466032
[[Variables]]
amplitude: 0.00518431 +/- 1.9615e-04 (3.78%) (init = 0.001236972)
exponent: -1.30151327 +/- 0.01797807 (1.38%) (init = -2.078565)
[[Correlations]] (unreported correlations are < 0.100)
C(amplitude, exponent) = 0.958
```
**powerlaw** output from power law fit of the power spectra of the response function,
where `sigma` is the standard deviation on the power `alpha`.
```toml
fit.power_law.alpha = 1.8239139860706672
fit.power_law.sigma = 0.5826287465799062
```
### All PSDs
```toml
[[Model]]
Model(powerlaw)
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 13
# data points = 889
# variables = 2
chi-square = 5.5977e-20
reduced chi-square = 6.3108e-23
Akaike info crit = -45441.1908
Bayesian info crit = -45431.6106
R-squared = 1.00000000
[[Variables]]
amplitude: 0.00109805 +/- 1.5096e-05 (1.37%) (init = 0.001195273)
exponent: -2.03713872 +/- 0.00465884 (0.23%) (init = -2.066896)
[[Correlations]] (unreported correlations are < 0.100)
C(amplitude, exponent) = -0.970
#######################################################################
[[Model]]
Model(powerlaw)
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 22
# data points = 889
# variables = 2
chi-square = 3.0428e-19
reduced chi-square = 3.4305e-22
Akaike info crit = -43936.0946
Bayesian info crit = -43926.5144
R-squared = 1.00000000
[[Variables]]
amplitude: 0.00847024 +/- 3.4066e-04 (4.02%) (init = 0.005694137)
exponent: -2.47469249 +/- 0.01655384 (0.67%) (init = -2.380111)
[[Correlations]] (unreported correlations are < 0.100)
C(amplitude, exponent) = -0.917
#######################################################################
[[Model]]
Model(powerlaw)
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 26
# data points = 889
# variables = 2
chi-square = 2.1253e-17
reduced chi-square = 2.3961e-20
Akaike info crit = -40161.1313
Bayesian info crit = -40151.5511
R-squared = 1.00000000
[[Variables]]
amplitude: 0.00544160 +/- 4.3852e-04 (8.06%) (init = 0.002388811)
exponent: -1.40405019 +/- 0.02461607 (1.75%) (init = -1.153645)
[[Correlations]] (unreported correlations are < 0.100)
C(amplitude, exponent) = -0.991
```
Sign in with Wallet
Connect another wallet