%% IRASA
% The most fascinating aspect of scale-free structures is what their name 
% suggest: they are scale-free. Let's start with a spatial example:
% Coastline paradox (https://en.wikipedia.org/wiki/Coastline_paradox). You
% can get higher and higher coastline length if you use smaller and smaller
% measurement lengths (see the figure in the link). This is due to the
% fractal-like property of the coastline. When you zoom into the coastline,
% you start to see that the zoomed in version is actually very similar to
% the zoomed-out one. In that sense, the structure of the coastline is
% independent from the spatial scale. 
% 
% Mathematicians were aware of the fractals since as early as Leibniz.
% Weird paradoxes stemming from conceptual structures such as Cantor's set
% (https://en.wikipedia.org/wiki/Cantor_set) was thought to be purely ideal
% ideas, with no applications to real world. In 1963, Lorenz did the first
% rigorous treatment of a chaotic system
% (https://doi.org/10.1175/1520-0469(1963)020%3C0130:DNF%3E2.0.CO;2).
% There, he found what is now called a strange attractor, a fractal that is
% neither a 3 dimensional nor a 2 dimensional surface, but has a
% dimensionality of somewhere in between. In a couple of decades, chaos and
% fractals were popping up everywhere. The work by Feigenbaum and
% Mandelbrot made the study rigorous and showed how so many different
% processes have these universal properties.
% 
% Our topic of interest is 1/f^alpha noise. When the PLE (alpha) is 
% approximately 1, we see a fractal property in time. The process is
% independent of the timescale (it has the same power spectrum in 1 day or
% 1 year or millions of centuries). Moreover, the power spectrum is also
% independent of the measurement length: sampling frequency. 
% 
% The neurophysiological time series are not exactly fractal. Alpha
% oscillations are an ubiquitous property of EEG / MEG PSD. An oscillation
% is not fractal, if you zoom into a peak, the whole thing would look like
% a gaussian instead of an oscillation. Imagine sampling only the peaks of
% an oscillation, the result would be a straight line. This property is
% what the irregular resampling autospectral analysis (IRASA) uses to
% seperate oscillations from the fractal component. (Wen and Liu, 2016;
% https://doi.org/10.1007/s10548-015-0448-0)
% 
% I did not write a function for IRASA since it is a bit complicated and
% already written. I could not find the code mentioned in the paper, the
% link seems broken but in the dynameas toolbox by Soren Wainio-Theberge,
% the IRASA package is in the dependencies
% (https://github.com/SorenWT/dynameas/tree/master/dependencies/measures/IRASA).
% Also, if one wants to do it in the fieldtrip, it is also implemented in
% the fieldtrip toolbox too
% (https://www.fieldtriptoolbox.org/example/irasa/#extracting-spectral-features-from-human-ecog-data).
% I will give a simple example below and talk about the method. However, if
% you are going to use this method, I urge you to check the Wen and Liu
% paper. 
% 
% Briefly, in the IRASA method, the signal is irregularly resampled. For
% example, a resampling of 1.5 means first upsampling the data by a
% factor of 3, then downsampling by a factor of 2. These irregular
% resamples change the position of oscillation in the power spectrum. The
% IRASA method does these enough times, then takes the geometric mean of
% the power spectra. This eliminates the oscillations and the resulting
% part is the fractal / scale-free component. 
% 
% Authored by Yasir Çatal
% catalyasir@gmail.com
%% Example
frange = [1 50];
fs = 100;

data = randn(5000, 5); % 5-channel data (time x channels)
data = data + sin(2*pi*10* linspace(1/fs, 10, 5000))'; % Add 10 hz oscillation


irasa_fractal = amri_sig_fractal(data, fs, 'frange', frange);
out = amri_sig_plawfit(irasa_fractal, frange);
fractal_ple = out.Beta;

% Plot the results for the first channel
loglog(irasa_fractal.freq, irasa_fractal.mixd(:, 1))
hold on
loglog(irasa_fractal.freq, irasa_fractal.frac(:, 1))
loglog(irasa_fractal.freq, irasa_fractal.osci(:, 1))
grid on
legend({'Mixed', 'Fractal','Oscillatory'})
