Add every computation files and data

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/data/10turns_60rmp_Icepap-FIR/res.txt

@@ -0,0 +1,110 @@
+Error X et Y
+
+LSCxtot =
+
+    0.3872
+
+
+LSCytot =
+
+    0.2419
+
+
+LSCxsync =
+
+    0.1577
+
+
+LSCysync =
+
+    0.1602
+
+
+LSCxasync =
+
+    0.2946
+
+
+LSCyasync =
+
+    0.1103
+
+
+LSCxytot =
+
+    0.3519
+
+
+Error X2 et Y2
+
+
+LSCxtot =
+
+    0.3354
+
+
+LSCytot =
+
+    0.3202
+
+
+LSCxsync =
+
+    0.1101
+
+
+LSCysync =
+
+    0.0808
+
+
+LSCxasync =
+
+    0.2588
+
+
+LSCyasync =
+
+    0.2791
+
+
+LSCxytot =
+
+    0.3642
+
+
+
+ErrorZ
+LSCxtot =
+
+    0.0775
+
+
+LSCytot =
+
+    0.0775
+
+
+LSCxsync =
+
+    0.0390
+
+
+LSCysync =
+
+    0.0390
+
+
+LSCxasync =
+
+    0.0617
+
+
+LSCyasync =
+
+    0.0617
+
+
+LSCxytot =
+
+    0.1000

Spindle/getAsynchronousError.m

@@ -0,0 +1,31 @@
+% getAsynchronousError
+
+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

Spindle/index.org

@@ -0,0 +1,621 @@
+#+TITLE: Spindle Analysis
+:drawer:
+#+STARTUP: overview
+
+#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/htmlize.css"/>
+#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/readtheorg.css"/>
+#+HTML_HEAD: <script src="../js/jquery.min.js"></script>
+#+HTML_HEAD: <script src="../js/bootstrap.min.js"></script>
+#+HTML_HEAD: <script src="../js/jquery.stickytableheaders.min.js"></script>
+#+HTML_HEAD: <script src="../js/readtheorg.js"></script>
+
+#+LATEX_CLASS: cleanreport
+#+LaTeX_CLASS_OPTIONS: [tocnp, secbreak, minted]
+
+#+PROPERTY: header-args:matlab  :session *MATLAB*
+#+PROPERTY: header-args:matlab+ :comments org
+#+PROPERTY: header-args:matlab+ :exports both
+#+PROPERTY: header-args:matlab+ :eval no-export
+#+PROPERTY: header-args:matlab+ :noweb yes
+#+PROPERTY: header-args:matlab+ :mkdirp yes
+#+PROPERTY: header-args:matlab+ :output-dir figs
+:end:
+
+* Data Processing
+#+begin_src matlab :results none :exports none
+  <<matlab-init>>
+#+end_src
+
+** Load Measurement Data
+#+begin_src matlab :results none :exports code
+  spindle_1rpm_table  = readtable('./data/10turns_1rpm_icepap.txt');
+  spindle_60rpm_table = readtable('./data/10turns_60rpm_IcepapFIR.txt');
+#+end_src
+
+#+begin_src matlab :results output :exports code
+  spindle_1rpm_table(1, :)
+#+end_src
+
+#+RESULTS:
+: spindle_1rpm_table(1, :)
+: {Undefined function or variable 'spindle_1rpm_table'.}
+
+#+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 [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
+
+* Time Domain Data
+** 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
+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]]
+
+** 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]]
+
+* Power Spectral Density
+** 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]]
+
+** Load the model of the spindle
+#+begin_src matlab :exports code :results none
+  load('./mat/spindle_model.mat', 'sys');
+  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
+#+begin_src matlab :results none :eval no :tangle getAsynchronousError.m
+  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/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/spindle_model.m

@@ -0,0 +1,88 @@
+%%
+clear; close all; clc;
+
+%% 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;
+
+%%
+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'};
+
+%% Generate the State Space Model
+% 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));
+
+%%
+sys = ss(A, B, C, D);
+sys.StateName = StateName;
+sys.StateUnit = StateUnit;
+
+sys.InputName = InputName;
+sys.InputUnit = InputUnit;
+
+sys.OutputName = OutputName;
+sys.OutputUnit = OutputUnit;
+
+%%
+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]');
+
+exportFig('spindle_f_to_d', 'normal-normal');
+
+%% Save the model
+save('./mat/spindle_model.mat', 'sys');

Spindle/spindle_psd.m

@@ -0,0 +1,130 @@
+%%
+clear; close all; clc;
+
+%% Load Spindle Data
+load('./mat/spindle_data.mat', 'spindle');
+
+%% 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);
+
+%%
+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)]);
+
+exportFig('spindle_psd_xyz_1rpm', 'normal-normal');
+
+%%
+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)]);
+
+exportFig('spindle_psd_xyz_60rpm', 'wide-normal');
+
+%% Load the model of the spindle
+load('./mat/spindle_model.mat', 'sys');
+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)]);
+
+exportFig('spindle_psd_f_60rpm', 'wide-normal');
+
+%% 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));
+
+% sisotool(tf(1), Wf) % SISOTOOL
+
+TWf = abs(squeeze(freqresp(Wf, f_60rpm, 'Hz')));
+
+%% PSD in [N]
+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');
+
+exportFig('spindle_psd_f_comp_60rpm', 'wide-normal');
+
+%% PSD in [m]
+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');
+
+exportFig('spindle_psd_d_comp_60rpm', 'normal-normal');
+
+%% 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('ASD [$nm/\sqrt{Hz}$]');
+xlim([f_60rpm(2), f_60rpm(end)]);
+legend('Location', 'northeast');
+
+exportFig('spindle_cps_d_comp_60rpm', 'normal-normal');
+
+%% 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('ASD [$nm/\sqrt{Hz}$]');
+xlim([f_1rpm(2), f_1rpm(end)]);
+legend('Location', 'northeast');
+
+exportFig('spindle_cps_d_comp_1rpm', 'normal-normal');
+
+%% Save
+save('./mat/spindle_Wf.mat', 'Wf');
+

Spindle/spindle_time_domain.m

@@ -0,0 +1,73 @@
+%%
+clear; close all; clc;
+
+%% Load Spindle Data
+load('./mat/spindle_data.mat', 'spindle');
+
+%% Plot X-Y-Z position with respect to Time
+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'});
+
+exportFig('spindle_xyz_1rpm', 'wide-normal');
+
+%% Plot X-Y-Z position with respect to Time
+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'});
+
+exportFig('spindle_xyz_60rpm', 'normal-normal');
+
+%%
+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'});
+
+exportFig('spindle_1rpm_sync_async', 'normal-normal');
+
+%%
+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'});
+
+exportFig('spindle_60rpm_sync_async', 'wide-normal');
+
+%% 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'});
+
+exportFig('spindle_xy_1_60rpm', 'normal-normal');
+
+%% Plot X against Y
+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'});
+
+exportFig('spindle_xy_1_60rpm_async', 'normal-normal');