[WIP] Breaking Change - Use Update

Folder name is changed, rework the html templates
Change the organisation.

spindle/Macros_ttt_spindle/Dual_Spindle_error.m

@@ -0,0 +1,152 @@
+function [Res] = Dual_Spindle_error(dataX, dataY, NbTurn,texte, path)
+    L= length(dataX);
+    res_per_rev = L/NbTurn;
+    P = 0: (res_per_rev * NbTurn)-1;
+    Pos =  P' * 360/res_per_rev;
+    Theta = degtorad(Pos)';
+
+    x1 = myfit2(Pos, dataX);
+    y1 = myfit2(Pos, dataY);
+
+    %Convert data to frequency domain and scale accordingly
+    X2 = 2/(res_per_rev*NbTurn)*fft(x1);
+    f2 = (0:L-1)./NbTurn;
+    Y2 = 2/(res_per_rev*NbTurn)*fft(y1);
+
+
+    % Separate the fft integers and not-integers
+    for i = 1 : length(f2)
+        if mod(f2(i), 1) == 0;
+            X2dec(i)= 0;
+            X2int(i)= X2(i);
+            Y2dec(i)= 0;
+            Y2int(i)= Y2(i);
+        else
+            X2dec(i)= X2(i);
+            X2int(i)= 0;
+            Y2dec(i)= Y2(i);
+            Y2int(i)= 0;
+        end
+    end
+
+    if mod(length(f2),2) ==1; % Case length(f2) is odd -> the mirror image of the FFT is reflected between 2 harmonique
+        for i = length(f2)/2 +1.5: length(f2)
+            if mod(f2(i-1), 1) == 0;
+                X2dec(i)= 0;
+                X2int(i)= X2(i);
+                Y2dec(i)= 0;
+                Y2int(i)= Y2(i);
+            else
+                X2dec(i)= X2(i);
+                X2int(i)= 0;
+                Y2dec(i)= Y2(i);
+                Y2int(i)= 0;
+            end
+        end
+    else % Case length(f2) is even -> the mirror image of the FFT is reflected at the Nyquist frequency
+        for i = length(f2)/2 +1: length(f2)
+            if mod(f2(i), 1) == 0;
+                X2dec(i)= 0;
+                X2int(i)= X2(i);
+                Y2dec(i)= 0;
+                Y2int(i)= Y2(i);
+            else
+                X2dec(i)= X2(i);
+                X2int(i)= 0;
+                Y2dec(i)= Y2(i);
+                Y2int(i)= 0;
+            end
+        end
+    end
+
+    X2int(1) = 0; %remove the data average/dc component
+    X2int(NbTurn+1) = 0; %Remove fondamental/eccentricity
+    %  X2int(length(f2)) = 0; %remove the data average/dc component
+    X2int(length(f2)-NbTurn+1) = 0; %Remove eccentricity
+
+    Y2int(1) = 0; %remove the data average/dc component
+    Y2int(NbTurn+1) = 0; %Remove fondamental/eccentricity
+    %  Y2int(length(f2)) = 0; %remove the data average/dc component
+    Y2int(length(f2)-NbTurn+1) = 0; %Remove eccentricity
+
+
+    % Extract the fondamentale-> exentricity
+    for i = 1:length(f2)
+        if i == NbTurn+1 || i==  length(f2)-NbTurn + 1
+            X2fond(i) = X2(i);
+            Y2fond(i) = Y2(i);
+        else
+            X2fond(i) = 0;
+            Y2fond(i) = 0;
+        end
+    end
+
+    X2tot= X2int + X2dec;
+    Y2tot= Y2int + Y2dec;
+
+    %Convert data to "time" domain and scale accordingly
+    Wxint = real((res_per_rev*NbTurn)/2 * ifft(X2int)) ;
+    Wxdec = real((res_per_rev*NbTurn)/2 * ifft(X2dec)) ;
+    Wxtot = real((res_per_rev*NbTurn)/2 * ifft(X2tot)) ;
+
+    %Convert data to "time" domain and scale accordingly
+    Wyint = real((res_per_rev*NbTurn)/2 * ifft(Y2int)) ;
+    Wydec = real((res_per_rev*NbTurn)/2 * ifft(Y2dec)) ;
+    Wytot = real((res_per_rev*NbTurn)/2 * ifft(Y2tot)) ;
+
+
+    %%
+    fig = figure();
+
+    % subplot(3, 2, 5); bar(f2(1:50*NbTurn),abs(W2int(1:50*NbTurn)),3) ;
+    % axis([0,50,0,max(abs(X2int(1:50*NbTurn)))]);
+    % title ('Fourier integer'  ); xlabel('UPR'); ylabel ('Microns')
+
+    % subplot(3, 2, 6); bar(f2(1:50*NbTurn),abs(W2dec(1:50*NbTurn)),2); title (' Fourier non-integer'  );
+    % axis([0,50,0,max(abs(X2dec(1:50*NbTurn)))]);
+    % title ('Fourier non-integer'  ); xlabel('UPR'); ylabel ('Microns')
+
+    % Sensitive pos dir (synchrone+asynchrone)
+    Wtot= Wytot.*cos(Theta)+Wxtot.*sin(Theta);
+    % Wtot= Wytot+Wxtot; = WTot = real((res_per_rev*NbTurn)/2 * ifft(W2tot)) ;
+
+    %Sensitive pos dir (synchrone+asynchrone)
+    Wint= Wyint.*cos(Theta)+Wxint.*sin(Theta);
+    %Sensitive pos dir (synchrone+asynchrone)
+    Wdec= Wydec.*cos(Theta)+Wxdec.*sin(Theta);
+
+
+    % total error motion
+    Total_Error = max(Wtot)- min(Wtot);
+
+    %lsc X synchronous
+    Synchronous_Error = max(Wint)- min(Wint);
+
+    %lsc X Asynchronous
+    var=reshape(Wxdec,length(Wxdec)/NbTurn,NbTurn);
+    for i=1:length(Wxdec)/NbTurn
+        Asynch(i) = max(var(i,:)) - min(var(i,:)) ;
+    end
+    Asynchronous_Error = max(Asynch)- min(Asynch);
+
+
+    % Raw Error Motion without Exentricity (sync +asynch)
+    subplot(2, 2, 2) ; polar2(Theta,Wtot, 'b'); title (' Total error'  );
+    % Residual Synchronous Error Motion without Exentricity (ie fondamental sync err motion)
+    subplot(2, 2, 3) ;polar2(Theta,Wint,'b'); title ( 'Residual synchronous error' );
+    % Asynchronous Error Motion
+    subplot(2, 2, 4) ;polar2(Theta,Wdec, 'b'); title ('Asynchronous error'  );
+
+
+    strmin1 = ['Total error = ', num2str(Total_Error*1000), ' nm'];
+    strmin2 = ['Residual synchronous error = ', num2str(Synchronous_Error*1000), ' nm' ];
+    strmin3 = ['Asynchronous error = ', num2str(Asynchronous_Error*1000), ' nm'];
+    dim0=[0.04 0.5 0.3 .3];%x y w h basgauche to hautdroite
+    dim1=[0.15 0.65 0.3 .3];
+    annotation('textbox',dim0, 'String',{ strmin1 , strmin2, strmin3}, 'FitBoxToText', 'on')
+    annotation('textbox',dim1, 'String',texte, 'FitBoxToText', 'on')
+
+    saveas(fig,fullfile(path,char(texte)),'jpg');
+    Res=1;
+    close all;
+end

spindle/Macros_ttt_spindle/Spindle.m

@@ -0,0 +1,32 @@
+clear; close all; clc;
+
+%%
+datapath = './data/';
+savepath = './Macros_ttt_spindle/';
+
+name_of_the_datafile = '10turns_60rpm_IcepapFIR.txt';
+NbTurn = 10;    % Be careful that this value match with the measure
+
+%%
+% [y, x, z, y2, x2] = importfile(fullfile(datapath,name_of_the_datafile),12, inf);
+
+spindle_table  = readtable([datapath, name_of_the_datafile]);
+
+data  = table2array(spindle_table);
+
+% data = [y x z y2 x2];
+data = data(:, 3:end);
+texte = {'Y Axis' 'X Axis' 'Z Axis' 'Y2 Axis' 'X2 Axis'};
+texte2 = {'Rotating radial Error XY' 'Rotating radial Error X2Y2' 'Tilt Error XX2' 'Tilt Error YY2' };
+
+%%
+for i = 1:5
+  Spindle_error(data(:,i), NbTurn, texte(:,i), savepath);
+end
+
+%%
+Dual_Spindle_error(data(:,2), data(:,1),NbTurn,texte2(:,1), savepath);
+Dual_Spindle_error(data(:,5), data(:,4),NbTurn,texte2(:,2), savepath);
+
+Tilt_Spindle_error(data(:,2), data(:,5),NbTurn,texte2(:,3), savepath);
+Tilt_Spindle_error(data(:,1), data(:,4),NbTurn,texte2(:,4), savepath);

spindle/Macros_ttt_spindle/Spindle_error.m

@@ -0,0 +1,131 @@
+function [Res] = Spindle_error(data, NbTurn, texte, path)
+    %%
+    L = length(data);
+    res_per_rev = L/NbTurn;
+
+    P = 0:(res_per_rev*NbTurn-1);
+    Pos = P' * 360/res_per_rev;
+    Theta = deg2rad(Pos)';
+
+    % Temperature correction
+    x1 = myfit2(Pos, data);
+
+    % Convert data to frequency domain and scale accordingly
+    X2 = 2/(res_per_rev*NbTurn)*fft(x1);
+    f2 = (0:L-1)./NbTurn; %upr -> once per revolution
+
+    %% Separate the fft integers and not-integers
+    for i = 1:length(f2)
+        if mod(f2(i), 1) == 0
+            X2dec(i) = 0;
+            X2int(i) = X2(i);
+        else
+            X2dec(i) = X2(i);
+            X2int(i) = 0;
+        end
+    end
+
+    %% Case length(f2) is odd -> the mirror image of the FFT is reflected between 2 harmonique
+    if mod(length(f2),2) == 1
+        for i = length(f2)/2+1.5:length(f2)
+            if mod(f2(i-1), 1) == 0
+                X2dec(i) = 0;
+                X2int(i) = X2(i);
+            else
+                X2dec(i) = X2(i);
+                X2int(i) = 0;
+            end
+        end
+    else % Case length(f2) is even -> the mirror image of the FFT is reflected at the Nyquist frequency
+        for i = length(f2)/2+1:length(f2)
+            if mod(f2(i), 1) == 0
+                X2dec(i) = 0;
+                X2int(i) = X2(i);
+            else
+                X2dec(i) = X2(i);
+                X2int(i) = 0;
+            end
+        end
+    end
+
+    %%
+    X2int(1) = 0; %remove the data average/dc component
+    X2int(NbTurn+1) = 0; %Remove fondamental/eccentricity
+    % X2int(length(f2)) = 0; %remove the data average/dc component
+    X2int(length(f2)-NbTurn+1) = 0; %Remove eccentricity
+
+
+    %% Extract the fondamental -> exentricity
+    for i = 1:length(f2)
+        if i == NbTurn+1 || i == length(f2)-NbTurn+1
+            X2fond(i) = X2(i);
+        else
+            X2fond(i) = 0;
+        end
+    end
+
+    X2tot = X2int + X2dec;
+
+    Wxfond = real((res_per_rev*NbTurn)/2 * ifft(X2fond)); % Convert data to "time" domain and scale accordingly
+    Wxint = real((res_per_rev*NbTurn)/2 * ifft(X2int));
+    Wxdec = real((res_per_rev*NbTurn)/2 * ifft(X2dec));
+    Wxtot = real((res_per_rev*NbTurn)/2 * ifft(X2tot));
+
+    %%
+    fig = figure();
+
+    subplot(3, 2, 5);
+    bar(f2(1:50*NbTurn),1000*abs(X2int(1:50*NbTurn)),3);
+    axis([0,50,0,1000*max(abs(X2int(1:50*NbTurn)))]);
+    title ('Fourier integer');
+    xlabel('UPR'); ylabel ('nm')
+
+    subplot(3, 2, 6);
+    bar(f2(1:50*NbTurn),1000*abs(X2dec(1:50*NbTurn)),2);
+    title('Fourier non-integer');
+    axis([0,50,0,1000*max(abs(X2dec(1:50*NbTurn)))]);
+    title('Fourier non-integer');
+    xlabel('UPR'); ylabel ('nm')
+
+    % Eccentricity
+    Eccentricity = max(Wxfond) - min(Wxfond);
+
+    % total error motion X
+    Total_Error = max(Wxtot) - min(Wxtot);
+
+    % lsc X synchronous
+    Synchronous_Error = max(Wxint) - min(Wxint);
+
+    % lsc X Asynchronous
+    var = reshape(Wxdec,length(Wxdec)/NbTurn,NbTurn);
+    for i = 1:length(Wxdec)/NbTurn
+        Asynch(i) = max(var(i,:)) - min(var(i,:));
+    end
+    Asynchronous_Error = max(Asynch) - min(Asynch);
+
+    % Raw Error Motion without Exentricity (sync +asynch)
+    subplot(3, 2, 2);
+    polar2(Theta,Wxtot, 'b');
+    title ('Total error');
+    % Residual Synchronous Error Motion without Exentricity (ie fondamental sync err motion)
+    subplot(3, 2, 3);
+    polar2(Theta,Wxint,'b');
+    title('Residual synchronous error');
+    % Asynchronous Error Motion
+    subplot(3, 2, 4);
+    polar2(Theta,Wxdec, 'b');
+    title('Asynchronous error');
+
+    strmin0 = ['Eccentricity= ', num2str(Eccentricity), ' \mum '];
+    strmin1 = ['Total error = ', num2str(Total_Error*1000), ' nm'];
+    strmin2 = ['Residual synchronous error = ', num2str(Synchronous_Error*1000), ' nm' ];
+    strmin3 = ['Asynchronous error = ', num2str(Asynchronous_Error*1000), ' nm'];
+    dim0 = [0.1 0.55 0.3 .3];%x y w h basgauche to hautdroite
+    dim1 = [0.15 0.65 0.3 .3];
+    annotation('textbox',dim0, 'String',{ strmin0 ,strmin1 , strmin2, strmin3}, 'FitBoxToText', 'on')
+    annotation('textbox',dim1, 'String',texte, 'FitBoxToText', 'on')
+
+    saveas(fig,fullfile(path,char(texte)),'jpg');
+    Res = 1;
+    close all;
+end

spindle/Macros_ttt_spindle/Tilt_Spindle_error.m

@@ -0,0 +1,146 @@
+function Res = Tilt_Spindle_error(dataX, dataX2,NbTurn, texte, path)
+    L = length(dataX);
+    res_per_rev = L/NbTurn;
+    P = 0: (res_per_rev * NbTurn)-1;
+    Pos = P'*360/res_per_rev;
+    Theta = deg2rad(Pos)';
+
+    D = 76.2; %distance entre les deux balls en milimetres
+
+    x1 = myfit2(Pos, dataX);
+    y1 = myfit2(Pos, dataX2);
+
+
+
+    %Convert data to frequency domain and scale accordingly
+    X2 = 2/(res_per_rev*NbTurn)*fft(x1);
+    f2 = (0:L-1)./NbTurn;
+    Y2 = 2/(res_per_rev*NbTurn)*fft(y1);
+
+
+    % Separate the fft integers and not-integers
+    for i = 1:length(f2)
+        if mod(f2(i), 1) == 0
+            X2dec(i) = 0;
+            X2int(i) = X2(i);
+            Y2dec(i) = 0;
+            Y2int(i) = Y2(i);
+        else
+            X2dec(i) = X2(i);
+            X2int(i) = 0;
+            Y2dec(i) = Y2(i);
+            Y2int(i) = 0;
+        end
+    end
+
+    if mod(length(f2),2) == 1 % Case length(f2) is odd -> the mirror image of the FFT is reflected between 2 harmonique
+        for i = length(f2)/2+1.5:length(f2)
+            if mod(f2(i-1), 1) == 0
+                X2dec(i) = 0;
+                X2int(i) = X2(i);
+                Y2dec(i) = 0;
+                Y2int(i) = Y2(i);
+            else
+                X2dec(i) = X2(i);
+                X2int(i) = 0;
+                Y2dec(i) = Y2(i);
+                Y2int(i) = 0;
+            end
+        end
+    else % Case length(f2) is even -> the mirror image of the FFT is reflected at the Nyquist frequency
+        for i = length(f2)/2+1:length(f2)
+            if mod(f2(i), 1) == 0;
+                X2dec(i) = 0;
+                X2int(i) = X2(i);
+                Y2dec(i) = 0;
+                Y2int(i) = Y2(i);
+            else
+                X2dec(i) = X2(i);
+                X2int(i) = 0;
+                Y2dec(i) = Y2(i);
+                Y2int(i) = 0;
+            end
+        end
+    end
+
+    X2int(1) = 0; %remove the data average/dc component
+    X2int(NbTurn+1) = 0; %Remove fondamental/eccentricity
+    %  X2int(length(f2)) = 0; %remove the data average/dc component
+    X2int(length(f2)-NbTurn+1) = 0; %Remove eccentricity
+
+    Y2int(1) = 0; %remove the data average/dc component
+    Y2int(NbTurn+1) = 0; %Remove fondamental/eccentricity
+    %  Y2int(length(f2)) = 0; %remove the data average/dc component
+    Y2int(length(f2)-NbTurn+1) = 0; %Remove eccentricity
+
+
+    % Extract the fondamentale-> exentricity
+    for i = 1:length(f2)
+        if i == NbTurn+1 || i==  length(f2)-NbTurn + 1
+            X2fond(i) = X2(i);
+            Y2fond(i) = Y2(i);
+        else
+            X2fond(i) = 0;
+            Y2fond(i) = 0;
+        end
+    end
+
+    X2tot = X2int + X2dec;
+    Y2tot = Y2int + Y2dec;
+
+    %Convert data to "time" domain and scale accordingly
+    Wxint = real((res_per_rev*NbTurn)/2*ifft(X2int));
+    Wxdec = real((res_per_rev*NbTurn)/2*ifft(X2dec));
+    Wxtot = real((res_per_rev*NbTurn)/2*ifft(X2tot));
+
+    %Convert data to "time" domain and scale accordingly
+    Wyint = real((res_per_rev*NbTurn)/2*ifft(Y2int));
+    Wydec = real((res_per_rev*NbTurn)/2*ifft(Y2dec));
+    Wytot = real((res_per_rev*NbTurn)/2*ifft(Y2tot));
+
+    Tint = atan((Wyint - Wxint)/(D*1000));
+    Tdec = atan((Wydec - Wxdec)/(D*1000));
+    Ttot = atan((Wytot - Wxtot)/(D*1000));
+
+    %%
+    fig = figure();
+
+    % total error motion
+    Total_Error = max(Ttot)- min(Ttot);
+
+    %lsc X synchronous
+    Synchronous_Error = max(Tint)- min(Tint);
+
+    %lsc X Asynchronous
+    var = reshape(Tdec,length(Tdec)/NbTurn,NbTurn);
+    for i = 1:length(Tdec)/NbTurn
+        Asynch(i) = max(var(i,:)) - min(var(i,:)) ;
+    end
+    Asynchronous_Error = max(Asynch)- min(Asynch);
+
+    % Raw Error Motion without Exentricity (sync +asynch)
+    subplot(2, 2, 2);
+    polar2(Theta,Ttot, 'b');
+    title('Total error');
+    % Residual Synchronous Error Motion without Exentricity (ie fondamental sync err motion)
+    subplot(2, 2, 3);
+    polar2(Theta,Tint,'b');
+    title('Residual synchronous error');
+    % Asynchronous Error Motion
+    subplot(2, 2, 4);
+    polar2(Theta,Tdec, 'b');
+    title ('Asynchronous error');
+
+    %%
+    strmin1 = ['Total error = ', num2str(Total_Error*1000000), ' \murad'];
+    strmin2 = ['Residual synchronous error = ', num2str(Synchronous_Error*1000000), ' \murad' ];
+    strmin3 = ['Asynchronous error = ', num2str(Asynchronous_Error*1000000), ' \murad'];
+    dim0 =[0.04 0.5 0.3 .3];%x y w h basgauche to hautdroite
+    dim1 =[0.15 0.65 0.3 .3];
+    annotation('textbox',dim0, 'String',{ strmin1 , strmin2, strmin3}, 'FitBoxToText', 'on')
+    annotation('textbox',dim1, 'String',texte, 'FitBoxToText', 'on')
+
+    saveas(fig,fullfile(path,char(texte)),'jpg');
+    Res = 1;
+    close all;
+end

spindle/Macros_ttt_spindle/exportFigure.m

@@ -0,0 +1,49 @@
+% This function exports the given figure with nice settings to the given
+% file path in the defined file format.
+% 
+% input parameter:
+%   - handle ... figure handle
+%   - path   ... file path
+%   - format ... file format ('pdf', 'meta', 'jpeg', 'png'), default 'pdf'
+
+function exportFigure(handle, path, format)
+
+% select figure
+figure(handle)
+
+% set output settings
+set(gcf, 'PaperPositionMode', 'manual', 'Renderer', 'painters', 'RendererMode', 'manual');
+
+
+% export figure
+print(handle, path, '-dpdf');
+
+% set export format
+% if(strcmp(format,'pdf') || strcmp(format,'jpeg') )
+%     dformat = ['-d' format];
+% else
+%     dformat = '-dpdf';
+% end
+% 
+% if(strcmp(format,'meta'))
+%     format = 'emf';
+% end
+% 
+% if(strcmp(format,'jpeg'))
+%     format = 'jpg';
+% end
+
+if(strcmp(format,'jpeg') || strcmp(format,'jpg') || ...
+        strcmp(format,'JPEG') || strcmp(format,'JPG'))
+    display(' --------------- Converting to JPG image ----------------- ');
+    system(['convert -density 300 ' path '.pdf -quality 95 ' path '.jpg']);
+    system(['rm ' path '.pdf']);
+end
+
+if(strcmp(format,'png') || strcmp(format,'PNG'))
+    display(' --------------- Converting to PNG image ----------------- ');
+    system(['convert -density 300 ' path '.pdf ' path '.png']);
+    system(['rm ' path '.pdf']);
+end
+
+end % end function

spindle/Macros_ttt_spindle/fitcircle2.m

@@ -0,0 +1,15 @@
+function [c,r] = fitcircle2(x,y)
+% Function to plot a linear least squares circle on a data set
+% x = vector of x-coordinates
+% y = vector of y-coordinates
+% c = center of the circle
+% r = radius of circle (not relevant)
+n=length(x);  xx=x.*x; yy=y.*y; xy=x.*y;
+A=[sum(x) sum(y) n;sum(xy) sum(yy) sum(y);sum(xx) sum(xy) sum(x)];
+B=[-sum(xx+yy) ; -sum(xx.*y+yy.*y) ; -sum(xx.*x+xy.*y)];
+a=A\B;
+xc = -.5*a(1);
+yc = -.5*a(2);
+r  =  sqrt((a(1)^2+a(2)^2)/4-a(3));
+c = [xc;yc];
+end

spindle/Macros_ttt_spindle/importfile.m

@@ -0,0 +1,64 @@
+function [mesure1,VarName4,VarName6,VarName8,VarName10] = importfile1(filename, startRow, endRow)
+%IMPORTFILE1 Import numeric data from a text file as column vectors.
+%   [MESURE1,VARNAME4,VARNAME6,VARNAME8,VARNAME10] = IMPORTFILE1(FILENAME)
+%   Reads data from text file FILENAME for the default selection.
+%
+%   [MESURE1,VARNAME4,VARNAME6,VARNAME8,VARNAME10] = IMPORTFILE1(FILENAME,
+%   STARTROW, ENDROW) Reads data from rows STARTROW through ENDROW of text
+%   file FILENAME.
+%
+% Example:
+%   [mesure1,VarName4,VarName6,VarName8,VarName10] = importfile1('mesure.txt',12, 550411);
+%
+%    See also TEXTSCAN.
+
+% Auto-generated by MATLAB on 2017/05/03 10:14:16
+
+%% Initialize variables.
+delimiter = '\t';
+if nargin<=2
+    startRow = 12;
+    endRow = inf;
+end
+
+%% Format string for each line of text:
+%   column2: double (%f)
+%	column4: double (%f)
+%   column6: double (%f)
+%	column8: double (%f)
+%   column10: double (%f)
+% For more information, see the TEXTSCAN documentation.
+formatSpec = '%*s%f%*s%f%*s%f%*s%f%*s%f%[^\n\r]';
+
+%% Open the text file.
+fileID = fopen(filename,'r');
+
+%% Read columns of data according to format string.
+% This call is based on the structure of the file used to generate this
+% code. If an error occurs for a different file, try regenerating the code
+% from the Import Tool.
+dataArray = textscan(fileID, formatSpec, endRow(1)-startRow(1)+1, 'Delimiter', delimiter, 'EmptyValue' ,NaN,'HeaderLines', startRow(1)-1, 'ReturnOnError', false);
+for block=2:length(startRow)
+    frewind(fileID);
+    dataArrayBlock = textscan(fileID, formatSpec, endRow(block)-startRow(block)+1, 'Delimiter', delimiter, 'EmptyValue' ,NaN,'HeaderLines', startRow(block)-1, 'ReturnOnError', false);
+    for col=1:length(dataArray)
+        dataArray{col} = [dataArray{col};dataArrayBlock{col}];
+    end
+end
+
+%% Close the text file.
+fclose(fileID);
+
+%% Post processing for unimportable data.
+% No unimportable data rules were applied during the import, so no post
+% processing code is included. To generate code which works for
+% unimportable data, select unimportable cells in a file and regenerate the
+% script.
+
+%% Allocate imported array to column variable names
+mesure1 = dataArray{:, 1};
+VarName4 = dataArray{:, 2};
+VarName6 = dataArray{:, 3};
+VarName8 = dataArray{:, 4};
+VarName10 = dataArray{:, 5};
+

spindle/Macros_ttt_spindle/myFit.m

@@ -0,0 +1,35 @@
+function [Y] = myFit(x, y, start, End, file)
+
+W = (y-y(1)) - (y(End-start)-y(1))/(x(End - start))*x;
+A = [x ones(size(x))]\y;
+
+a = A(1); b = A(2);
+Y = y - (a*x + b);
+
+fig = newFigure(16,20);
+subplot(2,1,1);
+plot (x, y, 'r');
+grid on
+xlabel ('Time (sec)', 'interpreter', 'latex')
+ylabel ('Straightness $( \mu{} \mathrm{m})$', 'interpreter', 'latex')
+title ('Straightness-Raw', 'interpreter', 'latex')
+
+subplot(2,1,2);
+plot (x, Y, 'color',[0 0.5 0]);
+grid on
+%ylim([-10 14]);
+xlabel ('Time (sec)', 'interpreter', 'latex')
+ylabel ('Straightness $( \mu{} \mathrm{m})$', 'interpreter', 'latex')
+title ('Straightness-Corrected', 'interpreter', 'latex')
+
+m= max(Y) - min(Y);
+text(1,10, 'max - min =', 'interpreter', 'latex')
+text(1,6, [num2str(m) ' $\mu{} \mathrm{m}$'], 'interpreter', 'latex')
+
+string = { file };
+Newstring = {string{1}(1:end-4)} ;
+p= char(Newstring);
+
+exportFigure(fig,p, 'png');
+end 
+

spindle/Macros_ttt_spindle/myfit2.m

@@ -0,0 +1,5 @@
+function Y =  myfit2(x,y)
+    A = [x ones(size(x))]\y;
+    a = A(1); b = A(2);
+    Y = y - (a*x + b);
+end

spindle/Macros_ttt_spindle/polar2.m

@@ -0,0 +1,270 @@
+function hpol = polar2(varargin)
+%POLAR2  Polar coordinate plot.
+%   POLA2R(THETA, RHO) makes a plot using polar coordinates of
+%   the angle THETA, in radians, versus the radius RHO.
+%   POLAR2(THETA,RHO,R) uses the radial limits specified by the two element
+%   vector R. If R is a 4-element vector, the second two elements are
+%   angular limits.
+%   POLAR2(THETA,RHO,S) uses the linestyle specified in string S.
+%   See PLOT for a description of legal linestyles.
+%   POLAR2(THETA,RHO,R,S) uses the linestyle specified in string S and the
+%   radial limits in R.
+%
+%   POLAR2(AX,...) plots into AX instead of GCA.
+%
+%   H = POLAR2(...) returns a handle to the plotted object in H.
+%
+%   Example:
+%      t = 0:.01:2*pi;
+%      polar2(t,sin(2*t).*cos(2*t),[-2 2], '--r');
+%
+%   See also PLOT, LOGLOG, SEMILOGX, SEMILOGY.
+%
+%   Copyright 2009-2012 The MathWorks, Inc.
+
+
+% Parse possible Axes input
+[cax,args,nargs] = axescheck(varargin{:});
+%error(nargchk(1,4,nargs,'struct'));
+
+if nargs < 1 || nargs > 4
+    error('MATLAB:polar:InvalidInput', 'Requires 2 to 4 data arguments.')
+elseif nargs == 2
+    theta = args{1};
+    rho = args{2};
+    if ischar(rho)
+        line_style = rho;
+        rho = theta;
+        [mr,nr] = size(rho);
+        if mr == 1
+            theta = 1:nr;
+        else
+            th = (1:mr)';
+            theta = th(:,ones(1,nr));
+        end
+    else
+        line_style = 'auto';
+    end
+    radial_limits = [];
+elseif nargs == 1
+    theta = args{1};
+    line_style = 'auto';
+    rho = theta;
+    [mr,nr] = size(rho);
+    if mr == 1
+        theta = 1:nr;
+    else
+        th = (1:mr)';
+        theta = th(:,ones(1,nr));
+    end
+    radial_limits = [];
+elseif nargs == 3
+    if ( ischar(args{3}) )
+        [theta,rho,line_style] = deal(args{1:3});
+        radial_limits = [];
+    else
+        [theta,rho,radial_limits] = deal(args{1:3});
+        line_style = 'auto';
+        if ( isempty(radial_limits) )
+      radial_limits = [];
+    elseif ( numel(radial_limits) == 2 )
+      radial_limits = [radial_limits(:); 0; 2*pi];
+    elseif ( ~(numel(radial_limits) == 4) )
+            error ( 'R must be a 2 element vector' );
+        end
+    end
+else %nargs == 4
+  [theta,rho,radial_limits,line_style] = deal(args{1:4});
+  if ( isempty(radial_limits) )
+    radial_limits = [];
+  elseif ( numel(radial_limits) == 2 )
+    radial_limits = [radial_limits(:); 0; 2*pi];
+  elseif ( ~(numel(radial_limits) == 4) )
+    error ( 'R must be a 2 element vector' );
+  end
+end
+
+if ischar(theta) || ischar(rho)
+    error('MATLAB:polar:InvalidInputType', 'Input arguments must be numeric.');
+end
+if ~isequal(size(theta),size(rho))
+    error('MATLAB:polar:InvalidInput', 'THETA and RHO must be the same size.');
+end
+
+% get hold state
+cax = newplot(cax);
+
+next = lower(get(cax,'NextPlot'));
+hold_state = ishold(cax);
+
+% get x-axis text color so grid is in same color
+tc = get(cax,'xcolor');
+ls = get(cax,'gridlinestyle');
+
+% Hold on to current Text defaults, reset them to the
+% Axes' font attributes so tick marks use them.
+fAngle  = get(cax, 'DefaultTextFontAngle');
+fName   = get(cax, 'DefaultTextFontName');
+fSize   = get(cax, 'DefaultTextFontSize');
+fWeight = get(cax, 'DefaultTextFontWeight');
+fUnits  = get(cax, 'DefaultTextUnits');
+set(cax, 'DefaultTextFontAngle',  get(cax, 'FontAngle'), ...
+    'DefaultTextFontName',   get(cax, 'FontName'), ...
+    'DefaultTextFontSize',   get(cax, 'FontSize'), ...
+    'DefaultTextFontWeight', get(cax, 'FontWeight'), ...
+    'DefaultTextUnits','data')
+
+% only do grids if hold is off
+if ~hold_state
+
+% make a radial grid
+    hold(cax,'on');
+    set(cax,'dataaspectratio',[1 1 1],'plotboxaspectratiomode','auto')
+
+    % ensure that Inf values don't enter into the limit calculation.
+    %arho = abs(rho(:));
+    if ( isempty(radial_limits) )
+        maxrho = max(rho(rho ~= Inf));
+        minrho = min(rho(rho ~= Inf));
+        hhh=line([minrho minrho maxrho maxrho],[minrho maxrho maxrho minrho],'parent',cax);
+        v = [get(cax,'xlim') get(cax,'ylim')];
+        ticks = numel(get(cax,'ytick'));
+        delete(hhh);
+        % check radial limits and ticks
+        rmin = v(1); rmax = v(4); rticks = max(ticks-1,2);
+        if rticks > 5   % see if we can reduce the number
+            if rem(rticks,2) == 0
+                rticks = rticks/2;
+            elseif rem(rticks,3) == 0
+                rticks = rticks/3;
+            end
+        end
+        rinc = (rmax-rmin)/rticks;
+
+  else
+
+    thmax = radial_limits(4);
+        thmin = radial_limits(3);
+    rmax = radial_limits(2);
+        rmin = radial_limits(1);
+
+        order = (10^floor(log10(rmax-rmin)));
+        firstDigit = floor((rmax-rmin)/order);
+        if ( firstDigit <= 1 )
+            step = 0.2*order;
+        elseif ( firstDigit <= 3 )
+            step = 0.5*order;
+        elseif ( firstDigit <= 7 )
+            step = order;
+        else
+            step = 2*order;
+        end
+        rinc = step;
+
+    %Also, make sure to drop any values outside the limits.
+    if ( rmin == 0 )
+      subset = rho >= -rmax & rho <= rmax & theta >= thmin & theta <= thmax;
+    else
+      subset = rho >= rmin & rho <= rmax & theta >= thmin & theta <= thmax;
+    end
+    theta = theta(subset);
+    rho = rho(subset);
+    end
+
+% define a circle
+    th = 0:pi/50:2*pi;
+    xunit = cos(th);
+    yunit = sin(th);
+% now really force points on x/y axes to lie on them exactly
+    inds = 1:(length(th)-1)/4:length(th);
+    xunit(inds(2:2:4)) = zeros(2,1);
+    yunit(inds(1:2:5)) = zeros(3,1);
+% plot background if necessary
+    if ~ischar(get(cax,'color')),
+       patch('xdata',xunit*(rmax-rmin),'ydata',yunit*(rmax-rmin), ...
+             'edgecolor',tc,'facecolor',get(cax,'color'),...
+             'handlevisibility','off','parent',cax);
+    end
+
+% draw radial circles
+    c82 = cos(82*pi/180);
+    s82 = sin(82*pi/180);
+    for i=(rmin+rinc):rinc:rmax
+        hhh = line(xunit*(i-rmin),yunit*(i-rmin),'linestyle',ls,'color',tc,'linewidth',1,...
+                   'handlevisibility','off','parent',cax);
+        text((i-rmin+rinc/20)*c82,(i-rmin+rinc/20)*s82, ...
+            ['  ' num2str(i)],'verticalalignment','bottom',...
+            'handlevisibility','off','parent',cax)
+    end
+    set(hhh,'linestyle',':') % Make outer circle solid
+
+% plot spokes
+    th = (1:6)*2*pi/12;
+    cst = cos(th); snt = sin(th);
+    cs = [-cst; cst];
+    sn = [-snt; snt];
+    line((rmax-rmin)*cs,(rmax-rmin)*sn,'linestyle',ls,'color',tc,'linewidth',1,...
+         'handlevisibility','off','parent',cax)
+
+% annotate spokes in degrees
+    rt = 1.1*(rmax-rmin);
+    for i = 1:length(th)
+        text(rt*cst(i),rt*snt(i),int2str(i*30),...
+             'horizontalalignment','center',...
+             'handlevisibility','off','parent',cax);
+        if i == length(th)
+            loc = int2str(0);
+        else
+            loc = int2str(180+i*30);
+        end
+        text(-rt*cst(i),-rt*snt(i),loc,'horizontalalignment','center',...
+             'handlevisibility','off','parent',cax)
+    end
+
+% set view to 2-D
+    view(cax,2);
+% set axis limits
+    axis(cax,(rmax-rmin)*[-1 1 -1.15 1.15]);
+
+    setappdata( cax, 'rMin', rmin );
+
+else
+    %Try to find the inner radius of the current axis.
+    if ( isappdata ( cax, 'rMin' ) )
+        rmin = getappdata( cax, 'rMin' );
+    else
+        rmin = 0;
+    end
+end
+
+% Reset defaults.
+set(cax, 'DefaultTextFontAngle', fAngle , ...
+    'DefaultTextFontName',   fName , ...
+    'DefaultTextFontSize',   fSize, ...
+    'DefaultTextFontWeight', fWeight, ...
+    'DefaultTextUnits',fUnits );
+
+% transform data to Cartesian coordinates.
+xx = (rho - rmin).*cos(theta);
+yy = (rho - rmin).*sin(theta);
+
+% plot data on top of grid
+if strcmp(line_style,'auto')
+    q = plot(xx,yy,'parent',cax);
+else
+    q = plot(xx,yy,line_style,'parent',cax);
+end
+
+if nargout == 1
+    hpol = q;
+end
+
+if ~hold_state
+    set(cax,'dataaspectratio',[1 1 1]), axis(cax,'off'); set(cax,'NextPlot',next);
+end
+set(get(cax,'xlabel'),'visible','on')
+set(get(cax,'ylabel'),'visible','on')
+
+if ~isempty(q) && ~isdeployed
+    makemcode('RegisterHandle',cax,'IgnoreHandle',q,'FunctionName','polar');
+end

spindle/figs/nohup.out

@@ -0,0 +1,3 @@
+invalid color string: rgba(50, 48, 47)
+
+(termite:2378): GLib-WARNING **: 15:04:00.329: GChildWatchSource: Exit status of a child process was requested but ECHILD was received by waitpid(). See the documentation of g_child_watch_source_new() for possible causes.

spindle/index.org

@@ -0,0 +1,757 @@
+#+TITLE: Spindle Analysis
+#+SETUPFILE: ../config.org
+
+The report made by the PEL is accessible [[file:documents/Spindle_report_test.pdf][here]].
+
+* Notes
+#+name: fig:setup_spindle
+#+caption: Measurement setup at the PEL lab
+[[file:./img/setup_spindle.png]]
+
+* Data Processing
+  :PROPERTIES:
+  :header-args:matlab+: :tangle matlab/spindle_data_processing.m
+  :header-args:matlab+: :comments org :mkdirp yes
+  :END:
+  <<sec:spindle_data_processing>>
+
+#+begin_src bash :exports none :results none
+  if [ matlab/spindle_data_processing.m -nt data/spindle_data_processing.zip ]; then
+    cp matlab/spindle_data_processing.m spindle_data_processing.m;
+    zip data/spindle_data_processing \
+        mat/10turns_1rpm_icepap.txt \
+        mat/10turns_60rpm_IcepapFIR.txt \
+        src/getAsynchronousError.m \
+        spindle_data_processing.m
+    rm spindle_data_processing.m;
+  fi
+#+end_src
+
+#+begin_note
+  All the files (data and Matlab scripts) are accessible [[file:data/spindle_data_processing.zip][here]].
+#+end_note
+
+** Matlab Init                                              :noexport:ignore:
+#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
+  <<matlab-dir>>
+#+end_src
+
+#+begin_src matlab :exports none :results silent :noweb yes
+  <<matlab-init>>
+#+end_src
+
+#+begin_src matlab :exports none :results silent :noweb yes
+  %% Add the getAsynchronousError to path
+  addpath('./src/');
+#+end_src
+
+** Load Measurement Data
+#+begin_src matlab :results none :exports code
+  spindle_1rpm_table  = readtable('./mat/10turns_1rpm_icepap.txt');
+  spindle_60rpm_table = readtable('./mat/10turns_60rpm_IcepapFIR.txt');
+#+end_src
+
+#+begin_src matlab :results output :exports code
+  spindle_1rpm_table(1, :)
+#+end_src
+
+#+begin_src matlab :results none :exports code
+  spindle_1rpm  = table2array(spindle_1rpm_table);
+  spindle_60rpm = table2array(spindle_60rpm_table);
+#+end_src
+
+** Convert Signals from [deg] to [sec]
+#+begin_src matlab :results none :exports code
+  speed_1rpm = 360/60; % [deg/sec]
+  spindle_1rpm(:, 1) = spindle_1rpm(:, 2)/speed_1rpm;  % From position [deg] to time [s]
+
+  speed_60rpm = 360/1; % [deg/sec]
+  spindle_60rpm(:, 1) = spindle_60rpm(:, 2)/speed_60rpm; % From position [deg] to time [s]
+#+end_src
+
+** Convert Signals
+#+begin_src matlab :results none :exports code
+  % scaling = 1/80000; % 80 mV/um
+  scaling = 1e-6; % [um] to [m]
+
+  spindle_1rpm(:, 3:end) = scaling*spindle_1rpm(:, 3:end); % [V] to [m]
+  spindle_1rpm(:, 3:end) = spindle_1rpm(:, 3:end)-mean(spindle_1rpm(:, 3:end)); % Remove mean
+
+  spindle_60rpm(:, 3:end) = scaling*spindle_60rpm(:, 3:end); % [V] to [m]
+  spindle_60rpm(:, 3:end) = spindle_60rpm(:, 3:end)-mean(spindle_60rpm(:, 3:end)); % Remove mean
+#+end_src
+
+** Ts and Fs for both measurements
+#+begin_src matlab :results none :exports code
+  Ts_1rpm = spindle_1rpm(end, 1)/(length(spindle_1rpm(:, 1))-1);
+  Fs_1rpm = 1/Ts_1rpm;
+
+  Ts_60rpm = spindle_60rpm(end, 1)/(length(spindle_60rpm(:, 1))-1);
+  Fs_60rpm = 1/Ts_60rpm;
+#+end_src
+
+** Find Noise of the ADC [$\frac{m}{\sqrt{Hz}}$]
+#+begin_src matlab :results none :exports code
+  data = spindle_1rpm(:, 5);
+  dV_1rpm = min(abs(data(1) - data(data ~= data(1))));
+  noise_1rpm = dV_1rpm/sqrt(12*Fs_1rpm/2);
+
+  data = spindle_60rpm(:, 5);
+  dV_60rpm = min(abs(data(50) - data(data ~= data(50))));
+  noise_60rpm = dV_60rpm/sqrt(12*Fs_60rpm/2);
+#+end_src
+
+** Save all the data under spindle struct
+#+begin_src matlab :results none :exports code
+  spindle.rpm1.time = spindle_1rpm(:, 1);
+  spindle.rpm1.deg  = spindle_1rpm(:, 2);
+  spindle.rpm1.Ts   = Ts_1rpm;
+  spindle.rpm1.Fs   = 1/Ts_1rpm;
+  spindle.rpm1.x    = spindle_1rpm(:, 3);
+  spindle.rpm1.y    = spindle_1rpm(:, 4);
+  spindle.rpm1.z    = spindle_1rpm(:, 5);
+  spindle.rpm1.adcn = noise_1rpm;
+
+  spindle.rpm60.time = spindle_60rpm(:, 1);
+  spindle.rpm60.deg  = spindle_60rpm(:, 2);
+  spindle.rpm60.Ts   = Ts_60rpm;
+  spindle.rpm60.Fs   = 1/Ts_60rpm;
+  spindle.rpm60.x    = spindle_60rpm(:, 3);
+  spindle.rpm60.y    = spindle_60rpm(:, 4);
+  spindle.rpm60.z    = spindle_60rpm(:, 5);
+  spindle.rpm60.adcn = noise_60rpm;
+#+end_src
+
+** Compute Asynchronous data
+#+begin_src matlab :results none :exports code
+  for direction = {'x', 'y', 'z'}
+      spindle.rpm1.([direction{1}, 'async']) = getAsynchronousError(spindle.rpm1.(direction{1}), 10);
+      spindle.rpm60.([direction{1}, 'async']) = getAsynchronousError(spindle.rpm60.(direction{1}), 10);
+  end
+#+end_src
+
+** Save data
+#+begin_src matlab
+  save('./mat/spindle_data.mat', 'spindle');
+#+end_src
+
+* Time Domain Data
+  :PROPERTIES:
+  :header-args:matlab+: :tangle matlab/spindle_time_domain.m
+  :header-args:matlab+: :comments org :mkdirp yes
+  :END:
+  <<sec:spindle_time_domain>>
+
+#+begin_src bash :exports none :results none
+  if [ matlab/spindle_time_domain.m -nt data/spindle_time_domain.zip ]; then
+    cp matlab/spindle_time_domain.m spindle_time_domain.m;
+    zip data/spindle_time_domain \
+        mat/spindle_data.mat \
+        spindle_time_domain.m
+    rm spindle_time_domain.m;
+  fi
+#+end_src
+
+#+begin_note
+  All the files (data and Matlab scripts) are accessible [[file:data/spindle_time_domain.zip][here]].
+#+end_note
+
+** Matlab Init                                              :noexport:ignore:
+#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
+<<matlab-dir>>
+#+end_src
+
+#+begin_src matlab :exports none :results none :noweb yes
+<<matlab-init>>
+#+end_src
+
+** Load the processed data
+#+begin_src matlab
+  load('./mat/spindle_data.mat', 'spindle');
+#+end_src
+
+** Plot X-Y-Z position with respect to Time - 1rpm
+#+begin_src matlab :results none :exports code
+  figure;
+  hold on;
+  plot(spindle.rpm1.time, spindle.rpm1.x);
+  plot(spindle.rpm1.time, spindle.rpm1.y);
+  plot(spindle.rpm1.time, spindle.rpm1.z);
+  hold off;
+  xlabel('Time [s]'); ylabel('Amplitude [m]');
+  legend({'tx - 1rpm', 'ty - 1rpm', 'tz - 1rpm'});
+#+end_src
+
+#+NAME: fig:spindle_xyz_1rpm
+#+HEADER: :tangle no :exports results :results raw :noweb yes
+#+begin_src matlab :var filepath="figs/spindle_xyz_1rpm.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:spindle_xyz_1rpm
+#+CAPTION: Raw time domain translation - 1rpm
+#+RESULTS: fig:spindle_xyz_1rpm
+[[file:figs/spindle_xyz_1rpm.png]]
+
+** Plot X-Y-Z position with respect to Time - 60rpm
+The measurements for the spindle turning at 60rpm are shown figure [[fig:spindle_xyz_60rpm]].
+
+#+begin_src matlab :results none :exports code
+  figure;
+  hold on;
+  plot(spindle.rpm60.time, spindle.rpm60.x);
+  plot(spindle.rpm60.time, spindle.rpm60.y);
+  plot(spindle.rpm60.time, spindle.rpm60.z);
+  hold off;
+  xlabel('Time [s]'); ylabel('Amplitude [m]');
+  legend({'tx - 60rpm', 'ty - 60rpm', 'tz - 60rpm'});
+#+end_src
+
+#+NAME: fig:spindle_xyz_60rpm
+#+HEADER: :tangle no :exports results :results raw :noweb yes
+#+begin_src matlab :var filepath="figs/spindle_xyz_60rpm.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:spindle_xyz_60rpm
+#+CAPTION: Raw time domain translation - 60rpm
+#+RESULTS: fig:spindle_xyz_60rpm
+[[file:figs/spindle_xyz_60rpm.png]]
+
+** Plot Synchronous and Asynchronous - 1rpm
+#+begin_src matlab :results none :exports code
+  figure;
+  hold on;
+  plot(spindle.rpm1.time, spindle.rpm1.x);
+  plot(spindle.rpm1.time, spindle.rpm1.xasync);
+  hold off;
+  xlabel('Time [s]'); ylabel('Amplitude [m]');
+  legend({'tx - 1rpm - Sync', 'tx - 1rpm - Async'});
+#+end_src
+
+#+NAME: fig:spindle_1rpm_sync_async
+#+HEADER: :tangle no :exports results :results raw :noweb yes
+#+begin_src matlab :var filepath="figs/spindle_1rpm_sync_async.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:spindle_1rpm_sync_async
+#+CAPTION: Comparison of the synchronous and asynchronous displacements - 1rpm
+#+RESULTS: fig:spindle_1rpm_sync_async
+[[file:figs/spindle_1rpm_sync_async.png]]
+
+** Plot Synchronous and Asynchronous - 60rpm
+The data is split into its Synchronous and Asynchronous part (figure [[fig:spindle_60rpm_sync_async]]). We then use the Asynchronous part for the analysis in the following sections as we suppose that we can deal with the synchronous part with feedforward control.
+
+#+begin_src matlab :results none :exports code
+  figure;
+  hold on;
+  plot(spindle.rpm60.time, spindle.rpm60.x);
+  plot(spindle.rpm60.time, spindle.rpm60.xasync);
+  hold off;
+  xlabel('Time [s]'); ylabel('Amplitude [m]');
+  legend({'tx - 60rpm - Sync', 'tx - 60rpm - Async'});
+#+end_src
+
+#+NAME: fig:spindle_60rpm_sync_async
+#+HEADER: :tangle no :exports results :results raw :noweb yes
+#+begin_src matlab :var filepath="figs/spindle_60rpm_sync_async.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:spindle_60rpm_sync_async
+#+CAPTION: Comparison of the synchronous and asynchronous displacements - 60rpm
+#+RESULTS: fig:spindle_60rpm_sync_async
+[[file:figs/spindle_60rpm_sync_async.png]]
+
+** Plot X against Y
+#+begin_src matlab :results none :exports code
+  figure;
+  hold on;
+  plot(spindle.rpm1.x, spindle.rpm1.y);
+  plot(spindle.rpm60.x, spindle.rpm60.y);
+  hold off;
+  xlabel('X Amplitude [m]'); ylabel('Y Amplitude [m]');
+  legend({'1rpm', '60rpm'});
+#+end_src
+
+#+NAME: fig:spindle_xy_1_60rpm
+#+HEADER: :tangle no :exports results :results raw :noweb yes
+#+begin_src matlab :var filepath="figs/spindle_xy_1_60rpm.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:spindle_xy_1_60rpm
+#+CAPTION: Synchronous x-y displacement
+#+RESULTS: fig:spindle_xy_1_60rpm
+[[file:figs/spindle_xy_1_60rpm.png]]
+
+** Plot X against Y - Asynchronous
+#+begin_src matlab :results none :exports code
+  figure;
+  hold on;
+  plot(spindle.rpm1.xasync, spindle.rpm1.yasync);
+  plot(spindle.rpm60.xasync, spindle.rpm60.yasync);
+  hold off;
+  xlabel('X Amplitude [m]'); ylabel('Y Amplitude [m]');
+  legend({'1rpm', '60rpm'});
+#+end_src
+
+#+NAME: fig:spindle_xy_1_60_rpm_async
+#+HEADER: :tangle no :exports results :results raw :noweb yes
+#+begin_src matlab :var filepath="figs/spindle_xy_1_60_rpm_async.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:spindle_xy_1_60_rpm_async
+#+CAPTION: Asynchronous x-y displacement
+#+RESULTS: fig:spindle_xy_1_60_rpm_async
+[[file:figs/spindle_xy_1_60_rpm_async.png]]
+
+* Model of the spindle
+  :PROPERTIES:
+  :header-args:matlab+: :tangle matlab/spindle_model.m
+  :header-args:matlab+: :comments org :mkdirp yes
+  :END:
+  <<sec:spindle_model>>
+
+#+begin_src bash :exports none :results none
+  if [ matlab/spindle_model.m -nt data/spindle_model.zip ]; then
+    cp matlab/spindle_model.m spindle_model.m;
+    zip data/spindle_model \
+        spindle_model.m
+    rm spindle_model.m;
+  fi
+#+end_src
+
+#+begin_note
+  All the files (data and Matlab scripts) are accessible [[file:data/spindle_model.zip][here]].
+#+end_note
+
+** Schematic of the model
+The model of the spindle used is shown figure [[fig:model_spindle]].
+
+$f$ is the perturbation force of the spindle and $d$ is the measured displacement.
+
+#+name: fig:model_spindle
+#+caption: Model of the Spindle
+#+attr_latex: :float t
+[[./figs/uniaxial-model-spindle.png]]
+
+** Matlab Init                                              :noexport:ignore:
+#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
+  <<matlab-dir>>
+#+end_src
+
+#+begin_src matlab :exports none :results silent :noweb yes
+  <<matlab-init>>
+#+end_src
+
+#+begin_src matlab :exports none :results silent :noweb yes
+  %% Add the getAsynchronousError to path
+  addpath('./src/');
+#+end_src
+
+** Parameters
+#+begin_src matlab :exports code :results none
+  mg = 3000; % Mass of granite [kg]
+  ms = 50;   % Mass of Spindle [kg]
+
+  kg = 1e8; % Stiffness of granite [N/m]
+  ks = 5e7; % Stiffness of spindle [N/m]
+#+end_src
+
+** Compute Mass and Stiffness Matrices
+#+begin_src matlab :exports code :results none
+  Mm = diag([ms, mg]);
+  Km = diag([ks, ks+kg]) - diag(ks, -1) - diag(ks, 1);
+#+end_src
+
+** Compute resonance frequencies
+#+begin_src matlab :exports code :results none
+  A = [zeros(size(Mm)) eye(size(Mm)) ; -Mm\Km zeros(size(Mm))];
+  eigA = imag(eigs(A))/2/pi;
+  eigA = eigA(eigA>0);
+  eigA = eigA(1:2);
+#+end_src
+
+** From model_damping compute the Damping Matrix
+#+begin_src matlab :exports code :results none
+  modal_damping = 1e-5;
+
+  ab = [0.5*eigA(1) 0.5/eigA(1) ; 0.5*eigA(2) 0.5/eigA(2)]\[modal_damping ; modal_damping];
+
+  Cm = ab(1)*Mm +ab(2)*Km;
+#+end_src
+
+** Define inputs, outputs and state names
+#+begin_src matlab :exports code :results none
+  StateName = {...
+      'xs', ... % Displacement of Spindle [m]
+      'xg', ... % Displacement of Granite [m]
+      'vs', ... % Velocity of Spindle [m]
+      'vg', ... % Velocity of Granite [m]
+              };
+  StateUnit = {'m', 'm', 'm/s', 'm/s'};
+
+  InputName = {...
+      'f' ... % Spindle Force [N]
+              };
+  InputUnit = {'N'};
+
+  OutputName = {...
+      'd' ... % Displacement [m]
+               };
+  OutputUnit = {'m'};
+#+end_src
+
+** Define A, B and C matrices
+#+begin_src matlab :exports code :results none
+  % A Matrix
+  A = [zeros(size(Mm)) eye(size(Mm)) ; ...
+       -Mm\Km          -Mm\Cm];
+
+  % B Matrix
+  B_low = zeros(length(StateName), length(InputName));
+  B_low(strcmp(StateName,'vs'), strcmp(InputName,'f')) =  1;
+  B_low(strcmp(StateName,'vg'), strcmp(InputName,'f')) = -1;
+  B = blkdiag(zeros(length(StateName)/2), pinv(Mm))*B_low;
+
+  % C Matrix
+  C = zeros(length(OutputName), length(StateName));
+  C(strcmp(OutputName,'d'), strcmp(StateName,'xs')) =  1;
+  C(strcmp(OutputName,'d'), strcmp(StateName,'xg')) = -1;
+
+  % D Matrix
+  D = zeros(length(OutputName), length(InputName));
+#+end_src
+
+** Generate the State Space Model
+#+begin_src matlab :exports code :results none
+  sys = ss(A, B, C, D);
+  sys.StateName = StateName;
+  sys.StateUnit = StateUnit;
+
+  sys.InputName = InputName;
+  sys.InputUnit = InputUnit;
+
+  sys.OutputName = OutputName;
+  sys.OutputUnit = OutputUnit;
+#+end_src
+
+** Bode Plot
+The transfer function from a disturbance force $f$ to the measured displacement $d$ is shown figure [[fig:spindle_f_to_d]].
+
+#+begin_src matlab :exports code :results none
+  freqs = logspace(-1, 3, 1000);
+
+  figure;
+  plot(freqs, abs(squeeze(freqresp(sys('d', 'f'), freqs, 'Hz'))));
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Amplitude [m/N]');
+#+end_src
+
+#+NAME: fig:spindle_f_to_d
+#+HEADER: :tangle no :exports results :results raw :noweb yes
+#+begin_src matlab :var filepath="figs/spindle_f_to_d.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:spindle_f_to_d
+#+CAPTION: Bode plot of the transfer function from $f$ to $d$
+#+RESULTS: fig:spindle_f_to_d
+[[file:figs/spindle_f_to_d.png]]
+
+** Save the model
+#+begin_src matlab
+  save('./mat/spindle_model.mat', 'sys');
+#+end_src
+
+* Frequency Domain Data
+  :PROPERTIES:
+  :header-args:matlab+: :tangle matlab/spindle_psd.m
+  :header-args:matlab+: :comments org :mkdirp yes
+  :END:
+  <<sec:spindle_psd>>
+
+#+begin_src bash :exports none :results none
+  if [ matlab/spindle_psd.m -nt data/spindle_psd.zip ]; then
+    cp matlab/spindle_psd.m spindle_psd.m;
+    zip data/spindle_psd \
+        mat/spindle_data.mat \
+        mat/spindle_model.mat \
+        spindle_psd.m
+    rm spindle_psd.m;
+  fi
+#+end_src
+
+#+begin_note
+  All the files (data and Matlab scripts) are accessible [[file:data/spindle_psd.zip][here]].
+#+end_note
+
+** Matlab Init                                              :noexport:ignore:
+#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
+  <<matlab-dir>>
+#+end_src
+
+#+begin_src matlab :exports none :results silent :noweb yes
+  <<matlab-init>>
+#+end_src
+
+** Load the processed data and the model
+#+begin_src matlab
+  load('./mat/spindle_data.mat', 'spindle');
+  load('./mat/spindle_model.mat', 'sys');
+#+end_src
+
+** Compute the PSD
+#+begin_src matlab :exports code :results none
+  n_av = 4; % Number of average
+
+  [pxx_1rpm, f_1rpm] = pwelch(spindle.rpm1.xasync, hanning(ceil(length(spindle.rpm1.xasync)/n_av)), [], [], spindle.rpm1.Fs);
+  [pyy_1rpm, ~]      = pwelch(spindle.rpm1.yasync, hanning(ceil(length(spindle.rpm1.yasync)/n_av)), [], [], spindle.rpm1.Fs);
+  [pzz_1rpm, ~]      = pwelch(spindle.rpm1.zasync, hanning(ceil(length(spindle.rpm1.zasync)/n_av)), [], [], spindle.rpm1.Fs);
+
+  [pxx_60rpm, f_60rpm] = pwelch(spindle.rpm60.xasync, hanning(ceil(length(spindle.rpm60.xasync)/n_av)), [], [], spindle.rpm60.Fs);
+  [pyy_60rpm, ~]       = pwelch(spindle.rpm60.yasync, hanning(ceil(length(spindle.rpm60.yasync)/n_av)), [], [], spindle.rpm60.Fs);
+  [pzz_60rpm, ~]       = pwelch(spindle.rpm60.zasync, hanning(ceil(length(spindle.rpm60.zasync)/n_av)), [], [], spindle.rpm60.Fs);
+#+end_src
+
+** Plot the computed PSD
+The Amplitude Spectral Densities of the displacement of the spindle for the $x$, $y$ and $z$ directions are shown figure [[fig:spindle_psd_xyz_60rpm]]. They correspond to the Asynchronous part shown figure [[fig:spindle_60rpm_sync_async]].
+
+#+begin_src matlab :exports code :results none
+  figure;
+  hold on;
+  plot(f_1rpm, (pxx_1rpm).^.5, 'DisplayName', '$P_{xx}$ - 1rpm');
+  plot(f_1rpm, (pyy_1rpm).^.5, 'DisplayName', '$P_{yy}$ - 1rpm');
+  plot(f_1rpm, (pzz_1rpm).^.5, 'DisplayName', '$P_{zz}$ - 1rpm');
+  % plot(f_1rpm, spindle.rpm1.adcn*ones(size(f_1rpm)), '--k', 'DisplayName', 'ADC - 1rpm');
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]');
+  legend('Location', 'northeast');
+  xlim([f_1rpm(2), f_1rpm(end)]);
+#+end_src
+#+NAME: fig:spindle_psd_xyz_1rpm
+#+HEADER: :tangle no :exports results :results raw :noweb yes
+#+begin_src matlab :var filepath="figs/spindle_psd_xyz_1rpm.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:spindle_psd_xyz_1rpm
+#+CAPTION: Power spectral density of the Asynchronous displacement - 1rpm
+#+RESULTS: fig:spindle_psd_xyz_1rpm
+[[file:figs/spindle_psd_xyz_1rpm.png]]
+
+#+begin_src matlab :exports code :results none
+  figure;
+  hold on;
+  plot(f_60rpm, (pxx_60rpm).^.5, 'DisplayName', '$P_{xx}$ - 60rpm');
+  plot(f_60rpm, (pyy_60rpm).^.5, 'DisplayName', '$P_{yy}$ - 60rpm');
+  plot(f_60rpm, (pzz_60rpm).^.5, 'DisplayName', '$P_{zz}$ - 60rpm');
+  % plot(f_60rpm, spindle.rpm60.adcn*ones(size(f_60rpm)), '--k', 'DisplayName', 'ADC - 60rpm');
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]');
+  legend('Location', 'northeast');
+  xlim([f_60rpm(2), f_60rpm(end)]);
+#+end_src
+#+NAME: fig:spindle_psd_xyz_60rpm
+#+HEADER: :tangle no :exports results :results raw :noweb yes
+#+begin_src matlab :var filepath="figs/spindle_psd_xyz_60rpm.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:spindle_psd_xyz_60rpm
+#+CAPTION: Power spectral density of the Asynchronous displacement - 60rpm
+#+RESULTS: fig:spindle_psd_xyz_60rpm
+[[file:figs/spindle_psd_xyz_60rpm.png]]
+
+** Compute the response of the model
+#+begin_src matlab :exports code :results none
+  Tfd = abs(squeeze(freqresp(sys('d', 'f'), f_60rpm, 'Hz')));
+#+end_src
+
+** Plot the PSD of the Force using the model
+#+begin_src matlab :exports code :results none
+  figure;
+  plot(f_60rpm, (pxx_60rpm.^.5)./Tfd, 'DisplayName', '$P_{xx}$');
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('ASD [$N/\sqrt{Hz}$]');
+  xlim([f_60rpm(2), f_60rpm(end)]);
+#+end_src
+#+NAME: fig:spindle_psd_f_60rpm
+#+HEADER: :tangle no :exports results :results raw :noweb yes
+#+begin_src matlab :var filepath="figs/spindle_psd_f_60rpm.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:spindle_psd_f_60rpm
+#+CAPTION: Power spectral density of the force - 60rpm
+#+RESULTS: fig:spindle_psd_f_60rpm
+[[file:figs/spindle_psd_f_60rpm.png]]
+
+** Estimated Shape of the PSD of the force
+#+begin_src matlab :exports code :results none
+  s = tf('s');
+
+  Wd_simple = 1e-8/(1+s/2/pi/0.5)/(1+s/2/pi/100);
+  Wf_simple = Wd_simple/tf(sys('d', 'f'));
+  TWf_simple = abs(squeeze(freqresp(Wf_simple, f_60rpm, 'Hz')));
+
+  % Wf = 0.48902*(s+327.9)*(s^2 + 109.6*s + 1.687e04)/((s^2 + 30.59*s + 8541)*(s^2 + 29.11*s + 3.268e04));
+  % Wf = 0.15788*(s+418.6)*(s+1697)^2*(s^2 + 124.3*s + 2.529e04)*(s^2 + 681.3*s + 9.018e05)/((s^2 + 23.03*s + 8916)*(s^2 + 33.85*s + 6.559e04)*(s^2 + 71.43*s + 4.283e05)*(s^2 + 40.64*s + 1.789e06));
+
+  Wf = (s+1697)^2*(s^2 + 114.5*s + 2.278e04)*(s^2 + 205.1*s + 1.627e05)*(s^2 + 285.8*s + 8.624e05)*(s+100)/((s+0.5)*3012*(s^2 + 23.03*s + 8916)*(s^2 + 17.07*s + 4.798e04)*(s^2 + 41.17*s + 4.347e05)*(s^2 + 78.99*s + 1.789e06));
+
+  TWf = abs(squeeze(freqresp(Wf, f_60rpm, 'Hz')));
+#+end_src
+
+** PSD in [N]
+Above $200Hz$, the Amplitude Spectral Density seems dominated by noise coming from the electronics (charge amplifier, ADC, ...). So we don't know what is the frequency content of the force above that frequency. However, we assume that $P_{xx}$ is decreasing with $1/f$ as it seems so be the case below $100Hz$ (figure [[fig:spindle_psd_xyz_60rpm]]).
+
+We then fit the PSD of the displacement with a transfer function (figure [[fig:spindle_psd_d_comp_60rpm]]).
+
+#+begin_src matlab :exports code :results none
+  figure;
+  hold on;
+  plot(f_60rpm, (pxx_60rpm.^.5)./Tfd, 'DisplayName', '$\sqrt{P_{xx}}/|T_{d/f}|$');
+  plot(f_60rpm, TWf, 'DisplayName', 'Wf');
+  plot(f_60rpm, TWf_simple, '-k', 'DisplayName', 'Wfs');
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('ASD [$N/\sqrt{Hz}$]');
+  xlim([f_60rpm(2), f_60rpm(end)]);
+  legend('Location', 'northeast');
+#+end_src
+#+NAME: fig:spindle_psd_f_comp_60rpm
+#+HEADER: :tangle no :exports results :results raw :noweb yes
+#+begin_src matlab :var filepath="figs/spindle_psd_f_comp_60rpm.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:spindle_psd_f_comp_60rpm
+#+CAPTION: Power spectral density of the force - 60rpm
+#+RESULTS: fig:spindle_psd_f_comp_60rpm
+[[file:figs/spindle_psd_f_comp_60rpm.png]]
+
+** PSD in [m]
+To obtain the PSD of the force $f$ that induce such displacement, we use the following formula:
+\[ \sqrt{PSD(d)} = |T_{d/f}| \sqrt{PSD(f)} \]
+
+And so we have:
+\[ \sqrt{PSD(f)} = |T_{d/f}|^{-1} \sqrt{PSD(d)} \]
+
+The obtain Power Spectral Density of the force is displayed figure [[fig:spindle_psd_f_comp_60rpm]].
+
+#+begin_src matlab :exports code :results none
+  figure;
+  hold on;
+  plot(f_60rpm, pxx_60rpm.^.5, 'DisplayName', '$\sqrt{P_{xx}}$');
+  plot(f_60rpm, TWf.*Tfd, 'DisplayName', '$|W_f|*|T_{d/f}|$');
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]');
+  xlim([f_60rpm(2), f_60rpm(end)]);
+  legend('Location', 'northeast');
+#+end_src
+
+#+NAME: fig:spindle_psd_d_comp_60rpm
+#+HEADER: :tangle no :exports results :results raw :noweb yes
+#+begin_src matlab :var filepath="figs/spindle_psd_d_comp_60rpm.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:spindle_psd_d_comp_60rpm
+#+CAPTION: Comparison of the power spectral density of the measured displacement and of the model
+#+RESULTS: fig:spindle_psd_d_comp_60rpm
+[[file:figs/spindle_psd_d_comp_60rpm.png]]
+
+** Compute the resulting RMS value [m]
+#+begin_src matlab :exports code :results none
+  figure;
+  hold on;
+  plot(f_60rpm, 1e9*cumtrapz(f_60rpm, (pxx_60rpm)).^.5, '--', 'DisplayName', 'Exp. Data');
+  plot(f_60rpm, 1e9*cumtrapz(f_60rpm, ((TWf.*Tfd).^2)).^.5, '--', 'DisplayName', 'Estimated');
+  hold off;
+  set(gca, 'XScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('CPS [$nm$ rms]');
+  xlim([f_60rpm(2), f_60rpm(end)]);
+  legend('Location', 'southeast');
+#+end_src
+
+#+NAME: fig:spindle_cps_d_comp_60rpm
+#+HEADER: :tangle no :exports results :results raw :noweb yes
+#+begin_src matlab :var filepath="figs/spindle_cps_d_comp_60rpm.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:spindle_cps_d_comp_60rpm
+#+CAPTION: Cumulative Power Spectrum - 60rpm
+#+RESULTS: fig:spindle_cps_d_comp_60rpm
+[[file:figs/spindle_cps_d_comp_60rpm.png]]
+
+** Compute the resulting RMS value [m]
+#+begin_src matlab :exports code :results none
+  figure;
+  hold on;
+  plot(f_1rpm, 1e9*cumtrapz(f_1rpm, (pxx_1rpm)), '--', 'DisplayName', 'Exp. Data');
+  plot(f_1rpm, 1e9*(f_1rpm(end)-f_1rpm(1))/(length(f_1rpm)-1)*cumsum(pxx_1rpm), '--', 'DisplayName', 'Exp. Data');
+  hold off;
+  set(gca, 'XScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('CPS [$nm$ rms]');
+  xlim([f_1rpm(2), f_1rpm(end)]);
+  legend('Location', 'southeast');
+#+end_src
+
+#+NAME: fig:spindle_cps_d_comp_1rpm
+#+HEADER: :tangle no :exports results :results raw :noweb yes
+#+begin_src matlab :var filepath="figs/spindle_cps_d_comp_1rpm.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:spindle_cps_d_comp_1rpm
+#+CAPTION: Cumulative Power Spectrum - 1rpm
+#+RESULTS: fig:spindle_cps_d_comp_1rpm
+[[file:figs/spindle_cps_d_comp_1rpm.png]]
+
+* Functions
+** getAsynchronousError
+  :PROPERTIES:
+  :header-args:matlab+: :tangle src/getAsynchronousError.m
+  :header-args:matlab+: :comments org :mkdirp yes
+  :END:
+  <<sec:getAsynchronousError>>
+
+This Matlab function is accessible [[file:src/getAsynchronousError.m][here]].
+
+#+begin_src matlab :eval no :results none
+  function [Wxdec] = getAsynchronousError(data, NbTurn)
+  %%
+      L = length(data);
+      res_per_rev = L/NbTurn;
+
+      P = 0:(res_per_rev*NbTurn-1);
+      Pos = P' * 360/res_per_rev;
+
+      % Temperature correction
+      x1 = myfit2(Pos, data);
+
+      % Convert data to frequency domain and scale accordingly
+      X2 = 2/(res_per_rev*NbTurn)*fft(x1);
+      f2 = (0:L-1)./NbTurn; %upr -> once per revolution
+
+      %%
+      X2dec = zeros(size(X2));
+      % Get only the non integer data
+      X2dec(mod(f2(:), 1) ~= 0) = X2(mod(f2(:), 1) ~= 0);
+
+      Wxdec = real((res_per_rev*NbTurn)/2 * ifft(X2dec));
+
+      %%
+      function Y =  myfit2(x,y)
+          A = [x ones(size(x))]\y;
+          a = A(1); b = A(2);
+          Y = y - (a*x + b);
+      end
+  end
+#+end_src

spindle/matlab/spindle_data_processing.m

@@ -0,0 +1,86 @@
+%% Clear Workspace and Close figures
+clear; close all; clc;
+
+%% Intialize Laplace variable
+s = zpk('s');
+
+%% Add the getAsynchronousError to path
+addpath('./src/');
+
+% Load Measurement Data
+
+spindle_1rpm_table  = readtable('./mat/10turns_1rpm_icepap.txt');
+spindle_60rpm_table = readtable('./mat/10turns_60rpm_IcepapFIR.txt');
+
+spindle_1rpm_table(1, :)
+
+spindle_1rpm  = table2array(spindle_1rpm_table);
+spindle_60rpm = table2array(spindle_60rpm_table);
+
+% Convert Signals from [deg] to [sec]
+
+speed_1rpm = 360/60; % [deg/sec]
+spindle_1rpm(:, 1) = spindle_1rpm(:, 2)/speed_1rpm;  % From position [deg] to time [s]
+
+speed_60rpm = 360/1; % [deg/sec]
+spindle_60rpm(:, 1) = spindle_60rpm(:, 2)/speed_60rpm; % From position [deg] to time [s]
+
+% Convert Signals
+
+% scaling = 1/80000; % 80 mV/um
+scaling = 1e-6; % [um] to [m]
+
+spindle_1rpm(:, 3:end) = scaling*spindle_1rpm(:, 3:end); % [V] to [m]
+spindle_1rpm(:, 3:end) = spindle_1rpm(:, 3:end)-mean(spindle_1rpm(:, 3:end)); % Remove mean
+
+spindle_60rpm(:, 3:end) = scaling*spindle_60rpm(:, 3:end); % [V] to [m]
+spindle_60rpm(:, 3:end) = spindle_60rpm(:, 3:end)-mean(spindle_60rpm(:, 3:end)); % Remove mean
+
+% Ts and Fs for both measurements
+
+Ts_1rpm = spindle_1rpm(end, 1)/(length(spindle_1rpm(:, 1))-1);
+Fs_1rpm = 1/Ts_1rpm;
+
+Ts_60rpm = spindle_60rpm(end, 1)/(length(spindle_60rpm(:, 1))-1);
+Fs_60rpm = 1/Ts_60rpm;
+
+% Find Noise of the ADC [$\frac{m}{\sqrt{Hz}}$]
+
+data = spindle_1rpm(:, 5);
+dV_1rpm = min(abs(data(1) - data(data ~= data(1))));
+noise_1rpm = dV_1rpm/sqrt(12*Fs_1rpm/2);
+
+data = spindle_60rpm(:, 5);
+dV_60rpm = min(abs(data(50) - data(data ~= data(50))));
+noise_60rpm = dV_60rpm/sqrt(12*Fs_60rpm/2);
+
+% Save all the data under spindle struct
+
+spindle.rpm1.time = spindle_1rpm(:, 1);
+spindle.rpm1.deg  = spindle_1rpm(:, 2);
+spindle.rpm1.Ts   = Ts_1rpm;
+spindle.rpm1.Fs   = 1/Ts_1rpm;
+spindle.rpm1.x    = spindle_1rpm(:, 3);
+spindle.rpm1.y    = spindle_1rpm(:, 4);
+spindle.rpm1.z    = spindle_1rpm(:, 5);
+spindle.rpm1.adcn = noise_1rpm;
+
+spindle.rpm60.time = spindle_60rpm(:, 1);
+spindle.rpm60.deg  = spindle_60rpm(:, 2);
+spindle.rpm60.Ts   = Ts_60rpm;
+spindle.rpm60.Fs   = 1/Ts_60rpm;
+spindle.rpm60.x    = spindle_60rpm(:, 3);
+spindle.rpm60.y    = spindle_60rpm(:, 4);
+spindle.rpm60.z    = spindle_60rpm(:, 5);
+spindle.rpm60.adcn = noise_60rpm;
+
+% Compute Asynchronous data
+
+for direction = {'x', 'y', 'z'}
+    spindle.rpm1.([direction{1}, 'async']) = getAsynchronousError(spindle.rpm1.(direction{1}), 10);
+    spindle.rpm60.([direction{1}, 'async']) = getAsynchronousError(spindle.rpm60.(direction{1}), 10);
+end
+
+% Save data
+
+save('./mat/spindle_data.mat', 'spindle');

spindle/matlab/spindle_model.m

@@ -0,0 +1,103 @@
+%% Clear Workspace and Close figures
+clear; close all; clc;
+
+%% Intialize Laplace variable
+s = zpk('s');
+
+%% Add the getAsynchronousError to path
+addpath('./src/');
+
+% Parameters
+
+mg = 3000; % Mass of granite [kg]
+ms = 50;   % Mass of Spindle [kg]
+
+kg = 1e8; % Stiffness of granite [N/m]
+ks = 5e7; % Stiffness of spindle [N/m]
+
+% Compute Mass and Stiffness Matrices
+
+Mm = diag([ms, mg]);
+Km = diag([ks, ks+kg]) - diag(ks, -1) - diag(ks, 1);
+
+% Compute resonance frequencies
+
+A = [zeros(size(Mm)) eye(size(Mm)) ; -Mm\Km zeros(size(Mm))];
+eigA = imag(eigs(A))/2/pi;
+eigA = eigA(eigA>0);
+eigA = eigA(1:2);
+
+% From model_damping compute the Damping Matrix
+
+modal_damping = 1e-5;
+
+ab = [0.5*eigA(1) 0.5/eigA(1) ; 0.5*eigA(2) 0.5/eigA(2)]\[modal_damping ; modal_damping];
+
+Cm = ab(1)*Mm +ab(2)*Km;
+
+% Define inputs, outputs and state names
+
+StateName = {...
+    'xs', ... % Displacement of Spindle [m]
+    'xg', ... % Displacement of Granite [m]
+    'vs', ... % Velocity of Spindle [m]
+    'vg', ... % Velocity of Granite [m]
+            };
+StateUnit = {'m', 'm', 'm/s', 'm/s'};
+
+InputName = {...
+    'f' ... % Spindle Force [N]
+            };
+InputUnit = {'N'};
+
+OutputName = {...
+    'd' ... % Displacement [m]
+             };
+OutputUnit = {'m'};
+
+% Define A, B and C matrices
+
+% A Matrix
+A = [zeros(size(Mm)) eye(size(Mm)) ; ...
+     -Mm\Km          -Mm\Cm];
+
+% B Matrix
+B_low = zeros(length(StateName), length(InputName));
+B_low(strcmp(StateName,'vs'), strcmp(InputName,'f')) =  1;
+B_low(strcmp(StateName,'vg'), strcmp(InputName,'f')) = -1;
+B = blkdiag(zeros(length(StateName)/2), pinv(Mm))*B_low;
+
+% C Matrix
+C = zeros(length(OutputName), length(StateName));
+C(strcmp(OutputName,'d'), strcmp(StateName,'xs')) =  1;
+C(strcmp(OutputName,'d'), strcmp(StateName,'xg')) = -1;
+
+% D Matrix
+D = zeros(length(OutputName), length(InputName));
+
+% Generate the State Space Model
+
+sys = ss(A, B, C, D);
+sys.StateName = StateName;
+sys.StateUnit = StateUnit;
+
+sys.InputName = InputName;
+sys.InputUnit = InputUnit;
+
+sys.OutputName = OutputName;
+sys.OutputUnit = OutputUnit;
+
+% Bode Plot
+% The transfer function from a disturbance force $f$ to the measured displacement $d$ is shown figure [[fig:spindle_f_to_d]].
+
+
+freqs = logspace(-1, 3, 1000);
+
+figure;
+plot(freqs, abs(squeeze(freqresp(sys('d', 'f'), freqs, 'Hz'))));
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Amplitude [m/N]');
+
+% Save the model
+
+save('./mat/spindle_model.mat', 'sys');

spindle/matlab/spindle_psd.m

@@ -0,0 +1,144 @@
+%% Clear Workspace and Close figures
+clear; close all; clc;
+
+%% Intialize Laplace variable
+s = zpk('s');
+
+% Load the processed data and the model
+
+load('./mat/spindle_data.mat', 'spindle');
+load('./mat/spindle_model.mat', 'sys');
+
+% Compute the PSD
+
+n_av = 4; % Number of average
+
+[pxx_1rpm, f_1rpm] = pwelch(spindle.rpm1.xasync, hanning(ceil(length(spindle.rpm1.xasync)/n_av)), [], [], spindle.rpm1.Fs);
+[pyy_1rpm, ~]      = pwelch(spindle.rpm1.yasync, hanning(ceil(length(spindle.rpm1.yasync)/n_av)), [], [], spindle.rpm1.Fs);
+[pzz_1rpm, ~]      = pwelch(spindle.rpm1.zasync, hanning(ceil(length(spindle.rpm1.zasync)/n_av)), [], [], spindle.rpm1.Fs);
+
+[pxx_60rpm, f_60rpm] = pwelch(spindle.rpm60.xasync, hanning(ceil(length(spindle.rpm60.xasync)/n_av)), [], [], spindle.rpm60.Fs);
+[pyy_60rpm, ~]       = pwelch(spindle.rpm60.yasync, hanning(ceil(length(spindle.rpm60.yasync)/n_av)), [], [], spindle.rpm60.Fs);
+[pzz_60rpm, ~]       = pwelch(spindle.rpm60.zasync, hanning(ceil(length(spindle.rpm60.zasync)/n_av)), [], [], spindle.rpm60.Fs);
+
+% Plot the computed PSD
+% The Amplitude Spectral Densities of the displacement of the spindle for the $x$, $y$ and $z$ directions are shown figure [[fig:spindle_psd_xyz_60rpm]]. They correspond to the Asynchronous part shown figure [[fig:spindle_60rpm_sync_async]].
+
+
+figure;
+hold on;
+plot(f_1rpm, (pxx_1rpm).^.5, 'DisplayName', '$P_{xx}$ - 1rpm');
+plot(f_1rpm, (pyy_1rpm).^.5, 'DisplayName', '$P_{yy}$ - 1rpm');
+plot(f_1rpm, (pzz_1rpm).^.5, 'DisplayName', '$P_{zz}$ - 1rpm');
+% plot(f_1rpm, spindle.rpm1.adcn*ones(size(f_1rpm)), '--k', 'DisplayName', 'ADC - 1rpm');
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]');
+legend('Location', 'northeast');
+xlim([f_1rpm(2), f_1rpm(end)]);
+
+
+
+% #+NAME: fig:spindle_psd_xyz_1rpm
+% #+CAPTION: Power spectral density of the Asynchronous displacement - 1rpm
+% #+RESULTS: fig:spindle_psd_xyz_1rpm
+% [[file:figs/spindle_psd_xyz_1rpm.png]]
+
+
+figure;
+hold on;
+plot(f_60rpm, (pxx_60rpm).^.5, 'DisplayName', '$P_{xx}$ - 60rpm');
+plot(f_60rpm, (pyy_60rpm).^.5, 'DisplayName', '$P_{yy}$ - 60rpm');
+plot(f_60rpm, (pzz_60rpm).^.5, 'DisplayName', '$P_{zz}$ - 60rpm');
+% plot(f_60rpm, spindle.rpm60.adcn*ones(size(f_60rpm)), '--k', 'DisplayName', 'ADC - 60rpm');
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]');
+legend('Location', 'northeast');
+xlim([f_60rpm(2), f_60rpm(end)]);
+
+% Compute the response of the model
+
+Tfd = abs(squeeze(freqresp(sys('d', 'f'), f_60rpm, 'Hz')));
+
+% Plot the PSD of the Force using the model
+
+figure;
+plot(f_60rpm, (pxx_60rpm.^.5)./Tfd, 'DisplayName', '$P_{xx}$');
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]'); ylabel('ASD [$N/\sqrt{Hz}$]');
+xlim([f_60rpm(2), f_60rpm(end)]);
+
+% Estimated Shape of the PSD of the force
+
+s = tf('s');
+
+Wd_simple = 1e-8/(1+s/2/pi/0.5)/(1+s/2/pi/100);
+Wf_simple = Wd_simple/tf(sys('d', 'f'));
+TWf_simple = abs(squeeze(freqresp(Wf_simple, f_60rpm, 'Hz')));
+
+% Wf = 0.48902*(s+327.9)*(s^2 + 109.6*s + 1.687e04)/((s^2 + 30.59*s + 8541)*(s^2 + 29.11*s + 3.268e04));
+% Wf = 0.15788*(s+418.6)*(s+1697)^2*(s^2 + 124.3*s + 2.529e04)*(s^2 + 681.3*s + 9.018e05)/((s^2 + 23.03*s + 8916)*(s^2 + 33.85*s + 6.559e04)*(s^2 + 71.43*s + 4.283e05)*(s^2 + 40.64*s + 1.789e06));
+
+Wf = (s+1697)^2*(s^2 + 114.5*s + 2.278e04)*(s^2 + 205.1*s + 1.627e05)*(s^2 + 285.8*s + 8.624e05)*(s+100)/((s+0.5)*3012*(s^2 + 23.03*s + 8916)*(s^2 + 17.07*s + 4.798e04)*(s^2 + 41.17*s + 4.347e05)*(s^2 + 78.99*s + 1.789e06));
+
+TWf = abs(squeeze(freqresp(Wf, f_60rpm, 'Hz')));
+
+% PSD in [N]
+% Above $200Hz$, the Amplitude Spectral Density seems dominated by noise coming from the electronics (charge amplifier, ADC, ...). So we don't know what is the frequency content of the force above that frequency. However, we assume that $P_{xx}$ is decreasing with $1/f$ as it seems so be the case below $100Hz$ (figure [[fig:spindle_psd_xyz_60rpm]]).
+
+% We then fit the PSD of the displacement with a transfer function (figure [[fig:spindle_psd_d_comp_60rpm]]).
+
+
+figure;
+hold on;
+plot(f_60rpm, (pxx_60rpm.^.5)./Tfd, 'DisplayName', '$\sqrt{P_{xx}}/|T_{d/f}|$');
+plot(f_60rpm, TWf, 'DisplayName', 'Wf');
+plot(f_60rpm, TWf_simple, '-k', 'DisplayName', 'Wfs');
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]'); ylabel('ASD [$N/\sqrt{Hz}$]');
+xlim([f_60rpm(2), f_60rpm(end)]);
+legend('Location', 'northeast');
+
+% PSD in [m]
+% To obtain the PSD of the force $f$ that induce such displacement, we use the following formula:
+% \[ \sqrt{PSD(d)} = |T_{d/f}| \sqrt{PSD(f)} \]
+
+% And so we have:
+% \[ \sqrt{PSD(f)} = |T_{d/f}|^{-1} \sqrt{PSD(d)} \]
+
+% The obtain Power Spectral Density of the force is displayed figure [[fig:spindle_psd_f_comp_60rpm]].
+
+
+figure;
+hold on;
+plot(f_60rpm, pxx_60rpm.^.5, 'DisplayName', '$\sqrt{P_{xx}}$');
+plot(f_60rpm, TWf.*Tfd, 'DisplayName', '$|W_f|*|T_{d/f}|$');
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]');
+xlim([f_60rpm(2), f_60rpm(end)]);
+legend('Location', 'northeast');
+
+% Compute the resulting RMS value [m]
+
+figure;
+hold on;
+plot(f_60rpm, 1e9*cumtrapz(f_60rpm, (pxx_60rpm)).^.5, '--', 'DisplayName', 'Exp. Data');
+plot(f_60rpm, 1e9*cumtrapz(f_60rpm, ((TWf.*Tfd).^2)).^.5, '--', 'DisplayName', 'Estimated');
+hold off;
+set(gca, 'XScale', 'log');
+xlabel('Frequency [Hz]'); ylabel('CPS [$nm$ rms]');
+xlim([f_60rpm(2), f_60rpm(end)]);
+legend('Location', 'southeast');
+
+% Compute the resulting RMS value [m]
+
+figure;
+hold on;
+plot(f_1rpm, 1e9*cumtrapz(f_1rpm, (pxx_1rpm)), '--', 'DisplayName', 'Exp. Data');
+plot(f_1rpm, 1e9*(f_1rpm(end)-f_1rpm(1))/(length(f_1rpm)-1)*cumsum(pxx_1rpm), '--', 'DisplayName', 'Exp. Data');
+hold off;
+set(gca, 'XScale', 'log');
+xlabel('Frequency [Hz]'); ylabel('CPS [$nm$ rms]');
+xlim([f_1rpm(2), f_1rpm(end)]);
+legend('Location', 'southeast');

spindle/matlab/spindle_time_domain.m

@@ -0,0 +1,75 @@
+%% Clear Workspace and Close figures
+clear; close all; clc;
+
+%% Intialize Laplace variable
+s = zpk('s');
+
+% Load the processed data
+
+load('./mat/spindle_data.mat', 'spindle');
+
+% Plot X-Y-Z position with respect to Time - 1rpm
+
+figure;
+hold on;
+plot(spindle.rpm1.time, spindle.rpm1.x);
+plot(spindle.rpm1.time, spindle.rpm1.y);
+plot(spindle.rpm1.time, spindle.rpm1.z);
+hold off;
+xlabel('Time [s]'); ylabel('Amplitude [m]');
+legend({'tx - 1rpm', 'ty - 1rpm', 'tz - 1rpm'});
+
+% Plot X-Y-Z position with respect to Time - 60rpm
+% The measurements for the spindle turning at 60rpm are shown figure [[fig:spindle_xyz_60rpm]].
+
+
+figure;
+hold on;
+plot(spindle.rpm60.time, spindle.rpm60.x);
+plot(spindle.rpm60.time, spindle.rpm60.y);
+plot(spindle.rpm60.time, spindle.rpm60.z);
+hold off;
+xlabel('Time [s]'); ylabel('Amplitude [m]');
+legend({'tx - 60rpm', 'ty - 60rpm', 'tz - 60rpm'});
+
+% Plot Synchronous and Asynchronous - 1rpm
+
+figure;
+hold on;
+plot(spindle.rpm1.time, spindle.rpm1.x);
+plot(spindle.rpm1.time, spindle.rpm1.xasync);
+hold off;
+xlabel('Time [s]'); ylabel('Amplitude [m]');
+legend({'tx - 1rpm - Sync', 'tx - 1rpm - Async'});
+
+% Plot Synchronous and Asynchronous - 60rpm
+% The data is split into its Synchronous and Asynchronous part (figure [[fig:spindle_60rpm_sync_async]]). We then use the Asynchronous part for the analysis in the following sections as we suppose that we can deal with the synchronous part with feedforward control.
+
+
+figure;
+hold on;
+plot(spindle.rpm60.time, spindle.rpm60.x);
+plot(spindle.rpm60.time, spindle.rpm60.xasync);
+hold off;
+xlabel('Time [s]'); ylabel('Amplitude [m]');
+legend({'tx - 60rpm - Sync', 'tx - 60rpm - Async'});
+
+% Plot X against Y
+
+figure;
+hold on;
+plot(spindle.rpm1.x, spindle.rpm1.y);
+plot(spindle.rpm60.x, spindle.rpm60.y);
+hold off;
+xlabel('X Amplitude [m]'); ylabel('Y Amplitude [m]');
+legend({'1rpm', '60rpm'});
+
+% Plot X against Y - Asynchronous
+
+figure;
+hold on;
+plot(spindle.rpm1.xasync, spindle.rpm1.yasync);
+plot(spindle.rpm60.xasync, spindle.rpm60.yasync);
+hold off;
+xlabel('X Amplitude [m]'); ylabel('Y Amplitude [m]');
+legend({'1rpm', '60rpm'});

spindle/spindle_measurements.m

@@ -0,0 +1,72 @@
+%%
+clear; close all; clc;
+
+%% Load Measurement Data
+spindle_1rpm_table  = readtable('../../Measurements/Spindle/10turns_1rpm_icepap.txt');
+spindle_60rpm_table = readtable('../../Measurements/Spindle/10turns_60rpm_IcepapFIR.txt');
+
+disp(spindle_1rpm_table(1, :));
+
+spindle_1rpm  = table2array(spindle_1rpm_table);
+spindle_60rpm = table2array(spindle_60rpm_table);
+
+%% Convert Signals from [deg] to [sec]
+speed_1rpm = 360/60; % [deg/sec]
+spindle_1rpm(:, 1) = spindle_1rpm(:, 2)/speed_1rpm;  % From position [deg] to time [s]
+
+speed_60rpm = 360/1; % [deg/sec]
+spindle_60rpm(:, 1) = spindle_60rpm(:, 2)/speed_60rpm; % From position [deg] to time [s]
+
+%% Convert Signals
+% scaling = 1/80000; % 80 mV/um
+scaling = 1e-6; % [um] to [m]
+
+spindle_1rpm(:, 3:end) = scaling*spindle_1rpm(:, 3:end); % [V] to [m]
+spindle_1rpm(:, 3:end) = spindle_1rpm(:, 3:end)-mean(spindle_1rpm(:, 3:end)); % Remove mean
+
+spindle_60rpm(:, 3:end) = scaling*spindle_60rpm(:, 3:end); % [V] to [m]
+spindle_60rpm(:, 3:end) = spindle_60rpm(:, 3:end)-mean(spindle_60rpm(:, 3:end)); % Remove mean
+
+%% Ts and Fs for both measurements
+Ts_1rpm = spindle_1rpm(end, 1)/(length(spindle_1rpm(:, 1))-1);
+Fs_1rpm = 1/Ts_1rpm;
+
+Ts_60rpm = spindle_60rpm(end, 1)/(length(spindle_60rpm(:, 1))-1);
+Fs_60rpm = 1/Ts_60rpm;
+
+%% Find Noise of the ADC [m/sqrt(Hz)]
+data = spindle_1rpm(:, 5);
+dV_1rpm = min(abs(data(1) - data(data ~= data(1))));
+noise_1rpm = dV_1rpm/sqrt(12*Fs_1rpm/2);
+
+data = spindle_60rpm(:, 5);
+dV_60rpm = min(abs(data(50) - data(data ~= data(50))));
+noise_60rpm = dV_60rpm/sqrt(12*Fs_60rpm/2);
+
+%% Save all the data under spindle struct
+spindle.rpm1.time = spindle_1rpm(:, 1);
+spindle.rpm1.deg  = spindle_1rpm(:, 2);
+spindle.rpm1.Ts   = Ts_1rpm;
+spindle.rpm1.Fs   = 1/Ts_1rpm;
+spindle.rpm1.x    = spindle_1rpm(:, 3);
+spindle.rpm1.y    = spindle_1rpm(:, 4);
+spindle.rpm1.z    = spindle_1rpm(:, 5);
+spindle.rpm1.adcn = noise_1rpm;
+
+spindle.rpm60.time = spindle_60rpm(:, 1);
+spindle.rpm60.deg  = spindle_60rpm(:, 2);
+spindle.rpm60.Ts   = Ts_60rpm;
+spindle.rpm60.Fs   = 1/Ts_60rpm;
+spindle.rpm60.x    = spindle_60rpm(:, 3);
+spindle.rpm60.y    = spindle_60rpm(:, 4);
+spindle.rpm60.z    = spindle_60rpm(:, 5);
+spindle.rpm60.adcn = noise_60rpm;
+
+%% Compute Asynchronous data
+for direction = {'x', 'y', 'z'}
+    spindle.rpm1.([direction{1}, 'async']) = getAsynchronousError(spindle.rpm1.(direction{1}), 10);
+    spindle.rpm60.([direction{1}, 'async']) = getAsynchronousError(spindle.rpm60.(direction{1}), 10);
+end
+
+%% Save data
+save('./mat/spindle_data.mat', 'spindle');

spindle/src/getAsynchronousError.m

@@ -0,0 +1,39 @@
+% getAsynchronousError
+%   :PROPERTIES:
+%   :header-args:matlab+: :tangle src/getAsynchronousError.m
+%   :header-args:matlab+: :comments org :mkdirp yes
+%   :END:
+%   <<sec:getAsynchronousError>>
+
+% This Matlab function is accessible [[file:src/getAsynchronousError.m][here]].
+
+
+function [Wxdec] = getAsynchronousError(data, NbTurn)
+%%
+    L = length(data);
+    res_per_rev = L/NbTurn;
+
+    P = 0:(res_per_rev*NbTurn-1);
+    Pos = P' * 360/res_per_rev;
+
+    % Temperature correction
+    x1 = myfit2(Pos, data);
+
+    % Convert data to frequency domain and scale accordingly
+    X2 = 2/(res_per_rev*NbTurn)*fft(x1);
+    f2 = (0:L-1)./NbTurn; %upr -> once per revolution
+
+    %%
+    X2dec = zeros(size(X2));
+    % Get only the non integer data
+    X2dec(mod(f2(:), 1) ~= 0) = X2(mod(f2(:), 1) ~= 0);
+
+    Wxdec = real((res_per_rev*NbTurn)/2 * ifft(X2dec));
+
+    %%
+    function Y =  myfit2(x,y)
+        A = [x ones(size(x))]\y;
+        a = A(1); b = A(2);
+        Y = y - (a*x + b);
+    end
+end