diff --git a/resultsanalysis.m b/resultsanalysis.m
new file mode 100644
index 0000000000000000000000000000000000000000..2173774443fe6c57b542cf0358738d6370209a18
--- /dev/null
+++ b/resultsanalysis.m
@@ -0,0 +1,44 @@
+function ResData = resultsanalysis(matFiles)
+    % RESULTSANALYSIS Analysis of the different results obtained
+    %   Todo
+
+    % ----------------------------------------------------------------------------------------------
+    % (c) Copyright 2022 University of Liege
+    % Author: Thomas Lambert <t.lambert@uliege.be>
+    % ULiege - Aeroelasticity and Experimental Aerodynamics
+    % MIT License
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    close all;
+    clc;
+
+    addpath(genpath('.'));
+
+    % Defaults and constants
+    persistent matFiles_
+    persistent ResData_
+
+    PLOT_RES = false;
+    FORCE_RELOAD = false;
+
+    if (isempty(matFiles_) || ~strcmp(matFiles_, matFiles)) && ~FORCE_RELOAD
+
+        % Load data
+        allFiles = dir(matFiles);
+        for iFile = 1:numel(allFiles)
+            resFile = fullfile(allFiles(iFile).folder, allFiles(iFile).name);
+            tmp = load(resFile);
+            ResData(iFile, :) = tmp.ExpData;
+        end
+        ResData = rmfield(ResData, {'arduStart', 'tunnelStart'});
+
+        % Cache values to prevent excessive reloads
+        matFiles_ = matFiles;
+        ResData_ = ResData;
+    else
+        ResData = ResData_;
+    end
+
+    idx = findwingphase(ResData(1, 1), 0);
+
+end
diff --git a/utils/analysis/findwingphase.m b/utils/analysis/findwingphase.m
new file mode 100644
index 0000000000000000000000000000000000000000..121631cf7f30eb843ce24a84e0f50657dd5354ef
--- /dev/null
+++ b/utils/analysis/findwingphase.m
@@ -0,0 +1,60 @@
+function idx = findwingphase(ResData, phase)
+    % FINDWINDPHASE Return the indexes where the wings are at a certain phase
+    %   Todo
+
+    % ----------------------------------------------------------------------------------------------
+    % (c) Copyright 2022 University of Liege
+    % Author: Thomas Lambert <t.lambert@uliege.be>
+    % ULiege - Aeroelasticity and Experimental Aerodynamics
+    % MIT License
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    % Find flapping frequencies and period from data
+    trueFreq(1) = getmainfrequency(ResData.front.arduShort(:, 1), ResData.front.arduShort(:, 2));
+    trueFreq(2) = getmainfrequency(ResData.aft.arduShort(:, 1), ResData.aft.arduShort(:, 2));
+
+    % Determine period offset that match phase
+
+    % Find peaks in the wing angles
+
+    % Find indexes with proper phase difference
+
+    idx = 0; % FIXME, proper return
+
+end
+
+function trueFreq = getmainfrequency(time, signal)
+    % GETFREQUENCIES Return the main signal frequencies
+
+    % Constants and defaults
+    DEBUG = false;
+
+    % Signal parameters
+    dt = diff(time);
+    sampling = 1 / dt(1);
+    if mod(length(signal), 2) ~= 0
+        signal = signal(1:end - 1);
+    end
+    nts = length(signal);
+
+    % Fourier transform
+    transform = fft(signal);
+
+    fT2 = abs(transform / nts); % Two-sided spectrum
+    fT1 = fT2(1:nts / 2 + 1); % Single-sided spectrum
+    fT1(2:end - 1) = 2 * fT1(2:end - 1);
+    freq = sampling * (0:(nts / 2)) / nts; % Frequency domain
+
+    % Get main frequency
+    [~, trueFreqIdx] = max(fT1);
+    trueFreq = freq(trueFreqIdx);
+
+    if DEBUG
+        figure;
+        plot(freq, fT1);
+        title("Single-Sided Amplitude Spectrum of S(t)");
+        xlabel("f (Hz)");
+        ylabel("|P1(f)|");
+    end
+
+end