|
| 1 | +""" Helper functions for constructing design regressors |
| 2 | +""" |
| 3 | +from __future__ import division |
| 4 | + |
| 5 | +import numpy as np |
| 6 | + |
| 7 | + |
| 8 | +def dct_ii_basis(volume_times, order=None, normcols=False): |
| 9 | + """ DCT II basis up to order `order` |
| 10 | +
|
| 11 | + See: https://en.wikipedia.org/wiki/Discrete_cosine_transform#DCT-II |
| 12 | +
|
| 13 | + By default, basis not normalized to length 1, and therefore, basis is not |
| 14 | + orthogonal. Normalize basis with `normcols` keyword argument. |
| 15 | +
|
| 16 | + Parameters |
| 17 | + ---------- |
| 18 | + volume_times : array-like |
| 19 | + Times of acquisition of each volume. Must be regular and continuous |
| 20 | + otherwise we raise an error. |
| 21 | + order : None or int, optional |
| 22 | + Order of DCT-II basis. If None, return full basis set. |
| 23 | + normcols : bool, optional |
| 24 | + If True, normalize columns to length 1, so return orthogonal |
| 25 | + `dct_basis`. |
| 26 | +
|
| 27 | + Returns |
| 28 | + ------- |
| 29 | + dct_basis : array |
| 30 | + Shape ``(len(volume_times), order)`` array with DCT-II basis up to |
| 31 | + order `order`. |
| 32 | +
|
| 33 | + Raises |
| 34 | + ------ |
| 35 | + ValueError |
| 36 | + If difference between successive `volume_times` values is not constant |
| 37 | + over the 1D array. |
| 38 | + """ |
| 39 | + N = len(volume_times) |
| 40 | + if order is None: |
| 41 | + order = N |
| 42 | + if not np.allclose(np.diff(np.diff(volume_times)), 0): |
| 43 | + raise ValueError("DCT basis assumes continuous regular sampling") |
| 44 | + n = np.arange(N) |
| 45 | + cycle = np.pi * (n + 0.5) / N |
| 46 | + dct_basis = np.zeros((N, order)) |
| 47 | + for k in range(0, order): |
| 48 | + dct_basis[:, k] = np.cos(cycle * k) |
| 49 | + if normcols: # Set column lengths to 1 |
| 50 | + lengths = np.ones(order) * np.sqrt(N / 2.) |
| 51 | + lengths[0:1] = np.sqrt(N) # Allow order=0 |
| 52 | + dct_basis /= lengths |
| 53 | + return dct_basis |
| 54 | + |
| 55 | + |
| 56 | +def dct_ii_cut_basis(volume_times, cut_period): |
| 57 | + """DCT-II regressors with periods >= `cut_period` |
| 58 | +
|
| 59 | + See: http://en.wikipedia.org/wiki/Discrete_cosine_transform#DCT-II |
| 60 | +
|
| 61 | + Parameters |
| 62 | + ---------- |
| 63 | + volume_times : array-like |
| 64 | + Times of acquisition of each volume. Must be regular and continuous |
| 65 | + otherwise we raise an error. |
| 66 | + cut_period: float |
| 67 | + Cut period (wavelength) of the low-pass filter (in time units). |
| 68 | +
|
| 69 | + Returns |
| 70 | + ------- |
| 71 | + cdrift: array shape (n_scans, n_drifts) |
| 72 | + DCT-II drifts plus a constant regressor in the final column. Constant |
| 73 | + regressor always present, regardless of `cut_period`. |
| 74 | + """ |
| 75 | + N = len(volume_times) |
| 76 | + hfcut = 1./ cut_period |
| 77 | + dt = volume_times[1] - volume_times[0] |
| 78 | + # Such that hfcut = 1/(2*dt) yields N |
| 79 | + order = int(np.floor(2 * N * hfcut * dt)) |
| 80 | + # Always return constant column |
| 81 | + if order == 0: |
| 82 | + return np.ones((N, 1)) |
| 83 | + basis = np.ones((N, order)) |
| 84 | + basis[:, :-1] = dct_ii_basis(volume_times, order, normcols=True)[:, 1:] |
| 85 | + return basis |
0 commit comments