The IAAFT algorithm
This page explains how the IAAFT algorithm works. Cloud scientists probably would like to visit the introduction on surrogate fields first. The figure below depicts the scheme of the IAAFT algorithm, and illustrates it with an Liquid Water Path (LWP) measurement as template. From this template time series (panel 1) the power spectrum is calculated and a sorted list is made of all values in the time series to be used in the amplitude conversions. For theoretical studies, the power spectrum and the values of the amplitude distribution can be predefined.
The algorithm starts with a random shuffle of the data points (panel 2) or with any other white noise time series. In each iteration, the Fourier spectrum is adjusted first (panel 3) and then the amplitudes (panel 4). To produce the desired power spectrum, at each iteration the Fourier transform of the iterated time series is calculated and its squared coefficients are replaced by those of the original time series. The phases are kept unaltered.
In Matlab code this reads:
spectrum = ifft(surrogate); % surrogate is the surrogate time series.
phase = angle(spectrum); % Angle is the Matlab function to calculate the phase from a complex number.
spectrum = fourierCoeff .* exp(i*phase); % fourierCoeff are the magnitudes of the Fourier coefficients of the template.
surrogate = fft(spectrum);
After this step the amplitudes of the iterated time series will no longer be the same. Therefore, in the second iterative step the amplitudes are adjusted by ranking their values and replacing them by the values of the original sorted list having the same ranking. For example, you can see in the figure that the maximum value of the time series in panel 3 is replaced by the maximum value measured.
In Matlab code this reads:
[dummy, index] = sort(surrogate); % We need only the indices. The first value of index points to the highest value, etc.
surrogate(index) = sortedValues; % sortedValues is the vector with the sorted values from the template.
Both iterative steps are repeated until a convergence threshold is reached (panel 5).
As a small variation on the above mentioned IAAFT algorithm, I have developed the Stochastic IAAFT (SIAAFT) algorithm. The difference is that the amplitude adaptation step in no longer done for all values, but only for a fraction of all values. Typically only 1 in 5 values of the time series are set to the value they should have according to their ranking. This way the algorithm becomes less greedy, the convergence rate becomes smaller, and the algorithm is better able to avoid local minima.
In Matlab code this reads:
noValues = numel(surrogate);
[dummy, index] = sort(surrogate); % We need only the indices. The first value of index points to the highest value, etc.
[dummy, index2index] = sort(rand(noValues )); % Create random numbers and sort them to get unique random indices.
index2index = index2index(1:floor(noValues * fraction)); % A random fraction of these unique indices is used.
surrogate(index(index2index)) = sortedValues(index2index); % Without the index2index all values would be adapted, now only a random fraction of the values.
victor.venema@uni-bonn.de
URL: http://www.meteo.uni-bonn.de/mitarbeiter/venema/themes/surrogates/iaaft/iaaft_algorithm.html
This document last modified on: Wednesday May 15, 2013.
