diff --git a/processall.m b/processall.m
index ea1b702c16c6d6e1733dc794a13dc549a1cc80af..987462d781a90361c7b67aa657253aa825d3a732 100644
--- a/processall.m
+++ b/processall.m
@@ -69,7 +69,7 @@ function ExpData = processall(idFile)
 
         % TODO
 
-        % Split datasets between front and aft
+        % Split datasets between front and aft [FIXME]: Get rid of old deprecated variables
         ExpData(iTest).front.arduFull = arduFull(:, [1:3, 6]);
         ExpData(iTest).aft.arduFull = arduFull(:, [1, 4:5, 7]);
         ExpData(iTest).front.arduShort = arduShort(:, [1:2, 4]);
@@ -77,6 +77,24 @@ function ExpData = processall(idFile)
         ExpData(iTest).front.tunnelData = tunnelData(:, 1:6);
         ExpData(iTest).aft.tunnelData = tunnelData(:, 7:end);
 
+        % Better variables
+        ExpData(iTest).time = arduShort(:, 1);
+        ExpData(iTest).Front.angles = arduShort(:, 2);
+        ExpData(iTest).Front.forces = tunnelData(:, 1:3);
+        ExpData(iTest).Front.moments = tunnelData(:, 4:6);
+        ExpData(iTest).Front.power = arduShort(:, 4);
+
+        ExpData(iTest).Aft.angles = arduShort(:, 3);
+        ExpData(iTest).Aft.forces = tunnelData(:, 6 + (1:3));
+        ExpData(iTest).Aft.moments = tunnelData(:, 6 + (4:6));
+        ExpData(iTest).Aft.power = arduShort(:, 5);
+
+        ExpData(iTest).Testcase.dx = idTable(iTest, :).dstX;
+        ExpData(iTest).Testcase.airspeed = idTable(iTest, :).v_air;
+        ExpData(iTest).Testcase.freq = idTable(iTest, :).freq;
+        ExpData(iTest).Testcase.temp = idTable(iTest, :).temp_air;
+        ExpData(iTest).Testcase.press = idTable(iTest, :).pres_air;
+
         if PLOT_CLEANED
             plottime(ExpData(iTest), 'arduShort', false);
         end
diff --git a/utils/analysis/findwingphase.m b/utils/analysis/findwingphase.m
index c8cb3d953f8e3d0ef77dd0c222a39b66c74bfc39..58023f681a62374cd040c5d6d55aed876777a63e 100644
--- a/utils/analysis/findwingphase.m
+++ b/utils/analysis/findwingphase.m
@@ -13,20 +13,20 @@ function peakRanges = findwingphase(ResData, phaseDeg)
     N_PHASE_PEAKS = 5; % Number of peaks in phase to return ideally
 
     % Abbreviations
-    sampling = 1 / diff(ResData.front.arduShort(1:2, 1));
+    sampling = 1 / diff(ResData.time(1:2));
 
     % 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));
+    trueFreq(1) = getmainfrequency(ResData.time, ResData.Front.angles);
+    trueFreq(2) = getmainfrequency(ResData.time, ResData.Aft.angles);
     period = 1 ./ [trueFreq(1), trueFreq(1)]; % Signal perdio, [s]
     nPeriodIdx = period * sampling; % Number of indexes in a period
 
     fprintf('    TrueFreq = %0.2f/%0.2fHz, theory = %0.2fHz\n', ...
-            trueFreq(1), trueFreq(2), ResData.freq);
+            trueFreq(1), trueFreq(2), ResData.Testcase.freq);
 
     % Determine the ranges when the two angles have the proper phase difference
-    peakRanges = findsyncedpeaks(ResData.front.arduShort(:, 2), ...
-                                 ResData.aft.arduShort(:, 2), ...
+    peakRanges = findsyncedpeaks(ResData.Front.angles, ...
+                                 ResData.Aft.angles, ...
                                  nPeriodIdx(1), phaseDeg, N_PHASE_PEAKS);
 
 end
@@ -83,6 +83,7 @@ end
 
 function signalRange = findsyncedpeaks(signal1, signal2, nPeriodIdx, offset, nPeaks)
     % FINDANGLEPEAKS Find peaks in the signals that are synced with respect to the offset desired
+    %   Returns signalRange : index of peak before and peak after we are sync
 
     % Constants and defaults
     DEBUG = false;
@@ -137,6 +138,7 @@ function signalRange = findsyncedpeaks(signal1, signal2, nPeriodIdx, offset, nPe
         hold off;
         grid on;
         legend('Front', 'Aft');
+        pause;
     end
 
 end
diff --git a/utils/analysis/forcetocoeff.m b/utils/analysis/forcetocoeff.m
new file mode 100644
index 0000000000000000000000000000000000000000..3d0148aadbd764612b4fdc4b3230b8616dd7e8ed
--- /dev/null
+++ b/utils/analysis/forcetocoeff.m
@@ -0,0 +1,23 @@
+function coeff = forcetocoeff(var, dens, speed, chord, span, type)
+    % CALCCOEFF Calculate force/moment coefficient
+    %   Todo
+
+    % ----------------------------------------------------------------------------------------------
+    % (c) Copyright 2022 University of Liege
+    % Author: Thomas Lambert <t.lambert@uliege.be>
+    % ULiege - Aeroelasticity and Experimental Aerodynamics
+    % MIT License
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    area = chord * span;
+
+    switch type
+        case 'force'
+            nondim = 1 / 2 * dens * area * speed^2;
+        case 'moment'
+            nondim = 1 / 2 * dens * area * speed^2 * chord;
+    end
+
+    coeff = var / nondim;
+
+end
diff --git a/utils/analysis/getboundingpeaks.m b/utils/analysis/getboundingpeaks.m
new file mode 100644
index 0000000000000000000000000000000000000000..4f5bf6a3484898c3681bd4bb6b45ae794e5bfdbc
--- /dev/null
+++ b/utils/analysis/getboundingpeaks.m
@@ -0,0 +1,22 @@
+function bounds = getboundingpeaks(signal, sampling)
+    % GETANGLESIGNALBOUNDS Return the bounds of the angle signal (correspond to first and last peak)
+    %   Todo
+
+    % ----------------------------------------------------------------------------------------------
+    % (c) Copyright 2022 University of Liege
+    % Author: Thomas Lambert <t.lambert@uliege.be>
+    % ULiege - Aeroelasticity and Experimental Aerodynamics
+    % MIT License
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    % Signal analysis
+    freq = getmainfrequency(signal, sampling);
+
+    % Get first and last peaks for wind-off conditions
+    [~, peakLocs] = findpeaks(signal, ...
+                              'MinPeakProminence', max(signal) / 2, ...
+                              'MinPeakDistance', ((1 / freq) * sampling) / 2);
+
+    bounds = [peakLocs(1), peakLocs(end)];
+
+end
diff --git a/utils/analysis/getmainfrequency.m b/utils/analysis/getmainfrequency.m
new file mode 100644
index 0000000000000000000000000000000000000000..f0d0d9fd578072f3da8bc7474fcf57d14226ea45
--- /dev/null
+++ b/utils/analysis/getmainfrequency.m
@@ -0,0 +1,30 @@
+function mainFreq = getmainfrequency(signal, sampling)
+    % GETMAINFREQUENCY Return the main signal frequency
+    %   Todo
+
+    % ----------------------------------------------------------------------------------------------
+    % (c) Copyright 2022 University of Liege
+    % Author: Thomas Lambert <t.lambert@uliege.be>
+    % ULiege - Aeroelasticity and Experimental Aerodynamics
+    % MIT License
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    % Signal parameters
+    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(2:end)); % Make sure it is not zero
+    mainFreq = freq(1 + trueFreqIdx);
+
+end
diff --git a/utils/analysis/removeinertia.m b/utils/analysis/removeinertia.m
new file mode 100644
index 0000000000000000000000000000000000000000..2d474acf1401d2cd593824a31d4c7608dbf06de7
--- /dev/null
+++ b/utils/analysis/removeinertia.m
@@ -0,0 +1,44 @@
+function ResData = removeinertia(ResData, ResDataWO, field)
+    % REMOVEINERTIA Remove intertia forces from ResData
+    %   Todo
+
+    % ----------------------------------------------------------------------------------------------
+    % (c) Copyright 2022 University of Liege
+    % Author: Thomas Lambert <t.lambert@uliege.be>
+    % ULiege - Aeroelasticity and Experimental Aerodynamics
+    % MIT License
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    % Calculate mean inertial forces and moments
+    [ResData.(field).inertiaF, ResData.(field).inertiaM] = inertialforces(ResDataWO, field);
+
+    % -------------------------------------------
+    % For all peaks found in phase, calculate mean forces and moments
+    for i = 1:length(ResData.Phase)
+        if ~isempty(ResData.Phase(i).idx)
+            for j = 1:size(ResData.Phase(i).idx, 1)
+                syncedRange = ResData.Phase(i).idx(1):ResData.Phase(i).idx(2);
+                meanF(j, :) = mean(ResDataWO.(field).forces(syncedRange));
+                meanM(j, :) = mean(ResDataWO.(field).moments(syncedRange));
+            end
+            ResData.Phase(i).(field).rawF = mean(meanF);
+            ResData.Phase(i).(field).rawM = mean(meanM);
+            ResData.Phase(i).(field).forces = mean(meanF) - ResData.(field).inertiaF;
+            ResData.Phase(i).(field).moments = mean(meanM) - ResData.(field).inertiaM;
+        end
+    end
+
+end
+
+function [meanForces, meanMom] = inertialforces(ResDataWO, field)
+
+    % Get first and last peaks for wind-off conditions
+    woSampling = 1 / diff(ResDataWO.time(1:2));
+    woBounds = getboundingpeaks(ResDataWO.(field).angles, woSampling);
+    woIdx = woBounds(1):woBounds(2);
+
+    % Calculate mean force and moments due to inertia
+    meanForces = mean(ResDataWO.(field).forces(woIdx));
+    meanMom = mean(ResDataWO.(field).moments(woIdx));
+
+end