function lambda = stabilization(u, dt, Nlim, omlim, tol, omvec, indi)
%STABILIZATION Experimental Modal Analysis Stabilization process
% LAMBDA = STABILIZATION(IR, DT, NLIM, OMLIM, TOL, W, I)
% performes the stabilization process in the experimental modal analysis.
% IR : Impulse response matrix (one IR per row)
% DT : Time increment of the impulse response data
% NLIM : Minimum and maximum of system order (2-element vector)
% OMLIM: Circular frequency band (2-element vector)
% I(W) : This indicator function is plotted in the stabilization diagram,
% where W is the circular frequency

% EMA Toolbox
% Peter Fiala, 2008. feb.

%% Diagram initialization
Nmin = Nlim(1); % minimum system order
Nmax = Nlim(2); % maximum system order
% limit the indicator to the selected frequency range
g = find(omvec > omlim(1) & omvec < omlim(2));
omvec = omvec(g);
indi = indi(g);
% scale the indicator to fit into the stabilization diagram
indi = (indi-min(indi))*(Nmax-Nmin)/(max(indi)-min(indi)) + Nmin;
f = figure;
hold on;
plot(omvec/2/pi, indi);
set(gca,'ylim',[Nmin-1 Nmax+1]);
set(gca, 'xlim', omlim/2/pi);
title('Stabilization diagram');
xlabel('Frequency [Hz]');
ylabel('# DOF');
% draw the horozontal lines into the diagram
X = repmat(omlim.'/2/pi,1,Nmax-Nmin+1);
Y = repmat((Nmin:Nmax)-0.5,2,1);
l = line(X,Y);
set(l, 'color', 'black');

%% Stabilization process
lambda1 = lsce(u, dt, Nmin); % initial pole fitting with minium system order
[om1, del1] = lambda2omd(lambda1); % modal frequency and damping
% progbar(Nmin, Nmax, Nmin);
for N = Nmin+1 : Nmax
% progbar(Nmin, Nmax, N);
lambda2 = lsce(u, dt, N); % new pole with higher order
[om2, del2] = lambda2omd(lambda2); % new frequency and damping
stable = zeros(N,2);
for n2 = 1 : length(om2)
% search for similar elements in the old and new om vectors
[omdiff, n1] = min(abs(om1-om2(n2))/om2(n2));
if omdiff < tol % if the minimal distance is smaller than tol
stable(n2,1) = 1; % then the frequency has been stabilized
if abs(del1(n1)-del2(n2))/del2(n2) < tol % same for the damping
stable(n2,2) = 1;
p = plot(om2(n2)/2/pi, N, 'rd'); % omega and delta stable
set(p, 'MarkerFaceColor', [1 0 0]);
else
p = plot(om2(n2)/2/pi, N, 'bd'); % omega stable, delta not
set(p, 'MarkerFaceColor', [0 0 1]);
end
set(p, 'MarkerSize', 6);
set(p, 'ButtonDownFcn', @(src,event)stabilpress(src));
else
p = plot(om2(n2)/2/pi, N, 'o'); % nothing is stable
end
d.lambda = lambda2(2*(n2-1)+(1:2));
d.selected = false;
set(p, 'UserData', d);
end
drawnow;
om1 = om2;
del1 = del2;
end

%% The user picks the modes with the mouse
fprintf(1, 'Select poles in the figure and press Enter when ready\n');
pause;

%% The user picks the modes with the mouse
c = get(get(f, 'CurrentAxes'), 'Children');
d = get(c, 'UserData');
lambda = [];
for n = 1 : length(d)
if ~isempty(d{n})
if d{n}.selected
lambda = [
lambda;
d{n}.lambda;
];
end
end
end
end

%%
function stabilpress(hObject)
d = get(hObject, 'userdata');
if ~d.selected
set(hObject, 'Marker', 'p');
set(hObject, 'MarkerSize', 10);
set(hObject, 'MarkerEdgeColor', 'black');
else
set(hObject, 'Marker', 'd');
set(hObject, 'MarkerSize', 6);
set(hObject, 'MarkerEdgeColor', get(hObject, 'MarkerFaceColor'));
end
d.selected = ~d.selected;
set(hObject, 'UserData', d);
end