feat(loadtunnel): [WIP] basic implementation of loadtunnel

loadtunnel.m

+function tunnelData = loadtunnel(tunnelFile, cutoff)
+    % LOADTUNNEL Load the Wind Tunnel data and clean it.
+    %   This function loads the WindTunnel data. Then it:
+    %       - applies a low-pass filter to clean the signal
+    %       - trim the instants before the flapping starts
+    % -----
+    % Output format:
+    %   The output array has one row per timestep. The columns are ordered as follows:
+    %       time, fxF, fyF, fzF, mxF, myF, mzF, fxA, fyA, fzA, mxA, myA, mzA
+    %
+    % Syntax:
+    %   tunnelData = loadtunnel(tunnelFile) load the wind tunnel data file, clean it and return a
+    %   cleaned version with only the relevant portion in form of an array.
+    %
+    % Inputs:
+    %   tunnelFile : Wind tunnel file
+    %   cutoff     : Cutoff frequency for the low-pass filter, [Hz]
+    %
+    % Outputs:
+    %   tunnelData: Array with the cleaned data.
+    %
+    % See also: readdata.
+
+    % ----------------------------------------------------------------------------------------------
+    % TODO: Check for inconsistencies in the table
+    % ----------------------------------------------------------------------------------------------
+    % (c) Copyright 2022 University of Liege
+    % Author: Thomas Lambert <t.lambert@uliege.be>
+    % ULiege - Aeroelasticity and Experimental Aerodynamics
+    % MIT License
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    % Defaults and constants
+    SAMPLING = 1000; % Sampling rate, [Hz]
+
+    % Load data
+    tunnelTable = table;
+    if ~isempty(tunnelFile)
+        tunnelTable = loadcsv(['wtFiles/', tunnelFile, '.txt']);
+    end
+
+    % Convert into table
+    tunnelRaw = table2array(tunnelTable);
+
+    % Filter
+    tunnelFiltered = filterwt(tunnelRaw, cutoff, SAMPLING, 'gaussianfft');
+    tess = fft(tunnelRaw(:, 3));
+    L = size(tess, 1);
+    P2 = abs(tess / L);
+    P1 = P2(1:L / 2 + 1);
+    P1(2:end - 1) = 2 * P1(2:end - 1);
+    f = SAMPLING * (0:(L / 2)) / L;
+    figure;
+    plot(f, P1);
+
+    for i = 1:size(tunnelFiltered, 2)
+        figure;
+        hold on;
+        plot(tunnelRaw(:, i));
+        plot(tunnelFiltered(:, i));
+        legend('Raw', 'Filtered');
+    end
+
+    % Calibrate
+    % Trim
+
+end
+
+function rawTable = loadcsv(file)
+    % LOADCSV Load the csv into a table
+
+    % Tweak options before import
+    opts = detectImportOptions(file);
+    opts.VariableNames = {'Fx_front', 'Fy_front', 'Fz_front', ...
+                          'Mx_front', 'My_front', 'Mz_front', ...
+                          'Fx_aft', 'Fy_aft', 'Fz_aft', ...
+                          'Mx_aft', 'My_aft', 'Mz_aft'};
+    opts.DataLines = [1 Inf];
+
+    % Import
+    rawTable = readtable(file, opts);
+
+end
+
+function calibrated = calibrate(raw, type)
+    % CALIBRATE Calibrate a force or a moment properly
+    %   My creates a parasitic force Fx that must be substracted from the real Fx
+
+    persistent kF kA % Calibration factors for My-Fx parasitic force
+
+    % Proceed to calibration using the calibration data
+    if isempty(kF) || isemepty(kA)
+        kF = 0; % [FIXME] Temporary placeholder
+        kA = 0; % [FIXME] Temporary placeholder
+    end
+
+    % Various calibrations
+    switch type
+        case 'fx'
+    end
+
+    calibrated = raw; % [FIXME] Temporary palceholder
+
+end
+
+function filtered = filterwt(raw, cutoff, sampling, filterFct)
+    % FILTERWT Apply a low-pass filter to the raw data
+
+    GAUSS_SIGMA = 8;
+    GAUSS_RANGE = 20;
+
+    switch filterFct
+        case 'lowpass'
+            filtered = lowpass(raw, cutoff, sampling);
+        case 'gaussianfft'
+            transf = fft(raw);
+            nts = size(raw, 1);
+
+            gauss = zeros(size(raw));
+            gauss(1:GAUSS_RANGE + 1) = exp(-(1:GAUSS_RANGE + 1).^2 / (2 * GAUSS_SIGMA^2));
+            gauss(end - GAUSS_RANGE + 1:end) = fliplr(gauss(2:GAUSS_RANGE + 1));
+            filtered = ifft(transf .* gauss, nts);
+
+        case 'fft'
+
+        otherwise
+            error('loadtunnel:UnknownFilter', ['Filter unknown. Please pick either ''lowpass'''...
+                                               'or ''fft''.']);
+    end
+
+end

processall.m

@@ -38,6 +38,7 @@ function exp = processall(idFile)
     addpath(genpath('.'));
 
     % Defaults and constants
+    CUTOFF_MULTIPLIER = 1.5; % Multiplier for the theorical frequency to use as a cutoff, [-]
 
     % Input checks and processing
     % [filePaths, fileNames] = getfiles(idFile);
@@ -47,16 +48,18 @@ function exp = processall(idFile)
     exp(size(idTable, 1)) = struct;
     for iTest = 1:size(idTable, 1)
         exp(iTest).arduFile = idTable(iTest, :).arduinoFile{:};
+        exp(iTest).tunnelFile = idTable(iTest, :).wtFile{:};
+
         arduData = loadardu(exp(iTest).arduFile);
+        tunnelData = loadtunnel(exp(iTest).tunnelFile, CUTOFF_MULTIPLIER * idTable(iTest, :).freq);
+        %
+        %   [ardData, dataTunnel] = syncdata(ardData, tunnelData);
 
         % Split datasets between front and aft
         size(arduData);
         exp(iTest).arduDataF = arduData(:, [1:3, 6]);
         exp(iTest).arduDataA = arduData(:, [1, 4:5, 7]);
 
-        %   tunnelData = loadtunnel(tunnelFile); (read csv, low pass filter)
-        %
-        %   [ardData, dataTunnel] = syncdata(ardData, tunnelData);
     end
 
     % loop on all rows of idTable
@@ -69,9 +72,11 @@ function exp = processall(idFile)
     interval = 1 / 10;
     figure;
     hold on;
-    for i = 1:4
-        plot(exp(iFile).arduData(1:round(end * interval), 1), ...
-             exp(iFile).arduData(1:round(end * interval), i + 1), 'linewidth', 1.5);
+    for i = 1:2
+        plot(exp(iFile).arduDataF(1:round(end * interval), 1), ...
+             exp(iFile).arduDataF(1:round(end * interval), i + 1), 'linewidth', 1.5);
+        plot(exp(iFile).arduDataA(1:round(end * interval), 1), ...
+             exp(iFile).arduDataA(1:round(end * interval), i + 1), 'linewidth', 1.5);
     end
     legend;