MataveID V17.0.0

MataveID is a basic system identification toolbox for both GNU Octave and MATLAB®. MataveID is based on the power of linear algebra and the library is easy to use. MataveID using the classical realization and polynomal theories to identify state space models from data. There are lots of subspace methods in the "old" folder and the reason why I'm not using these files is because they can't handle noise quite well.

I'm building this library because I feel that the commercial libraries are just for theoretical experiments. I'm focusing on real practice and solving real world problems.

Papers:

MataveID contains realization identification, polynomal algorithms and subspace algorithms. They can be quite hard to understand, so I highly recommend to read papers in the "reports" folder about the algorithms if you want to understand how they work, or read the literature.

Literature:

I have been using these books for creating the .m files. All these books have different audience. Some techniques are meant for researchers and some are meant for practical engineering.

Applied System Identification

This book include techniques for linear mechanical systems such as vibrating beams, damping, structural mechanics etc. These techniques comes from NASA and the techniques are created by Jer-Nan Juang. This is a very practical book. The book uses the so called realization theory methods for identify dynamical models from data.

Advantages:

Easy to read and very practical
Include mechanical model buildning
Include impulse, frequency, stochastic, closed loop and recursive identification
These techniques are applied onto Hubble Telescope, Space Shuttle Discovery and Galileo spacecraft

Disadvantages:

Do not include nonlinear system identification and subspace methods
Do not include filtering
MATLAB files from this book is export controlled from NASA = Difficult to download
This book is not produced anymore. I have the PDF.

System Modeling & Identification

This book covering techniques for all types of systems, linear and nonlinear, but it's more a general book for system identfication. Professor Rolf Johansson book contains lots of practice, but also theory as well. More theory and less practice compared to Applied System Identification from Jer-Nan Juang. This book uses both the realization theory methods and subspace methods for identify dynamical systems from data. Also this book includes filters as well such as Uncented Kalman Filter. Can be purchased from https://kfsab.se/sortiment/system-modeling-and-identification/

Advantages:

Easy to read and somtimes practical
Include filtering, statistics and other types of modeling techniques
Include impulse, frequency, stochastic, closed loop, nonlinear and recursive identification
Include both realization theory, subspace and nonlinear system identification methods

Disadvantages:

Do not include closed loop identification
Some methods are difficult to understand how to apply with MATLAB-code. Typical univerity literature for students

Subspace Methods For System Identification

This book include techniques for all types of linear systems. It's a general book of linear system identification. The advantages of this book is that it include modern system identification techniques. The disadvantages about this book is that it contains only theory and no practice, but Professor Tohru Katayama, have made a great work for collecting all these subspace methods. Use this book if you want to have knowledge about the best subspace identification methods.

Advantages:

Include MATLAB code examples and lots of step by step examples
Include stochastic and closed identification
Include the latest methods for linear system identification
Include both realization theory and subspace system identification methods

Disadvantages:

Difficult to read and understand
Does not include impulse, frequency and nonlinear identification
Does not include filtering, statistics and other types of modeling techniques

Adaptive Control

This book is only for adaptive control. But there is one algorithm that are very useful - Recursive Least Squares. This is a very pratical book for applied adaptive control. It's uses the legacy SISO adaptive techniques such as pole placement, Self Tuning Regulator(STR) and Model Reference Adaptive Systems(MRAS) combined with Recursive Least Squares(RLS). If you wonder why only SISO and not MIMO, it's because adaptive control is very difficult to apply in practice and create a reliable controller for all types of systems. The more difficult problem is to solve, the more simplier technique need to be used.

Advantages:

The authors of the book explains which chapters are for pratcial engineering and theoretical researchers
Easy to read
Include both advanced and simple methods depending on which type of problem to solve

Disadvantages:

Only one system identification algorithm is taught
Only SISO model are applied
This book is made for adaptive control and have only one chapter that contains system identification

Examples

Multivariable Output-Error State Space

MOESP is an algorithm that identify a linear state space model. It was invented in 1992. It can both identify SISO and MISO models. Try MOESP or N4SID. They give the same result, but sometimes MOESP can be better than N4SID. It all depends on the data.

[sysd] = mi.moesp(u, y, k, sampleTime, delay, systemorder); % k = Integer tuning parameter such as 10, 20, 25, 32, 47 etc.

Example MOESP

MataveID/examples/moespExample.m

Lines 1 to 22 in 39c9b1b

    
           clc; clear close all; 
        
           [u, t] = mc.gensig('square', 10, 10, 100); 
        
           G = mc.tf(1, [1 0.8 3]); % Model 
        
           y = mc.lsim(G, u, t); % Simulation 
        
           y = y + 0.4*rand(1, length(t)); 
        
           close 
        
           k = 30; 
        
           sampleTime = t(2) - t(1); 
        
           systemorder = 3; 
        
           delay = 0; 
        
           ktune = 0.01; 
        
           [sysd, K] = mi.moesp(u, y, k, sampleTime, ktune, delay, systemorder); % This example works better with MOESP, rather than N4SID 
        
           % Create the observer 
        
           observer = mc.ss(sysd.delay, sysd.A - K*sysd.C, [sysd.B K], sysd.C, [sysd.D sysd.D*0]); 
        
           observer.sampleTime = sysd.sampleTime; 
        
           % Check observer 
        
           [yf, tf] = mc.lsim(observer, [u; y], t); 
        
           close 
        
           plot(tf, yf, t, y) 
        
           grid on

RLS - Recursive Least Squares

RLS is an algorithm that creates a SISO model from data. Here you can select if you want to estimate an ARX, OE model or an ARMAX model, depending on the number of zeros in the polynomal "nze". Select number of error-zeros-polynomal "nze" to 1, and you will get a ARX model or select "nze" equal to model poles "np", you will get an ARMAX model that also includes a kalman gain matrix K. I recommending that. This algorithm can handle data with noise. This algorithm was invented 1821 by Carl Friedrich Gauss, but it was until 1950 when it got its attention in adaptive control.

Use this algorithm if you have data from a open/close loop system and you want to apply that algorithm into embedded system that have low RAM and low flash memory. RLS is very suitable for system that have a lack of memory.

There is a equivalent C-code for RLS algorithm here. Works on ALL embedded systems. https://github.com/DanielMartensson/CControl

[sysd, K] = mi.rls(u, y, np, nz, nze, sampleTime, delay, forgetting);

Notice that there are sevral functions that simplify the use of rls.m

[sysd, K] = mi.oe(u, y, np, nz, sampleTime, delay, forgetting);
[sysd, K] = mi.arx(u, y, np, nz, sampleTime, ktune, delay, forgetting);
[sysd, K] = mi.armax(u, y, np, nz, nze, sampleTime, ktune, delay, forgetting);

Example RLS

This is a hanging load of a hydraulic system. This system is a linear system due to the hydraulic cylinder that lift the load. Here I create two linear first order models. One for up lifting up and one for lowering down the weight. I'm also but a small orifice between the outlet and inlet of the hydraulic cylinder. That's create a more smooth behavior. Notice that this RLS algorithm also computes a Kalman gain matrix.

MataveID/examples/rlsExample.m

Lines 1 to 60 in 2014b74

    
           % Load data 
        
           file = fullfile('..','data','HangingLoad.csv'); 
        
           X = csvread(file); 
        
           t = X(:, 1)'; % Time 
        
           r = X(:, 2)'; % Reference 
        
           y = X(:, 3)'; % Output position 
        
           u = X(:, 4)'; % Input signal from P-controller with gain 3 
        
           sampleTime = 0.02; 
        
           % Do identification of the first data set 
        
           l = length(r) + 2000; % This is half data 
        
           % Poles and zeros 
        
           np = 1; % Number of poles for A(q) 
        
           nz = 1; % Number of zeros of B(q) 
        
           nze = 1; % Number of zeros of C(q) 
        
           % Model up 
        
           u_up = r(1:l/2); 
        
           e_up = randn(1, length(u_up)); % Noise 
        
           y_up = y(1:l/2) + e_up; 
        
           [sysd, K] = mi.rls(u_up, y_up, np, nz, nze, sampleTime); 
        
           % Observer 
        
           sysd_up = mc.ss(0, sysd.A - K*sysd.C, [sysd.B K], sysd.C, [sysd.D 0]); 
        
           sysd_up.sampleTime = sysd.sampleTime; 
        
           % Model down 
        
           u_down = r(l/2+1:end); 
        
           e_down = randn(1, length(u_down)); % Noise 
        
           y_down = y(l/2+1:end) + e_down; 
        
           [sysd, K] = mi.rls(u_down, y_down, np, nz, nze, sampleTime); 
        
           % Observer 
        
           sysd_down = mc.ss(0, sysd.A - K*sysd.C, [sysd.B K], sysd.C, [sysd.D 0]); 
        
           sysd_down.sampleTime = sysd.sampleTime; 
        
           % Simulate model up 
        
           time_up = t(1:l/2); 
        
           [~,~,x] = mc.lsim(sysd_up, [u_up; e_up], time_up); 
        
           hold on 
        
           % Simulate model down 
        
           time_down = t(l/2+1:end); 
        
           x0 = x(:, end); % Initial state 
        
           mc.lsim(sysd_down, [u_down; e_down], time_down, x0); 
        
           % Place legend, title, labels for the signals 
        
           subplot(2, 1, 1) 
        
           legend('Up model', 'Down model', 'Measured'); 
        
           title('Hanging load - Hydraulic system') 
        
           xlabel('Time [s]') 
        
           ylabel('Hanging load position'); 
        
           % Place legend, title, labels for the noise 
        
           subplot(2, 1, 2) 
        
           legend('Noise up', 'Noise down'); 
        
           title('Hanging load - Hydraulic system') 
        
           xlabel('Time [s]') 
        
           ylabel('Noise');

Here we can se that the first model follows the measured position perfect. The "down-curve" should be measured a little bit longer to get a perfect linear model.

SINDy - Sparse Identification of Nonlinear Dynamics

This is a new identification technique made by Eurika Kaiser from University of Washington. It extends the identification methods of grey-box modeling to a much simplier way. This is a very easy to use method, but still powerful because it use least squares with sequentially thresholded least squares procedure. I have made it much simpler because now it also creates the formula for the system. In more practical words, this method identify a nonlinear ordinary differential equations from time domain data.

This is very usefull if you have heavy nonlinear systems such as a hydraulic orifice or a hanging load.

[dx] = mi.sindy(inputs, outputs, degree, lambda, sampleTime);

SINDy Example

This example is a real world example with noise and nonlinearities. Here I set up a hydraulic motor in a test bench and measure it's output and the current to the valve that gives the motor oil. The motor have two nonlinearities - Hysteresis and the input signal is not propotional to the output signal. By using two nonlinear models, we can avoid the hysteresis.

MataveID/examples/sindyExample.m

Lines 1 to 45 in 2014b74

    
           clc; clear; close all; 
        
           % Load CSV data 
        
           file = fullfile('..','data','MotorRotation.csv'); 
        
           X = csvread(file); % Can be found in the folder "data" 
        
           t = X(:, 1); 
        
           u = X(:, 2); 
        
           y = X(:, 3); 
        
           sampleTime = 0.02; 
        
           % Do filtering of y 
        
           y = mi.filtfilt(y', t', 0.1)'; 
        
           % Sindy - Sparce identification Dynamics 
        
           degree = 5; 
        
           lambda = 0.05; 
        
           l = length(u); 
        
           h = floor(l/2); 
        
           s = ceil(l/2); 
        
           fx_up = mi.sindy(u(1:h), y(1:h), degree, lambda, sampleTime); % We go up 
        
           fx_down = mi.sindy(u(s:end), y(s:end), degree, lambda, sampleTime); % We go down 
        
           % Simulation up 
        
           x0 = y(1); 
        
           u_up = u(1:h); 
        
           u_up = u_up(1:100:end)'; 
        
           stepTime = 1.2; 
        
           [x_up, t] = mc.nlsim(fx_up, u_up, x0, stepTime, 'ode15s'); 
        
           % Simulation down 
        
           x0 = y(s); 
        
           u_down = u(s:end); 
        
           u_down = u_down(1:100:end)'; 
        
           stepTime = 1.2; 
        
           [x_down, t] = mc.nlsim(fx_down, u_down, x0, stepTime, 'ode15s'); 
        
           % Compare 
        
           figure 
        
           plot([x_up x_down]) 
        
           hold on 
        
           plot(y(1:100:end)); 
        
           legend('Simulation', 'Measurement') 
        
           ylabel('Rotation') 
        
           xlabel('Time') 
        
           grid on

Square Root Uncented Kalman Filter for parameter estimation

This is Uncented Kalman Filter that using cholesky update method (more stable), instead of cholesky decomposition. This algorithm can estimate parameters to very a complex function if data is available. This method is reqursive and there is a C code version in CControl as well. Use this when you need to estimate parameters to a function if you have data that are generated from that function. It can be for example an object that you have measured data and you know the mathematical formula for that object. Use the measured data with this algorithm and find the parameters for the formula.

[Sw, what] = mi.sr_ukf_parameter_estimation(d, what, Re, x, G, lambda_rls, Sw, alpha, beta, L);

Square Root Uncented Kalman Filter for parameter estimation example

MataveID/examples/sr_ukf_parameter_estimationExample.m

Lines 1 to 62 in 2014b74

    
           clc; clear; close all; 
        
           % Initial parameters 
        
           L = 3;                  % How many states we have 
        
           e = 0.1;                % Tuning factor for noise 
        
           alpha = 0.1;            % Alpha value - A small number like 0.01 -> 1.0 
        
           beta = 2.0;             % Beta value - Normally 2 for gaussian noise 
        
           Re = e*eye(L);          % Initial noise covariance matrix - Recommended to use identity matrix 
        
           Sw = eye(L);            % Initial covariance matrix - Recommended to use identity matrix 
        
           what = zeros(L, 1);     % Estimated parameter vector 
        
           d = zeros(L, 1);        % This is our measurement 
        
           x = [4.4; 6.2; 1.0];    % State vector 
        
           lambda_rls = 1.0;       % RLS forgetting parameter between 0.0 and 1.0, but very close to 1.0 
        
           % Our transition function - This is the orifice equation Q = a*sqrt(P2 - P1) for hydraulics 
        
           G = @(x, w) [w(1)*sqrt(x(2) - x(1)); 
        
                       % We only need to use w(1) so we assume that w(2) and w(3) will become close to 1.0  
        
                        w(2)*x(2); 
        
                        w(3)*x(3)]; 
        
           % Start clock time 
        
           tic 
        
           % Declare arrays  
        
           samples = 100; 
        
           WHAT = zeros(samples, L); 
        
           E = zeros(samples, L); 
        
           % Do SR-UKF for parameter estimation 
        
           for i = 1:samples 
        
             % Assume that this is our measurement  
        
             d(1) = 5 + e*randn(1,1); 
        
             % This is just to make sure w(2) and w(3) becomes close to 1.0 
        
             d(2) = x(2); 
        
             d(3) = x(3); 
        
             % SR-UKF 
        
             [Sw, what] = mi.sr_ukf_parameter_estimation(d, what, Re, x, G, lambda_rls, Sw, alpha, beta, L); 
        
             % Save the estimated parameter  
        
             WHAT(i, :) = what'; 
        
             % Measure the error 
        
             E(i, :) = abs(d - G(x, what))'; 
        
           end 
        
           % Stop the clock 
        
           toc  
        
           % Print the data 
        
           [M, N] = size(WHAT); 
        
           for k = 1:N                                  
        
             subplot(3,1,k); 
        
             plot(1:M, WHAT(:,k), '-', 1:M, E(:, k), '--'); 
        
             title(sprintf('Parameter estimation for parameter w%i', k)); 
        
             ylabel(sprintf('w%i', k)); 
        
             grid on 
        
             legend('Estimated parameters', 'Parameter error') 
        
           end

Square Root Uncented Kalman Filter for state estimation

This is Uncented Kalman Filter that using cholesky update method (more stable), instead of cholesky decomposition. This algorithm can estimate states from a very complex model. This method is reqursive and there is a C code version in CControl as well. Use this when you need to estimate state to a model if you have data that are generated from that function. It can be for example an object that you have measured data and you know the mathematical formula for that object. Use the measured data with this algorithm and find the states for the model.

[S, xhat] = mi.sr_ukf_state_estimation(y, xhat, Rn, Rv, u, F, S, alpha, beta, L);

Square Root Uncented Kalman Filter for state estimation example

MataveID/examples/sr_ukf_state_estimationExample.m

Lines 1 to 68 in 2014b74

    
           clc; clear; close all; 
        
           % Initial parameters 
        
           L = 3;                  % How many states we have 
        
           r = 1.5;                % Tuning factor for noise 
        
           q = 0.2;                % Tuning factor for disturbance 
        
           alpha = 0.1;            % Alpha value - A small number like 0.01 -> 1.0 
        
           beta = 2.0;             % Beta value - Normally 2 for gaussian noise 
        
           Rv = q*eye(L);          % Initial disturbance covariance matrix - Recommended to use identity matrix 
        
           Rn = r*eye(L);          % Initial noise covariance matrix - Recommended to use identity matrix 
        
           S = eye(L);             % Initial covariance matrix - Recommended to use identity matrix 
        
           xhat = [0; 0; 0];       % Estimated state vector 
        
           y = [0; 0; 0];          % This is our measurement 
        
           u = [0; 0; 0];          % u is not used in this example due to the transition function not using an input signal 
        
           x = [0; 0; 0];          % State vector for the system (unknown in reality) 
        
           % Our transition function 
        
           F = @(x, u) [x(2); 
        
                        x(3); 
        
                        0.05*x(1)*(x(2) - x(3))]; 
        
           % Start clock time 
        
           tic 
        
           % Declare arrays  
        
           samples = 200; 
        
           X = zeros(samples, L); 
        
           XHAT = zeros(samples, L); 
        
           Y = zeros(samples, L); 
        
           phase = [90;180;140]; 
        
           amplitude = [1.5;2.5;3.5]; 
        
           % Do SR-UKF for state estimation 
        
           for i = 1:samples 
        
             % Create measurement  
        
             y = x + r*randn(L, 1); 
        
             % Save measurement  
        
             Y(i, :) = y'; 
        
             % Save actual state 
        
             X(i, :) = x'; 
        
             % SR-UKF 
        
             [S, xhat] = mi.sr_ukf_state_estimation(y, xhat, Rn, Rv, u, F, S, alpha, beta, L); 
        
             % Save the estimated parameter  
        
             XHAT(i, :) = xhat'; 
        
             % Update process 
        
             x = F(x, u) + q*amplitude.*sin(i-1 + phase); 
        
           end 
        
           % Stop the clock 
        
           toc  
        
           % Print the data 
        
           [M, N] = size(XHAT); 
        
           for k = 1:N                                  
        
             subplot(3,1,k); 
        
             plot(1:M, Y(:,k), '-g', 1:M, XHAT(:, k), '-r', 1:M, X(:, k), '-b'); 
        
             title(sprintf('State estimation for state x%i', k)); 
        
             ylabel(sprintf('x%i', k)); 
        
             grid on 
        
             legend('y', 'xhat', 'x') 
        
           end

Numerical algorithm for Subspace State Space System IDentification.

N4SID is an algoritm that identify a linear state space model. Use this if you got regular data from a dynamical system. This algorithm can handle both SISO and MISO. N4SID algorithm was invented 1994. If you need a nonlinear state space model, check out the SINDy algorithm. Try N4SID or MOESP. They give the same result, but sometimes N4SID can be better than MOESP. It all depends on the data.

[sysd, K] = mi.n4sid(u, y, k, sampleTime, ktune, delay, systemorder); % k = Integer tuning parameter such as 10, 20, 25, 32, 47 etc. ktune = kalman filter tuning such as 0.1, 0.01 etc

Example N4SID

Here I programmed a Beijer PLC that controls the multivariable cylinder system. It's a nonlinear system, but N4SID can handle it because it's not so nonlinear as a hydraulic motor. Cylinder 0 and Cylinder 1 affecting each other when the propotional control valves opens.

MataveID/examples/n4sidExample.m

Lines 1 to 45 in 4c18a13

    
           clc; clear; close all; 
        
           % Load the data 
        
           file = fullfile('..','data','MultivariableCylinders.csv'); 
        
           X = csvread(file); 
        
           t = X(:, 1); 
        
           r0 = X(:, 2); 
        
           r1 = X(:, 3); 
        
           y0 = X(:, 4); 
        
           y1 = X(:, 5); 
        
           sampleTime = 0.1; 
        
           % Transpose the CSV data 
        
           u = [r0';r1']; 
        
           y = [y0';y1']; 
        
           t = t'; 
        
           % Create the model 
        
           k = 10; 
        
           sampleTime = t(2) - t(1); 
        
           ktune = 0.01; % Kalman filter tuning 
        
           % This won't result well with MOESP and system order = 2 
        
           [sysd, K] = mi.n4sid(u, y, k, sampleTime, ktune); % Delay argment is default 0. Select model order = 2 when n4sid ask you 
        
           % Create the observer 
        
           observer = mc.ss(sysd.delay, sysd.A - K*sysd.C, [sysd.B K], sysd.C, [sysd.D sysd.D*0]); 
        
           observer.sampleTime = sysd.sampleTime; 
        
           % Do simulation 
        
           [outputs, T, x] = mc.lsim(observer, [u; y], t); 
        
           close 
        
           plot(T, outputs(1, :), t, y(1, :)) 
        
           title('Cylinder 0'); 
        
           xlabel('Time'); 
        
           ylabel('Position'); 
        
           grid on 
        
           legend('Identified', 'Measured'); 
        
           ylim([0 12]); 
        
           figure 
        
           plot(T, outputs(2, :), t, y(2, :)) 
        
           title('Cylinder 1'); 
        
           xlabel('Time'); 
        
           ylabel('Position'); 
        
           grid on 
        
           legend('Identified', 'Measured'); 
        
           ylim([0 12]);

OCID - Observer Controller Identification

This is an extention from OKID. The idea is the same, but OCID creates a LQR contol law as well. This algorithm works only for closed loop data. It have its orgin from NASA around 1992 when NASA wanted to identify a observer, model and a LQR control law from closed loop data that comes from an actively controlled aircraft wing in a wind tunnel at NASA Langley Research Center. This algorithm works for both SISO and MIMO models.

Use this algorithm if you want to extract a LQR control law, kalman observer and model from a running dynamical system. Or if your open loop system is unstable and it requries some kind of feedback to stabilize it. Then OCID is the perfect choice.

This OCID algorithm have a particle filter that estimates the markov parameters.

[sysd, K, L] = mi.ocid(r, uf, y, sampleTime, alpha, regularization, systemorder);

OCID Example

MataveID/examples/ocidExample.m

Lines 1 to 59 in 2014b74

    
           clc; close all; clear; 
        
           %% Matrix A 
        
           A = [0 1  0  0; 
        
               -7 -5 0  1; 
        
                0 0  0  1; 
        
                0 1 -8 -5]; 
        
           %% Matrix B 
        
           B = [0 0; 
        
                1 0; 
        
                0 0; 
        
                0 1]; 
        
           %% Matrix C 
        
           C = [1 0 0 0; 
        
                0 0 0 1]; 
        
           %% Model and signals 
        
           delay = 0; 
        
           sys = mc.ss(delay, A, B, C); 
        
           t = linspace(0, 20, 1000); 
        
           r = [linspace(5, -11, 100) linspace(7, 3, 100) linspace(-6, 9, 100) linspace(-7, 1, 100) linspace(2, 0, 100) linspace(6, -9, 100) linspace(4, 1, 100) linspace(0, 0, 100) linspace(10, 17, 100) linspace(-30, 0, 100)]; 
        
           r = [r;2*r]; % MIMO 
        
           %% Feedback 
        
           Q = sys.C'*sys.C; 
        
           R = [1 0; 0 1]; 
        
           L = mc.lqr(sys, Q, R); 
        
           [feedbacksys] = mc.reg(sys, L); 
        
           yf = mc.lsim(feedbacksys, r, t); 
        
           close 
        
           %% Add noise 
        
           v = randn(1, 1000); 
        
           for i = 1:length(yf) 
        
             noiseSigma = 0.10*yf(:, i); 
        
             noise = noiseSigma*v(i); % v = noise, 1000 samples -1 to 1 
        
             yf(:, i) = yf(:, i) + noise; 
        
           end 
        
           %% Identification 
        
           uf = yf(3:4, :); % Input feedback signals 
        
           y = yf(1:2, :); % Output feedback signals 
        
           regularization = 10000; 
        
           modelorder = 6; 
        
           sampleTime = t(2) - t(1); 
        
           alpha = 20; % Filtering integer parameter 
        
           [sysd, K, L] = mi.ocid(r, uf, y, sampleTime, alpha, regularization, modelorder); 
        
           %% Validation 
        
           u = -uf + r; % Input signal %u = -Lx + r = -uf + r 
        
           yt = mc.lsim(sysd, u, t); 
        
           close 
        
           %% Check 
        
           plot(t, yt(1:2, 1:2:end), t, yf(1:2, :)) 
        
           legend("Identified 1", "Identified 2", "Data 1", "Data 2", 'location', 'northwest') 
        
           grid on

Oriented FAST Rotated Pattern

Use this algorithm if you want to convert keypoints into one large binary matrix for image classification and detection.

[data, X1, X2, G, corners, scores] = mi.orp(X, sigma1, sigma2, threshold_sobel, threshold_fast, fast_method);

Example

MataveID/examples/orpExample.m

Lines 1 to 37 in a300c9e

    
           % Clear 
        
           close all 
        
           clear all 
        
           clc 
        
           % Load image 
        
           X = imread(fullfile('..', 'data', 'hus.jpg')); 
        
           % Make image greyscale 
        
           if size(X, 3) > 1 
        
             X =  rgb2gray(X); 
        
           end 
        
           % Compute Oriented FAST Rotated Pattern 
        
           sigma1 = 1;             % Backgroun filtering 
        
           sigma2 = 6;             % Filtering for the descriptors 
        
           threshold_sobel = 127;  % Threshold for the corners 
        
           threshold_fast = 50;    % Threshold for the corners 
        
           fast_method = 9;        % FAST method: 9, 10, 11, 12 
        
           [data, X1, X2, G, corners] = mi.orp(X, sigma1, sigma2, threshold_sobel, threshold_fast, fast_method); 
        
           % Plot 
        
           figure 
        
           imshow(uint8(X)) 
        
           hold on 
        
           plot(corners(:, 1), corners(:, 2), 'r.'); 
        
           title('Corner detection') 
        
           hold off 
        
           figure 
        
           imshow(uint8(X1)) 
        
           title('Gaussian filter - Background') 
        
           figure 
        
           imshow(uint8(X2)) 
        
           title('Gaussian filter - For descriptor') 
        
           figure 
        
           imshow(uint8(G)); 
        
           title('Sobel filter - For corner/edge detection')

Results

Sobel filter

Use this filter if you want to find the gradients and the orientations inside an image

[G, O] = mi.sobel(image);

Example

MataveID/examples/sobelExample.m

Lines 1 to 23 in 1fb0316

    
           % Close all 
        
           close all 
        
           % Read image 
        
           image = imread(fullfile('..', 'data', 'happy.gif')); 
        
           % Turn to grey scale 
        
           if(size(image, 3) == 3) 
        
             image = rgb2gray(image); 
        
           end 
        
           % Compute gradients and orientations 
        
           [G, O] = mi.sobel(image); 
        
           % Show original image 
        
           subplot(1, 2, 1); 
        
           imshow(image); 
        
           title('Original image'); 
        
           % Show gradient image 
        
           subplot(1, 2, 2); 
        
           imshow(uint8(G)); 
        
           title('Gradient image');

Result

Particle Filter - Nonlinear filter

A particle filter is another estimation filter such as Square Root Uncented Kalman Filter (SR-UKF), but SR-UKF assume that the noise is gaussian (normally distributed) and SR-UKF requries a dynamical model. The particle filter does not require the user to specify a dynamical model and the particle filter assume that the noise can be non-gaussian or gaussian, nonlinear in other words.

The particle filter is using Kernel Density Estimation algorithm to create the internal probability model, hence the user only need to specify one parameter with the following example. If you don't have a model that describes the dynamical behaviour, this filter is the right choice for you then.

[xhat, horizon, k, noise] = mi.pf(x, xhatp, k, horizon, noise);

Particle Filter example 1

MataveID/examples/particleFilterExample1.m

Lines 1 to 44 in 817443d

    
           clc; clear; close all; 
        
           % Create inputs 
        
           N = 200; 
        
           u = linspace(1, 1, N); 
        
           u = [5*u 10*u -4*u 3*u 5*u 0*u -5*u 0*u]; 
        
           % Create time 
        
           t = linspace(0, 100, length(u)); 
        
           % Create second order model 
        
           G = mc.tf(1, [1 0.8 3]); 
        
           % Simulate outputs 
        
           y = mc.lsim(G, u, t); 
        
           close 
        
           % Add noise 
        
           e = 0.1*randn(1, length(u)); 
        
           y = y + e; 
        
           % Do particle filtering - Tuning parameters 
        
           p = 4;                            % Length of the horizon (Change this) 
        
           % Particle filter - No tuning 
        
           [m, n] = size(y);                  % Dimension of the output state and length n 
        
           yf = zeros(m, n);                  % Filtered outputs 
        
           horizon = zeros(m, p);             % Horizon matrix 
        
           xhatp = zeros(m, 1);               % Past estimated state 
        
           xhatp(1) = y(1);                   % First estimation is a real measurement 
        
           k = 1;                             % Horizon counting (will be counted to p. Do not change this) 
        
           noise = rand(m, p);                % Random noise, not normal distributed 
        
           % Particle filter - Simulation 
        
           for i = 1:n 
        
             x = y(:, i);                     % Get the state 
        
             [xhat, horizon, k, noise] = mi.pf(x, xhatp, k, horizon, noise); 
        
             yf(:, i) = xhat;                 % Estimated state 
        
             xhatp = xhat;                    % This is the past estimated state 
        
           end 
        
           % Plot restult 
        
           plot(t, y, t, yf, '-r') 
        
           grid on

Particle Filter example 2

MataveID/examples/particleFilterExample2.m

Lines 1 to 32 in 817443d

    
           clc; clear; close all; 
        
           file = fullfile('..','data','ParticleFilterDataRaw.csv'); 
        
           X = dlmread(file,';',1,0); 
        
           t = X(:, 1)'; 
        
           y = X(:, 2)'; 
        
           % Do particle filtering - Tuning parameters 
        
           p = 14;                            % Length of the horizon (Change this) 
        
           % Particle filter - No tuning 
        
           [m, n] = size(y);                  % Dimension of the output state and length n 
        
           yf = zeros(m, n);                  % Filtered outputs 
        
           horizon = zeros(m, p);             % Horizon matrix 
        
           xhatp = zeros(m, 1);               % Past estimated state 
        
           xhatp(1) = y(1);                   % First estimation is a real measurement 
        
           k = 1;                             % Horizon counting (will be counted to p. Do not change this) 
        
           noise = rand(m, p);                % Random noise, not normal distributed 
        
           % Particle filter - Simulation 
        
           for i = 1:n 
        
             x = y(:, i);                     % Get the state 
        
             [xhat, horizon, k, noise] = mi.pf(x, xhatp, k, horizon, noise); 
        
             yf(:, i) = xhat;                 % Estimated state 
        
             xhatp = xhat;                    % This is the past estimated state 
        
           end 
        
           % Plot restult 
        
           plot(t, y) 
        
           hold on 
        
           plot(t, yf, '-r') 
        
           grid on

BJ - Box-Jenkins

Box-Jenkins is a special case when a system model sysd and a disturbance model sysh need to be found. The disturbance is difficult to know and with this Box-Jenkins algorithm, then the user can identify the disturbance model and create an observer of it by using the kalman gain matrices K1, K2. Notice that this Box-Jenkins algorithm using subspace methods, instead of classical polynomial methods.

The disturbance model can be used for:

Create a disturbance simulation with feedback control
Create filtering for sensors

[sysd, K1, sysh, K2] = mi.bj(u, y, k, sampleTime, ktune, delay, systemorder_sysd, systemorder_sysh);

Example

MataveID/examples/bjExample.m

Lines 1 to 50 in 328f7f5

    
           % Clear all 
        
           clear all 
        
           close all 
        
           % Create system model 
        
           G = mc.tf(1, [1 1.5 2]); 
        
           % Create disturbance model 
        
           H = mc.tf([2 3], [1 5 6]); 
        
           % Create input signal 
        
           [u, t] = mc.gensig('square', 10, 10, 100); 
        
           u = [u*5 u*2 -u 10*u -2*u]; 
        
           t = linspace(0, 30, length(u)); 
        
           % Create disturbance signal 
        
           e = randn(1, length(t)); 
        
           % Simulate with noise 
        
           d = mc.lsim(H, e, t); 
        
           y = mc.lsim(G, u, t) + d 
        
           close 
        
           % Use Box-Jenkins to find the system model and the disturbance model 
        
           k = 50; 
        
           sampleTime = t(2) - t(1); 
        
           ktune = 0.5; 
        
           delay = 0; 
        
           systemorder_sysd = 2; 
        
           systemorder_sysh = 2; 
        
           [sysd, K1, sysh, K2] = mi.bj(u, y, k, sampleTime, ktune, delay, systemorder_sysd, systemorder_sysh); 
        
           % Plot sysd 
        
           [sysd_y, sysd_t] = mc.lsim(sysd, u, t); 
        
           close all 
        
           plot(t, y, sysd_t, sysd_y); 
        
           legend('Measurement', 'Identified') 
        
           grid on 
        
           title('System model', 'FontSize', 20) 
        
           % Plot sysh 
        
           figure(2) 
        
           [sysh_y, sysh_t] = mc.lsim(sysh, e, t); 
        
           close(2) 
        
           figure(2) 
        
           plot(t, d, sysh_t, sysh_y); 
        
           legend('Measurement', 'Identified') 
        
           grid on 
        
           title('Disturbance model', 'FontSize', 20)

Results

Canny filter

Use this filter if you want to find the edges inside an image

[E] = mi.canny(image);

Example

MataveID/examples/cannyExample.m

Lines 1 to 23 in 210c93d

    
           % Close all 
        
           close all 
        
           % Read image 
        
           image = imread(fullfile('..', 'data', 'way.jpg')); 
        
           % Turn to grey scale 
        
           if(size(image, 3) == 3) 
        
             image = rgb2gray(image); 
        
           end 
        
           % Compute gradients and orientations 
        
           E = mi.canny(image); 
        
           % Show original image 
        
           subplot(1, 2, 1); 
        
           imshow(image); 
        
           title('Original image'); 
        
           % Show gradient image 
        
           subplot(1, 2, 2); 
        
           imshow(uint8(E)); 
        
           title('Edge image');

Result

Notice that Canny is quite slow, but gives very thin edges, which is good. But if you only want to have the edges and you don't care how thick they are. Then Sobel is the right solution for you because Sobel is much faster than Canny.

>> G = mi.sobel(imread('way.jpg'));
>> G(G < 255) = 0; % Every pixel that are not white is going to be black
>> imshow(uint8(G));

Result

Principal Component Analysis

Principal Component Analysis can be used for dimension reduction and projection on maximum variance between classes.

[P, W, mu] = mi.pca(X, c);

Principal Component Analysis example

MataveID/examples/pcaExample.m

Lines 1 to 60 in 48b6eaf

    
           % First clear and close all figures 
        
           clear 
        
           close all 
        
           % Create data 
        
           X = [-3244 5324 1345; 
        
               10 30 100; 
        
               20 40 93; 
        
               30 60 163; 
        
               60 100 126; 
        
               55 13 134; 
        
               306 34 104; 
        
               316 56 120; 
        
               326 74 127; 
        
               337 80 128; 
        
               347 89 131; 
        
               358 103 139; 
        
               -31 -56 -120; 
        
               -32 -74 -127; 
        
               -33 -80 -128; 
        
               -34 -89 -131; 
        
               -35 -103 -139; 
        
               700 600 500; 
        
               1000 1000 1000; 
        
               -3231 4345 -4352; 
        
               -2342 4356 3453; 
        
               -2364 4326 3353; 
        
               658 143 1692]; 
        
           % Noisy data 
        
           X = [X; 10*randn(10, 3)]; 
        
           % Remove noise 
        
           y = mi.dbscan(X, 50, 2); 
        
           X = X(y ~= 0, :); 
        
           % Plot original data 
        
           cmap = jet(length(X)); 
        
           scatter3(X(:, 1), X(:, 2), X(:, 3), 50,cmap) 
        
           title('Original 3D data', 'FontSize', 20) 
        
           % Do PCA 
        
           c = 2; 
        
           [P, W] = mi.pca(X, c); 
        
           % Reconstruct 
        
           X_reconstructed = W * P; 
        
           reconstruction_error = norm(X - X_reconstructed, 'fro') / norm(X, 'fro'); 
        
           figure 
        
           switch(c) 
        
           case 3 
        
             scatter3(X_reconstructed(:, 1), X_reconstructed(:, 2), X_reconstructed(:, 3)); 
        
           case 2 
        
             scatter(X_reconstructed(:, 1), X_reconstructed(:, 2)); 
        
           case 1 
        
             scatter(X_reconstructed(:, 1), 0*X_reconstructed(:, 1)); 
        
           end 
        
           grid on 
        
           title(sprintf('Reconstruction %iD with error %f', c, reconstruction_error), 'FontSize', 20)

CCA - Canonical Correlation Analysis

If N4SID won't work for you due to high noise measurement, then CCA is an alternative method to use. CCA returns a state space model and a kalman gain matrix K.

[sysd, K] = mi.cca(u, y, k, sampleTime, delay); % k = Integer tuning parameter such as 10, 20, 25, 32, 47 etc.

Example CCA

MataveID/examples/ccaExample.m

Lines 1 to 38 in 2014b74

    
           clc; clear; close all; 
        
           % Create model 
        
           G = mc.tf(1, [1 1.5 1]); 
        
           % Create control inputs 
        
           [u, t] = mc.gensig('square', 10, 10, 100); 
        
           u = [u*5 u*2 -u 10*u -2*u]; 
        
           t = linspace(0, 50, length(u)); 
        
           % Simulate 
        
           y = mc.lsim(G, u, t); 
        
           close 
        
           % Add noise 
        
           yn = y + randn(1, length(y)); 
        
           % Identify the model 
        
           sampleTime = t(2) - t(1); 
        
           delay = 0; 
        
           systemOrder = 2; 
        
           k = 30; 
        
           [sysd, K] = mi.cca(u, yn, k, sampleTime, delay, systemOrder); 
        
           % Create an observer 
        
           delay = sysd.delay; 
        
           A = sysd.A; 
        
           B = sysd.B; 
        
           C = sysd.C; 
        
           D = sysd.D; 
        
           observer = mc.ss(delay, A - K*C, [B K], C, [D 0]); 
        
           observer.sampleTime = sysd.sampleTime; 
        
           % Simulate the observer 
        
           [yobs, tobs] = mc.lsim(observer, [u; yn], t); 
        
           close 
        
           plot(t, yn, '-r', tobs, yobs, '-b'); 
        
           grid on

Density-based spatial clustering of applications with noise

This is a cluster algorithm that can identify the amount of clusters. This algorithm requries two tuning parameters, epsilon and min_pts, which stands for radius and minimum points. This algorithm does not work if you have varying densities, else this algorithm is considered to be one of the best clustering algorithms. So, make sure that all your classes have the same amount of variance before you are using this algorithm due to its robustness against noise/outliers.

It exist an equivalent C-code dbscan inside CControl repository.

[idx] = mi.dbscan(X, epsilon, min_pts);

Example Density-based spatial clustering of applications with noise

MataveID/examples/dbscanExample.m

Lines 1 to 110 in 7e1b828

    
           clc; 
        
           clear; 
        
           close all; 
        
           % Clusters 
        
           X = [-0.251521,   1.045117,  -1.281658, 
        
                -1.974109,   0.278170,  -1.023392, 
        
                -0.957729,  -0.977450,   0.477872, 
        
                -0.449159,  -1.016680,   0.095610, 
        
                -1.785787,  -1.403543,   0.483454, 
        
                 1.366889,  -0.762590,  -1.162454, 
        
                 2.129839,   0.358568,  -2.118250, 
        
                 0.751071,  -1.766582,   0.178434, 
        
                -1.980847,  -1.320933,  -0.457778, 
        
                -0.478030,   0.606917,  -1.630624, 
        
                 3.674916,   0.088957,   0.877373, 
        
                 0.637213,   0.079176,   0.342038, 
        
                 1.142329,   0.629997,   0.311134, 
        
                -0.878974,   0.042527,   0.736522, 
        
                 1.751637,  -1.434299,  -1.325140, 
        
                 1.110682,   1.091970,   1.434869, 
        
                -0.504482,  -2.504821,  -1.245315, 
        
                -0.102915,  -0.203266,  -0.849767, 
        
                -0.822834,   1.158801,  -0.405579, 
        
                -1.278287,   0.391306,   0.857077, 
        
                10.6772,   10.7365,    9.9264, 
        
                 8.7785,   11.1680,    9.5915, 
        
                 8.6872,    9.6464,   10.3801, 
        
                10.0142,    8.8311,    9.2021, 
        
                 8.4179,    9.8572,   11.6356, 
        
                 9.8487,   10.4979,   10.8620, 
        
                10.0957,    9.7878,   12.2653, 
        
                11.4528,   11.5186,   10.3050, 
        
                10.9284,    9.9654,   10.4562, 
        
                 8.5272,   10.7451,    9.8355, 
        
                10.1508,   10.2318,   10.2417, 
        
                10.7342,   10.0689,    9.9918, 
        
                10.4784,    9.2032,   10.6060, 
        
                10.1309,    9.4392,   10.9674, 
        
                10.6971,   10.3347,   11.0447, 
        
                 7.9489,    9.4566,    9.5258, 
        
                10.4827,   10.3030,   10.5582, 
        
                10.4496,   10.3880,   11.1661, 
        
                11.0291,   10.0233,    9.9280, 
        
                 9.0638,    9.3650,    9.3670 
        
                -34.233,  -30.841,  -31.720, 
        
                -32.629,  -31.786,  -31.290, 
        
                -31.466,  -31.984,  -33.254, 
        
                -31.878,  -33.052,  -31.761, 
        
                -33.528,  -30.921,  -32.836, 
        
                -31.793,  -32.082,  -30.453, 
        
                -31.812,  -32.417,  -31.874, 
        
                -32.127,  -32.599,  -32.806, 
        
                -32.979,  -32.096,  -31.754, 
        
                -31.759,  -31.925,  -31.313, 
        
                -30.531,  -31.838,  -31.179, 
        
                -32.168,  -31.928,  -30.649, 
        
                -31.049,  -32.092,  -31.408, 
        
                -33.006,  -31.753,  -31.961, 
        
                -32.092,  -32.391,  -31.501, 
        
                -31.184,  -31.634,  -32.802, 
        
                -30.658,  -31.616,  -31.493, 
        
                -31.958,  -31.694,  -31.425, 
        
                -33.114,  -32.029,  -31.459, 
        
                -31.081,  -34.486,  -32.020, 
        
                % Outliers 
        
                45, 43, 0, 
        
                23, -3, 2, 
        
                32, 54, 23, 
        
                23, 51, 77, 
        
                -22, -34, 53 
        
                ]; 
        
           % Run dbscan clustering algorithm 
        
           epsilon = 3; 
        
           min_pts = 3; 
        
           [index] = mi.dbscan(X, epsilon, min_pts) 
        
           % Get number of clusters 
        
           clusters = max(index); 
        
           legends = {}; 
        
           for i = 0:clusters 
        
             % Get the clusters 
        
             Xi = X(find(index == i), :); 
        
             % Create legends 
        
             if i~=0 
        
               legends{end+1} = ['Cluster #' num2str(i)]; 
        
              else 
        
               if ~isempty(Xi) 
        
                 legends{end+1} = 'Noise'; 
        
               end 
        
             end 
        
             % Plot 
        
             if ~isempty(Xi) 
        
               scatter3(Xi(:,1),Xi(:,2), Xi(:, 3)); 
        
             end 
        
             hold on; 
        
           end 
        
           grid on; 
        
           legend(legends); 
        
           legend('Location', 'NorthEastOutside');

MataveID/examples/dbscanNonGaussianExample.m

Lines 1 to 1035 in 837fe5a

    
           clc; 
        
           clear; 
        
           close all; 
        
           X = [-2.2594e-02  -2.8426e-01 
        
             -6.1004e-02  -1.2529e+00 
        
              1.8012e+00  -4.3953e-02 
        
             -2.3199e-01  -4.3215e-01 
        
              4.5101e-01   5.0945e-01 
        
             -3.3231e-01   3.5719e-01 
        
              1.2116e+00  -9.8216e-01 
        
             -3.1336e-01   1.1489e+00 
        
             -9.2645e-01   4.3461e-01 
        
              1.0152e+00  -1.0830e+00 
        
             -3.2503e-01   2.4147e-01 
        
             -3.3034e-02  -3.4398e-01 
        
             -6.7431e-01   7.8376e-01 
        
             -3.2854e-01  -5.3393e-01 
        
             -9.9127e-01  -1.2883e+00 
        
             -4.4060e-01   9.6685e-01 
        
             -3.8547e-01  -8.6362e-01 
        
              4.6876e-01  -1.3097e-01 
        
             -9.4030e-02  -3.6845e-01 
        
              2.7813e-01   1.6698e-02 
        
              7.5687e-01  -2.0824e+00 
        
             -1.0325e+00  -1.4041e+00 
        
             -3.3536e-02   4.5628e-01 
        
              9.0834e-01   1.6325e-01 
        
             -1.1218e+00   2.1863e-02 
        
              4.0391e-01  -1.2450e+00 
        
              1.6359e+00   5.2046e-01 
        
              8.0423e-01  -1.5942e+00 
        
             -7.9504e-02  -5.6721e-01 
        
             -4.2685e-01   4.9866e-02 
        
              1.6710e-01   1.4084e-01 
        
              2.3298e-01  -1.1611e-01 
        
             -6.2643e-01   9.0413e-02 
        
              1.3640e+00   8.3066e-01 
        
              6.5098e-01  -1.0465e+00 
        
             -1.1512e-01  -7.7194e-01 
        
              4.0594e-04   7.4381e-01 
        
             -1.5233e+00   1.4047e-01 
        
              1.3478e+00  -7.9848e-01 
        
             -8.1220e-01   5.7723e-02 
        
              1.3301e+00  -1.2374e+00 
        
             -1.3062e+00   6.4647e-01 
        
              5.0703e-01   4.4456e-01 
        
             -3.2293e-02  -5.5105e-01 
        
              5.4644e-01   2.6747e-01 
        
             -1.1859e-01   1.2764e+00 
        
             -6.2604e-01  -9.5314e-02 
        
              4.5376e-01   2.8721e-01 
        
              9.6966e-01   5.3815e-01 
        
              1.2529e+00   7.2463e-01 
        
             -5.9314e-01   1.3138e+00 
        
             -1.4222e+00   8.9346e-01 
        
             -9.2827e-02   1.4011e+00 
        
              1.4324e+00   5.8978e-01 
        
              4.4863e-01   3.6086e-01 
        
             -1.3215e-02  -3.3732e-01 
        
              1.6081e-01   5.4614e-01 
        
             -5.8436e-01  -1.5588e+00 
        
             -6.1729e-01   6.9923e-01 
        
              5.1493e-01  -1.4936e+00 
        
              1.3864e-02  -6.4726e-01 
        
              1.9678e+00  -5.7283e-02 
        
              5.1616e-01  -7.9715e-01 
        
              1.2910e+00  -1.4818e+00 
        
              1.0457e+00   6.2221e-01 
        
             -5.4134e-01   1.0099e+00 
        
              2.6611e-01  -1.0964e+00 
        
             -3.2171e-01   6.3858e-01 
        
              1.4057e+00   1.9043e+00 
        
             -7.8970e-01   1.0897e+00 
        
              1.3863e-01  -8.1524e-01 
        
             -8.9715e-01   7.1834e-01 
        
             -1.0474e+00   1.0906e+00 
        
             -1.2458e-01   3.1919e-01 
        
              5.0943e-01  -1.7801e-01 
        
             -2.4136e-01  -1.5820e+00 
        
             -3.9367e-01  -7.6766e-01 
        
              1.4186e+00  -1.7849e-02 
        
              6.7347e-01  -1.4257e+00 
        
             -8.9084e-01   6.3767e-01 
        
              4.2386e-01   4.7502e-01 
        
             -2.5999e-01  -1.0623e-02 
        
             -1.3009e-01   1.9707e-01 
        
              1.1956e+00   3.6417e-01 
        
             -1.9497e+00  -9.5261e-01 
        
              2.9175e-01   6.6693e-01 
        
              1.5844e+00   1.1732e+00 
        
             -6.4014e-01  -7.3549e-01 
        
              6.0286e-01   6.2266e-01 
        
             -5.0203e-01  -2.3024e-01 
        
              7.0439e-01  -2.3294e-01 
        
             -1.2900e+00  -3.0045e-01 
        
             -1.1341e+00  -4.9062e-01 
        
             -7.0234e-01  -8.5860e-02 
        
              8.7111e-01  -1.0650e+00 
        
              8.9724e-01  -4.3898e-01 
        
             -5.2737e-01   3.9492e-02 
        
             -1.1465e+00   1.2651e+00 
        
              5.6792e-01  -3.5536e-01 
        
              1.1580e+00  -7.7697e-01 
        
             -9.3490e-01  -5.2147e-01 
        
             -1.2923e-01  -1.3970e+00 
        
              5.0877e-01  -5.9896e-01 
        
             -1.9880e-01  -6.0254e-01 
        
             -2.4738e-01   3.0433e-01 
        
              8.7637e-01   5.3424e-01 
        
              4.3255e-01  -2.4342e-01 
        
             -1.3875e+00   1.2751e+00 
        
             -4.4636e-01  -2.2641e+00 
        
              1.7463e+00  -1.3649e+00 
        
             -5.1376e-01   2.3474e-01 
        
              7.9510e-01   7.1752e-01 
        
              2.6842e-01   6.2655e-01 
        
              1.8058e-01   1.2069e+00 
        
             -7.0412e-01   6.6662e-01 
        
              7.4482e-01  -7.1173e-01 
        
             -5.0069e-01  -1.3617e+00 
        
              1.9838e+00  -2.5069e-02 
        
             -6.3467e-01   3.1359e-01 
        
              9.4120e-02  -1.2647e-01 
        
             -4.2874e-01  -1.6351e+00 
        
              1.2893e-01   1.0993e-01 
        
              3.8425e-01   1.1330e+00 
        
             -1.1372e+00  -1.7900e+00 
        
              1.7557e-02  -6.4322e-01 
        
              2.1162e-01   9.8959e-01 
        
             -7.4307e-01  -1.4591e+00 
        
              8.3507e-01   4.6402e-01 
        
              4.7851e-01  -4.3382e-01 
        
              8.1312e-01   2.1421e-01 
        
              4.7624e-01   1.4632e+00 
        
             -9.5593e-02   6.2127e-01 
        
             -1.2288e+00  -1.1578e-01 
        
             -1.5070e-01   3.3992e-01 
        
             -9.8036e-01  -2.9709e-01 
        
              1.6111e-01  -5.2987e-02 
        
             -4.0171e-01   6.7162e-02 
        
             -2.9232e-02  -4.1439e-01 
        
              4.3478e-02  -4.5269e-02 
        
              1.7912e+00   8.1039e-01 
        
             -5.0055e-01   5.1394e-01 
        
              8.8572e-01  -3.8679e-01 
        
             -2.8815e-01  -1.0590e+00 
        
             -1.0037e+00  -3.9236e-01 
        
              9.2463e-01  -8.0038e-01 
        
             -1.3819e+00  -6.8772e-01 
        
              4.0937e-01   1.2109e+00 
        
             -8.1253e-01   6.1003e-02 
        
              2.0177e-01  -1.6219e+00 
        
              1.8623e+00  -1.4742e+00 
        
              1.6828e+00   6.4318e-01 
        
              4.8329e-01   1.8501e-01 
        
              4.6985e-01   4.3149e-01 
        
             -2.3324e-01   1.3160e+00 
        
             -3.0353e-01  -4.0624e-01 
        
              6.9906e-01  -2.0681e-01 
        
             -7.5593e-01  -7.1684e-01 
        
             -1.3730e-01  -1.0454e+00 
        
              1.8132e+00  -7.9082e-01 
        
              6.6350e-02   3.4364e-01 
        
             -3.5959e-01   7.4675e-01 
        
             -1.6664e+00  -4.2597e-01 
        
             -1.9181e-01  -1.1711e+00 
        
             -3.1683e-01   5.6108e-01 
        
             -3.3534e-01  -1.6527e+00 
        
             -3.7552e-01  -5.4841e-01 
        
              1.8554e-02  -1.1342e+00 
        
              1.0121e+00  -6.1726e-01 
        
             -8.7643e-01  -2.0309e-01 
        
              1.0830e-01   1.3623e+00 
        
              2.1879e-03   6.9546e-01 
        
             -9.3998e-01  -3.5847e-01 
        
              1.3028e+00   2.6779e-03 
        
              1.1960e+00  -1.1866e+00 
        
              5.9888e-01  -1.6982e-01 
        
             -2.2740e-01   1.0280e-01 
        
             -8.7580e-01  -5.9513e-01 
        
              4.4730e-01  -2.8122e-01 
        
              5.1289e-01   5.7525e-01 
        
              4.4592e-02  -1.1267e+00 
        
             -1.1201e+00   9.1469e-01 
        
              1.3613e+00   1.4609e+00 
        
             -1.1503e+00  -5.0065e-01 
        
              3.0240e-01   1.9900e-01 
        
             -2.2838e-01  -1.8653e-02 
        
             -2.2641e+00  -1.5467e-01 
        
              9.8325e-02  -8.9050e-01 
        
             -7.4354e-01  -3.5522e-01 
        
              7.7523e-01   5.2970e-01 
        
              5.2140e-01   6.0336e-01 
        
             -1.3399e+00  -7.5159e-01 
        
             -3.5416e-01   1.0219e-01 
        
              1.3219e-01   4.3036e-01 
        
              1.6526e-01   1.0786e+00 
        
              1.1448e+00   4.0254e-01 
        
             -7.8574e-01  -1.6350e+00 
        
              3.0571e-02  -5.3033e-01 
        
              1.1864e-01  -5.7182e-01 
        
             -7.0881e-01  -3.5339e-01 
        
              9.0962e-01  -3.6339e-01 
        
             -1.1097e+00   8.1711e-01 
        
              4.6581e-01   6.5476e-02 
        
             -6.9856e-01   1.5500e+00 
        
              3.9624e-01   3.3903e-01 
        
             -7.6366e-01  -1.9503e-01 
        
              2.4550e-01   6.4408e-01 
        
             -1.7021e+00  -4.7724e-01 
        
             -1.1912e+00  -4.1702e-01 
        
             -5.2656e-01  -1.1103e-01 
        
              1.3801e-01  -1.0644e+00 
        
              1.1659e-01  -3.4802e-01 
        
             -3.2305e-02  -7.7581e-01 
        
             -7.0772e-01  -6.4922e-01 
        
             -3.8910e-01  -4.2625e-01 
        
             -1.1980e+00  -1.2898e+00 
        
             -1.0149e+00  -4.1357e-01 
        
             -2.1268e+00   3.6212e-03 
        
              1.4508e-01   1.2375e+00 
        
              1.1223e-01   1.1133e+00 
        
              1.8048e-01  -2.0847e-01 
        
             -9.6554e-01  -1.6057e-01 
        
              2.5574e-01   9.6502e-01 
        
              2.8972e-01   2.0669e-01 
        
              1.4599e-01  -9.2372e-01 
        
              1.9805e+00  -3.0612e-01 
        
             -7.2656e-01  -2.6484e-01 
        
             -1.5556e-01  -1.1237e+00 
        
             -1.0872e+00  -2.1061e-01 
        
              9.8903e-01   1.2864e+00 
        
              1.3381e+00   6.4767e-01 
        
             -1.2577e+00   9.7166e-01 
        
             -1.0540e-01   6.9013e-02 
        
              1.5760e+00   1.4252e-02 
        
              1.4796e-01   8.9120e-01 
        
              4.4674e-02  -2.7585e-01 
        
              4.9260e-01   2.5560e-01 
        
              5.4535e-01  -6.3176e-01 
        
             -8.8599e-01  -1.5476e+00 
        
              2.0921e-01  -4.4986e-01 
        
             -4.2817e-01   1.4654e-01 
        
             -1.1914e+00   7.9870e-01 
        
             -9.4725e-01   1.0820e-01 
        
              3.1469e-02  -1.0509e+00 
        
             -2.3799e-01   6.3451e-01 
        
             -9.0205e-01   8.4549e-03 
        
              4.3721e-01  -5.2694e-01 
        
             -6.5788e-01  -5.1944e-01 
        
              9.3045e-01  -9.8878e-01 
        
              7.0144e-01   2.2506e-01 
        
             -2.8446e-01   5.1133e-01 
        
              7.2686e-01   2.4091e-01 
        
              1.8543e-01   5.5498e-01 
        
             -1.3096e-01   8.8532e-02 
        
              3.4370e-01   6.3674e-02 
        
             -2.1059e-01  -2.8671e-01 
        
             -1.0803e+00   1.3199e+00 
        
              6.9216e-01  -8.9556e-01 
        
              2.8898e-01   3.8980e-01 
        
             -5.1632e-01   1.0011e+00 
        
              2.2772e-01   6.5038e-01 
        
              3.0817e-02   8.0929e-01 
        
             -1.2333e+00  -4.6792e-02 
        
             -1.7284e+00  -4.5820e-01 
        
             -1.6015e-01  -3.1036e-02 
        
             -1.9803e+00   9.2639e-01 
        
              1.6141e-01   5.6222e-01 
        
             -1.1604e+00  -4.2928e-01 
        
              2.4349e-01  -1.7070e+00 
        
             -1.0081e+00   1.1214e-01 
        
             -9.8005e-01   1.5151e+00 
        
              1.0284e+00   6.1577e-01 
        
             -2.7440e-02   1.4227e+00 
        
             -8.9793e-01   6.1860e-01 
        
              1.4029e-01  -1.2717e+00 
        
             -3.1327e-01  -3.5770e-01 
        
             -3.6321e-01   2.3765e-01 
        
             -2.5930e-01  -6.6566e-01 
        
             -1.7959e+00   4.5196e-02 
        
             -1.3197e+00  -1.0249e+00 
        
             -1.0818e+00   3.4267e-01 
        
              4.0395e-01   1.0319e+00 
        
             -7.7854e-02  -9.0068e-01 
        
              2.5850e-01   2.8964e-01 
        
             -6.5539e-01   4.6649e-01 
        
             -1.3545e+00  -4.1875e-02 
        
              3.1113e-01  -1.8785e-01 
        
              9.0605e-01  -6.1710e-01 
        
             -1.2463e-01  -1.3233e-01 
        
              9.2519e-01  -5.7300e-01 
        
              1.4881e-01   2.9275e-01 
        
              9.8515e-01  -1.2500e+00 
        
             -5.8100e-01  -1.4843e+00 
        
             -5.0102e-01   1.6918e+00 
        
              3.7155e-01  -4.3769e-01 
        
              9.8810e-01  -8.7995e-02 
        
              1.2130e-01  -1.9874e-01 
        
              8.8764e-01  -7.8846e-01 
        
             -4.5774e-02   1.3533e+00 
        
             -1.4075e+00   7.6904e-01 
        
             -1.0125e+00  -3.5512e-01 
        
              5.4788e-01   5.0982e-01 
        
              6.1117e-01   6.1312e-01 
        
             -1.1951e-02  -7.4847e-01 
        
             -2.3306e-01   1.2249e+00 
        
              3.9031e-01   3.3457e-01 
        
             -6.0342e-01   8.3258e-01 
        
             -1.4474e-02   4.3668e-02 
        
             -1.0383e+00   4.7550e-01 
        
              3.0181e-01  -1.2014e+00 
        
             -1.2241e+00   1.1161e+00 
        
              2.8633e-01  -1.7207e+00 
        
              1.1247e+00  -1.2828e+00 
        
              5.7747e-01   6.3723e-02 
        
             -5.0130e-01  -7.1998e-01 
        
             -6.2614e-01   6.8299e-02 
        
             -2.6887e-01   7.0889e-01 
        
              2.0675e+00   6.7224e-01 
        
              1.4923e+00   3.5369e-01 
        
              2.5688e-01   5.8178e-01 
        
              1.7251e-01   6.6159e-01 
        
              3.0943e-01   6.7498e-01 
        
              1.6170e+00   6.2428e-01 
        
              2.9097e-01   7.5892e-01 
        
             -2.5750e-01   5.3641e-01 
        
             -1.2159e+00   5.4265e-01 
        
              1.3026e+00   1.3931e-01 
        
             -4.5703e-02   1.2355e-02 
        
              2.9373e-01   5.7735e-02 
        
              6.2447e-01  -1.6189e-01 
        
             -1.2530e+00   1.8869e+00 
        
              7.2625e-01   2.0691e+00 
        
              3.2277e-01  -5.4754e-01 
        
              1.4995e+00   6.1965e-01 
        
             -4.6988e-01  -4.1966e-01 
        
             -5.7082e-01  -1.0154e-01 
        
              2.0436e+00   5.9040e-01 
        
              2.5329e-01   5.8436e-01 
        
              4.0538e-01   1.0455e+00 
        
              6.1089e-01  -1.4342e+00 
        
             -9.1801e-01  -1.4673e+00 
        
              6.8906e-01   9.7031e-01 
        
              8.8317e-01  -1.5323e-01 
        
             -8.5984e-01   5.2180e-01 
        
              7.1558e-01   1.2908e+00 
        
             -1.2936e+00  -3.3510e-01 
        
              2.6768e-01  -4.0868e-01 
        
              4.1302e-01  -1.1250e-01 
        
             -1.5059e+00  -2.8917e-01 
        
             -2.0200e-01   1.2922e+00 
        
             -9.5919e-01  -7.9332e-01 
        
              4.1857e-01  -9.2348e-01 
        
              1.9825e-01   1.8656e-01 
        
              1.2912e-03  -2.8916e-02 
        
             -2.8855e-01  -5.2092e-01 
        
             -5.7201e-01  -7.3726e-02 
        
             -7.1932e-01   1.2260e+00 
        
             -1.3339e-01  -1.3236e-01 
        
             -5.5982e-01   1.3592e+00 
        
             -1.7448e+00  -4.2976e-01 
        
             -1.0092e+00   1.1888e+00 
        
             -4.5201e-01   4.8121e-02 
        
             -1.0593e+00   1.1674e+00 
        
              4.1674e-01   4.4604e-01 
        
             -4.8155e-01   1.3528e+00 
        
              4.8446e-01  -1.8226e-01 
        
             -6.0829e-01   4.4278e-01 
        
             -1.8806e-01   7.1140e-01 
        
             -2.0274e-01   4.2813e-01 
        
              3.4594e-01  -1.2564e+00 
        
             -2.1356e-01  -6.1439e-01 
        
              7.9506e-01   9.3215e-01 
        
             -2.0346e+00  -1.1513e+00 
        
              3.9213e-01  -2.8696e+00 
        
              1.5285e+00  -6.9886e-01 
        
              2.3145e-01   8.5767e-01 
        
              2.5474e-02  -8.8660e-01 
        
             -1.1214e-01   7.9744e-01 
        
              4.8624e-01   3.0845e-02 
        
              1.5264e-01   2.9106e-01 
        
              6.9704e-01   1.4676e-01 
        
              9.2192e-01  -1.5356e+00 
        
             -1.0916e+00  -7.2503e-01 
        
              4.6123e-01  -7.0263e-01 
        
             -7.6231e-01   1.8866e-01 
        
             -2.6786e-02   3.3630e-01 
        
             -3.1738e-01   4.1458e-01 
        
              2.5785e-01  -1.1331e+00 
        
              1.5854e+00   9.3968e-01 
        
              2.3971e-01   2.9838e-01 
        
             -8.7588e-01   9.3150e-01 
        
             -1.2594e-01   1.6150e-01 
        
             -1.1880e+00  -5.9392e-01 
        
             -3.4828e-01   1.5762e+00 
        
             -5.5121e-01   2.3123e+00 
        
              1.5885e+00  -1.1694e+00 
        
             -3.1883e-01   3.4189e-01 
        
              3.0167e-01  -2.1803e-01 
        
              8.0542e-01   1.8468e-01 
        
             -6.4677e-01   3.0375e-01 
        
             -2.8966e-01  -1.0468e+00 
        
             -1.4670e+00   3.3995e-01 
        
             -8.5143e-01  -1.0992e-01 
        
              4.0036e-01  -2.7837e-01 
        
              5.5388e-01  -1.2908e+00 
        
             -2.9182e-01   2.1048e-01 
        
             -6.1801e-01  -1.1581e-01 
        
              1.4329e-01  -3.9157e-01 
        
              1.2600e-01   3.3449e-01 
        
             -1.4600e-01   9.1032e-01 
        
              6.4881e-01  -3.2057e-01 
        
             -6.9968e-01  -1.1923e-01 
        
             -1.5710e-01  -9.3134e-01 
        
              8.0183e-02   4.5401e-01 
        
              1.9205e+00   3.7203e-01 
        
              9.7946e-01   1.1894e-01 
        
             -7.7945e-04   6.5408e-01 
        
             -3.3258e-01  -1.4887e-02 
        
             -9.0232e-01   7.6812e-01 
        
              2.6022e-01   7.4557e-01 
        
             -7.0804e-01   1.3025e+00 
        
             -1.5161e-01  -6.9582e-02 
        
             -2.3222e+00  -7.2860e-01 
        
              6.0072e-01  -1.1950e+00 
        
             -7.5329e-02  -3.2731e-01 
        
              1.2730e-01   3.6225e-01 
        
             -2.0631e-01   2.2414e-01 
        
              1.0894e+00  -1.1487e+00 
        
             -1.0216e-01   1.1211e+00 
        
              4.3208e-01  -1.6927e+00 
        
              1.3122e-01   2.6871e-01 
        
              1.7542e-01   1.9413e+00 
        
              2.5970e-01   2.3125e-02 
        
             -5.3065e-01   9.7859e-01 
        
              9.0815e-01   1.4329e-01 
        
             -1.7835e+00  -5.2359e-01 
        
              2.4299e-01  -3.0839e-01 
        
             -2.3602e-01  -3.4265e-01 
        
              2.8698e-01   9.8544e-01 
        
              1.0736e-01  -1.2778e+00 
        
             -2.9109e-01  -7.8224e-01 
        
             -1.4884e+00   6.1048e-01 
        
             -2.9043e-01   4.4939e-01 
        
             -1.2425e-01   5.2429e-01 
        
             -1.1091e+00  -9.2264e-01 
        
             -2.6662e-01   4.3608e-02 
        
              6.7411e-02  -1.3036e+00 
        
              7.3121e-01  -2.9942e-01 
        
             -7.8314e-02   9.8487e-01 
        
             -6.7113e-02   1.3408e-01 
        
             -3.6407e-01  -7.2287e-01 
        
              1.5144e+00  -7.7194e-02 
        
             -9.9624e-01  -2.4709e-01 
        
              7.8199e-01  -6.3529e-01 
        
             -4.1514e-01   1.0290e+00 
        
              2.6198e-01   2.2724e-01 
        
              1.0518e-01  -1.3066e+00 
        
              1.3108e-01   1.1973e+00 
        
             -5.7243e-01  -5.1774e-01 
        
              7.6119e-01  -1.5153e-01 
        
             -7.0079e-01  -1.2582e-01 
        
             -7.1052e-01  -1.1474e+00 
        
              2.7173e-01   7.7234e-01 
        
             -1.0089e+00   1.5093e+00 
        
              6.3435e-01  -5.4633e-01 
        
             -9.6893e-01   9.7580e-01 
        
             -2.1307e-01   4.0638e-01 
        
             -7.6041e-01  -4.6260e-01 
        
              2.2764e-02   4.9912e-01 
        
             -9.5916e-01  -1.2737e+00 
        
             -1.8180e-01  -3.1700e-01 
        
              1.9270e+00  -2.4544e-01 
        
             -4.6769e-01  -1.0749e+00 
        
             -4.8108e-01   4.8027e-02 
        
              2.5861e+00  -2.2196e-01 
        
              5.1651e-01  -9.1977e-01 
        
             -9.3873e-01  -4.9621e-01 
        
             -2.0005e+00  -8.6770e-01 
        
             -3.0755e-01   5.0605e-01 
        
              8.9145e-01   8.6706e-01 
        
              2.2512e-02   5.9619e-01 
        
             -2.8439e-03   9.3896e-01 
        
              5.4244e-01   9.3050e-01 
        
             -3.9861e-01  -8.3621e-01 
        
              6.2413e-01  -8.9209e-01 
        
              1.7433e-01  -3.9924e-01 
        
              5.7852e-01   2.2729e+00 
        
              3.6801e-01   9.7082e-01 
        
             -7.4364e-01   4.0424e-01 
        
             -6.2995e-02   1.3242e+00 
        
             -4.0151e-01   4.7026e-02 
        
             -1.5476e+00  -3.1665e-01 
        
              1.5461e+00  -2.3309e-01 
        
              5.1846e-01   3.2604e-01 
        
             -4.8491e-02   5.6833e-01 
        
             -4.6390e-01  -2.6539e-01 
        
             -4.6638e-01  -1.4529e-01 
        
              1.6006e-01   1.9965e-01 
        
             -5.7754e-01  -3.4843e-01 
        
              2.5064e+00  -1.8994e+00 
        
              1.8796e+00   2.7897e+00 
        
              2.1753e+00  -3.1012e+00 
        
              3.2423e+00   2.2498e+00 
        
              2.4233e+00   1.9034e+00 
        
              3.2928e+00   2.0893e+00 
        
              3.4098e+00   1.6821e+00 
        
              1.3400e+00   3.4924e+00 
        
              5.0302e-01  -3.2448e+00 
        
              3.3161e-01  -3.6767e+00 
        
              1.8220e-01  -3.0915e+00 
        
              3.8594e+00  -9.7802e-01 
        
              2.4418e+00   3.0404e+00 
        
              5.9453e-01   3.7994e+00 
        
              3.5970e+00  -6.1518e-01 
        
              1.1209e+00  -3.5091e+00 
        
              6.0344e-01   3.4087e+00 
        
              1.3407e+00   3.5889e+00 
        
              3.5929e+00  -4.4246e-01 
        
              3.8790e+00  -1.6489e-01 
        
              2.6651e+00   1.5094e+00 
        
              9.9814e-01   3.5185e+00 
        
              3.9643e+00   9.3112e-02 
        
              2.1964e+00   2.1655e+00 
        
              1.8028e+00  -2.5982e+00 
        
              2.6869e+00   2.4722e+00 
        
              3.2196e+00   3.7027e-01 
        
              1.5385e+00  -3.4240e+00 
        
              3.2676e+00   3.1447e-01 
        
              3.2073e+00   1.1570e+00 
        
              3.4418e+00  -4.9699e-01 
        
              2.5102e+00  -1.6744e+00 
        
              3.6915e+00  -2.6637e-02 
        
              3.5055e+00   6.7888e-01 
        
              1.2939e+00   3.7114e+00 
        
              1.7398e+00  -2.8125e+00 
        
              2.8134e+00  -1.1297e+00 
        
              2.1141e+00   2.8490e+00 
        
              3.5201e+00   7.4635e-02 
        
              2.2945e+00  -3.0125e+00 
        
              2.2593e+00   2.8871e+00 
        
              1.6613e+00   3.6236e+00 
        
              5.1963e-01  -3.6549e+00 
        
              3.3865e+00  -1.0652e+00 
        
              1.1877e+00  -3.4293e+00 
        
              3.1282e-01   3.9580e+00 
        
              2.6355e+00  -1.6796e+00 
        
              2.8577e+00   1.6887e+00 
        
              3.3582e+00   1.7489e-01 
        
              1.3943e-01   3.7145e+00 
        
              2.6274e+00  -2.7427e+00 
        
              2.2097e+00   2.3777e+00 
        
              2.8432e+00   1.5836e+00 
        
              4.9565e-01  -3.9420e+00 
        
              3.8645e-02  -3.9705e+00 
        
              1.2480e+00  -2.8542e+00 
        
              3.8542e+00  -9.6652e-01 
        
              2.1074e-01   3.8111e+00 
        
              1.4573e+00   3.3015e+00 
        
              3.2221e+00   6.2539e-01 
        
              3.3170e+00   1.7089e+00 
        
              2.9395e+00  -7.4254e-01 
        
              1.6952e+00  -3.3422e+00 
        
              2.5714e+00   2.9874e+00 
        
              2.7133e+00   2.2495e+00 
        
              2.3265e+00  -2.8315e+00 
        
              3.5208e+00  -5.8560e-01 
        
              2.5236e+00  -2.9332e+00 
        
              3.2719e-01   3.2804e+00 
        
              8.5447e-01  -2.8869e+00 
        
              2.8434e+00  -1.9127e+00 
        
              1.7349e+00  -2.7292e+00 
        
              3.3296e+00   8.9364e-01 
        
              3.4242e+00   8.7185e-01 
        
              2.6763e+00  -2.3934e+00 
        
              1.4823e+00  -3.3920e+00 
        
              2.5875e+00   1.5550e+00 
        
              2.0689e+00  -2.3129e+00 
        
              4.5555e-01  -3.6125e+00 
        
              4.9937e-01  -3.8371e+00 
        
              4.1292e-01   3.2048e+00 
        
              3.2482e+00  -1.1052e+00 
        
              3.9099e+00   7.8330e-01 
        
              1.8258e-02  -3.1769e+00 
        
              1.3208e-01  -3.7136e+00 
        
              1.9657e+00  -2.3062e+00 
        
              1.4741e-01   3.1537e+00 
        
              3.3000e+00  -1.5788e+00 
        
              1.0503e+00  -3.4513e+00 
        
              3.2239e+00   1.7308e+00 
        
              1.8663e+00   2.9013e+00 
        
              1.4533e+00  -2.8478e+00 
        
              1.0689e+00   3.0239e+00 
        
              3.3468e+00   1.0380e-01 
        
              2.7156e+00   2.5194e+00 
        
              3.0932e+00  -2.1815e-01 
        
              3.5523e+00  -3.9567e-01 
        
              6.6928e-01   3.7763e+00 
        
              1.4625e+00  -2.6852e+00 
        
              2.4077e+00   1.9012e+00 
        
              3.0820e+00   4.5203e-01 
        
              2.1356e+00  -2.1662e+00 
        
              2.1501e+00   2.4851e+00 
        
              3.2060e+00  -4.0972e-01 
        
              2.8445e+00  -2.6934e+00 
        
              1.4754e-01   3.3714e+00 
        
              2.7770e+00   2.3291e+00 
        
              2.8692e+00  -1.9879e+00 
        
              2.9093e+00   2.3058e+00 
        
              3.5688e+00   1.6988e+00 
        
              1.6169e+00   2.7511e+00 
        
              1.5958e+00   3.2149e+00 
        
              2.8247e+00   2.4592e+00 
        
              3.1434e+00  -1.0674e+00 
        
              7.8656e-01  -3.6465e+00 
        
              1.8007e+00  -3.1528e+00 
        
              1.5054e-01   3.8977e+00 
        
              2.3431e+00  -2.4848e+00 
        
              6.9042e-01   3.2131e+00 
        
              1.2733e+00  -3.7317e+00 
        
              4.8481e-01   3.2917e+00 
        
              3.3422e-01  -3.9832e+00 
        
              3.2679e+00   2.2584e+00 
        
              2.7642e+00  -2.5446e+00 
        
              1.8973e+00  -2.5222e+00 
        
              1.0259e+00   3.6717e+00 
        
              3.1450e+00  -2.4362e+00 
        
              2.5200e+00  -2.2837e+00 
        
              3.1098e+00   6.4202e-02 
        
              3.4358e+00   8.9585e-01 
        
              3.3312e-01  -3.6223e+00 
        
              7.3089e-01  -3.7199e+00 
        
              4.3466e-02   3.0581e+00 
        
              3.1891e+00   1.0840e+00 
        
              1.9716e+00   3.3469e+00 
        
              3.3584e+00  -1.5314e+00 
        
              2.4616e+00   2.5843e+00 
        
              3.6928e+00   1.3721e+00 
        
              1.4344e+00   2.9788e+00 
        
              8.5398e-01   3.7019e+00 
        
              3.3877e+00   1.7259e-01 
        
              2.6197e+00   2.9497e+00 
        
              1.9873e+00   2.5815e+00 
        
              3.3271e+00   1.2904e+00 
        
              3.1860e+00   2.2036e+00 
        
              3.2937e+00  -1.6439e+00 
        
              7.2128e-01  -3.1683e+00 
        
              3.2746e+00   1.4682e+00 
        
              1.7367e+00   2.8936e+00 
        
              3.9239e+00  -7.3084e-01 
        
              8.9442e-01  -3.6034e+00 
        
              8.3828e-01  -2.9414e+00 
        
              1.1226e+00  -3.2601e+00 
        
              3.7244e+00  -4.3417e-01 
        
              1.5585e+00   2.6746e+00 
        
              2.1768e+00  -2.4645e+00 
        
              2.7365e+00   2.1054e+00 
        
              3.5530e+00   1.5495e+00 
        
              2.9085e+00   2.5095e+00 
        
              2.2721e+00  -2.7614e+00 
        
              3.3224e+00  -1.8044e+00 
        
              3.3618e+00   1.3507e+00 
        
              3.2811e+00  -6.3843e-01 
        
              3.0514e+00  -1.6722e-01 
        
              2.6975e+00  -2.6324e+00 
        
              7.3917e-01   3.2963e+00 
        
              2.0482e+00   2.3292e+00 
        
              3.0462e+00  -1.9496e+00 
        
              3.2474e+00  -2.7507e-02 
        
              2.4865e+00   3.1105e+00 
        
              1.7824e+00  -3.4810e+00 
        
              3.3282e+00   9.5392e-01 
        
              2.2863e+00  -2.0887e+00 
        
              1.0727e+00  -3.6550e+00 
        
              3.5058e+00   2.4701e-01 
        
              3.1500e+00   2.1377e+00 
        
              2.9694e+00   5.8342e-01 
        
              3.6459e+00   2.6234e-01 
        
              2.7560e+00  -1.9036e+00 
        
              3.8862e+00  -2.8632e-01 
        
              1.8027e+00   3.4033e+00 
        
              1.0073e+00   2.9809e+00 
        
              1.9742e+00  -2.5839e+00 
        
              3.4128e+00   1.6383e+00 
        
              2.3616e+00  -2.5208e+00 
        
              3.5011e+00   1.4275e+00 
        
              1.1735e+00   3.6123e+00 
        
              3.2781e+00  -1.6373e-02 
        
              3.7625e+00  -8.7912e-01 
        
              7.8129e-01   3.9224e+00 
        
              1.0445e-01  -3.8023e+00 
        
              2.7070e+00   1.6086e+00 
        
              3.2320e+00   1.4832e+00 
        
              2.6423e+00  -1.9589e+00 
        
              2.9778e+00  -2.4839e+00 
        
              1.9220e+00  -3.2348e+00 
        
              3.5437e+00   1.2368e+00 
        
              8.8257e-01   3.8494e+00 
        
              2.0639e+00   2.9133e+00 
        
              2.5943e+00  -1.7020e+00 
        
              2.6024e+00  -2.4817e+00 
        
              2.8694e+00  -2.6930e+00 
        
              3.9516e+00  -3.0726e-01 
        
              3.6030e+00  -6.1647e-01 
        
              2.0183e+00   3.2818e+00 
        
              3.5283e+00  -3.9476e-02 
        
              3.1129e+00   4.2291e-01 
        
              2.1280e+00  -3.1000e+00 
        
              3.3317e+00   1.6095e+00 
        
              3.2273e+00   5.9118e-01 
        
              2.4380e+00  -2.4406e+00 
        
              2.4699e+00  -2.9463e+00 
        
              2.2945e+00   2.4440e+00 
        
              9.1557e-02  -3.9608e+00 
        
              1.3583e+00  -2.8671e+00 
        
              4.0359e-01  -3.3651e+00 
        
              3.2635e+00   1.9620e+00 
        
              2.4965e+00  -2.1700e+00 
        
              2.7656e+00   1.3538e+00 
        
              1.9656e+00  -3.2959e+00 
        
              2.7649e+00  -2.3874e+00 
        
              3.2658e-01  -3.9351e+00 
        
              9.9548e-01   3.0591e+00 
        
              8.1040e-01  -3.5035e+00 
        
              8.5208e-01  -2.9588e+00 
        
              3.8297e+00   1.0967e-01 
        
              1.8124e+00  -3.2696e+00 
        
              2.0394e+00  -2.6714e+00 
        
              3.0270e+00   9.1670e-01 
        
              2.8553e-02  -3.8588e+00 
        
              2.5301e+00   1.6696e+00 
        
              3.4365e+00  -6.0862e-01 
        
              1.1126e+00  -3.0184e+00 
        
              7.2068e-01   3.1825e+00 
        
              2.9699e+00   7.1126e-01 
        
              2.2811e+00  -1.9724e+00 
        
              7.2441e-02  -3.7541e+00 
        
              3.2042e+00  -1.6911e+00 
        
              1.9610e+00   2.9393e+00 
        
              3.5195e+00  -1.4951e+00 
        
              1.5605e+00   2.8595e+00 
        
              1.0510e+00   3.6493e+00 
        
              3.5461e+00   1.0493e+00 
        
              2.8276e+00   1.0332e+00 
        
              3.7304e+00   1.3005e+00 
        
              9.6439e-01   3.2598e+00 
        
              1.4995e+00   2.6443e+00 
        
              1.7180e+00   2.9554e+00 
        
              3.3394e+00  -3.9563e-01 
        
              1.4999e+00  -3.3477e+00 
        
              2.9676e+00   1.9276e+00 
        
              1.0324e+00  -3.7712e+00 
        
              7.1990e-01  -3.7543e+00 
        
              2.6063e+00   1.7879e+00 
        
              3.2948e+00   2.3942e-01 
        
              8.0931e-01   3.0524e+00 
        
              3.3258e+00  -4.7174e-01 
        
              1.9011e+00   2.3970e+00 
        
              3.4902e+00  -1.3928e+00 
        
              3.8354e+00  -1.1240e+00 
        
              8.8554e-01   3.5922e+00 
        
              2.3322e+00  -2.3546e+00 
        
              5.0670e-01  -3.3023e+00 
        
              3.2252e-01   3.8157e+00 
        
              1.1884e+00   2.9691e+00 
        
              1.9961e+00   2.5945e+00 
        
              2.9979e+00   1.5755e+00 
        
              1.4056e+00  -3.3139e+00 
        
              3.2239e+00  -5.9160e-01 
        
              2.9049e-02  -3.0027e+00 
        
              3.4695e+00   1.5536e+00 
        
              2.6766e+00   1.3715e+00 
        
              1.6287e+00  -2.7554e+00 
        
              8.4711e-01   3.4302e+00 
        
              1.5241e+00  -2.8215e+00 
        
              3.8257e+00   1.0343e+00 
        
              2.9502e+00   1.8943e+00 
        
              9.0620e-01   2.9154e+00 
        
              3.3868e+00   2.0231e+00 
        
              3.2486e+00   6.0980e-01 
        
              2.7731e-01  -3.4781e+00 
        
              2.4410e+00   2.5972e+00 
        
              3.6845e+00  -1.2832e+00 
        
              3.6547e+00  -5.9191e-01 
        
              2.3732e+00  -2.7207e+00 
        
              1.7777e+00   2.4539e+00 
        
              9.3953e-01   3.4290e+00 
        
              2.7743e+00   2.1890e+00 
        
              9.7148e-02   3.7860e+00 
        
              1.5750e+00  -2.7746e+00 
        
              2.5651e+00   2.3856e+00 
        
              1.9695e+00  -3.2011e+00 
        
              2.2635e+00  -2.0502e+00 
        
              2.4135e+00  -2.2269e+00 
        
              3.5860e-01  -3.4283e+00 
        
              3.4557e+00  -1.3020e+00 
        
              1.5625e-01  -3.1063e+00 
        
              3.5074e+00  -1.8139e+00 
        
              3.1180e+00   5.1707e-02 
        
              3.3736e+00   3.5003e-01 
        
              3.1920e+00  -1.0053e+00 
        
              1.4406e+00  -2.9575e+00 
        
              1.6117e+00   3.2349e+00 
        
              1.5411e+00   3.4321e+00 
        
              5.2396e-01  -3.2167e+00 
        
              4.6195e-01   3.3701e+00 
        
              3.5753e+00  -1.5894e+00 
        
              7.0689e-03   3.0558e+00 
        
              3.8859e+00   5.3956e-01 
        
              1.7532e+00   2.4984e+00 
        
              2.3481e+00   2.9376e+00 
        
              3.0615e+00   1.7630e+00 
        
              9.2096e-01  -2.8662e+00 
        
              1.5273e+00  -3.5901e+00 
        
              3.0061e+00  -5.9953e-01 
        
              2.0830e+00   2.5680e+00 
        
              2.9102e+00  -2.6058e+00 
        
              2.6168e+00  -1.5470e+00 
        
              9.1578e-01   3.5705e+00 
        
              3.2631e+00  -1.8946e-01 
        
              1.9386e+00  -2.5159e+00 
        
              1.3709e+00   3.4482e+00 
        
              3.1303e+00   1.3016e+00 
        
              3.2657e+00  -9.2408e-01 
        
              5.5065e-01   3.5064e+00 
        
              8.2463e-01   3.3932e+00 
        
              4.5411e-01   3.5487e+00 
        
              2.9182e+00  -1.8959e+00 
        
              3.1259e+00   1.0038e-01 
        
              6.9314e-01  -3.7521e+00 
        
              2.8149e+00   1.1239e+00 
        
              2.3563e+00  -2.2201e+00 
        
              2.7031e+00  -1.8030e+00 
        
              1.1006e+00  -3.7385e+00 
        
              1.5132e+00  -3.1552e+00 
        
              1.5252e+00  -2.9858e+00 
        
              3.2626e+00  -2.1708e+00 
        
              1.4854e+00   3.1607e+00 
        
              3.1205e+00   1.5873e+00 
        
              2.3807e+00   2.5499e+00 
        
              1.4191e+00  -3.3969e+00 
        
              2.1637e+00   2.3969e+00 
        
              3.5766e+00   3.1294e-01 
        
              1.5803e+00   3.1561e+00 
        
              3.1436e+00   1.6469e+00 
        
              3.4251e+00  -7.8622e-01 
        
              2.2618e+00   2.8274e+00 
        
              2.4409e+00  -2.7737e+00 
        
              2.7358e+00   2.7505e+00 
        
              3.3272e+00  -7.2573e-01 
        
              2.4268e+00   2.6684e+00 
        
              1.5434e+00  -2.8522e+00 
        
              3.0459e+00  -1.7342e+00 
        
              2.9497e+00   2.4843e+00 
        
              2.8535e+00  -1.0976e+00 
        
              9.0001e-01  -3.5770e+00 
        
              4.7459e-01  -3.2130e+00 
        
              2.0746e+00  -2.3376e+00 
        
              2.0328e+00   2.3414e+00 
        
              2.1990e+00   2.8137e+00 
        
              3.3143e+00   1.6692e+00 
        
              2.5338e+00   2.5043e+00 
        
              2.4769e+00   1.8473e+00 
        
              2.8596e+00  -1.5264e+00 
        
              2.0440e+00   3.0064e+00 
        
              2.7127e+00  -2.8693e+00 
        
              3.6661e+00   5.2810e-01 
        
              3.2483e+00  -1.0199e+00 
        
              3.1396e+00   2.8018e-01 
        
              3.2530e-01   3.1358e+00 
        
              3.6527e+00  -1.4630e-01 
        
              9.8621e-01  -2.8706e+00 
        
              3.3159e+00  -7.0739e-01 
        
              8.4547e-01   3.5865e+00 
        
              6.9317e-01   3.8918e+00 
        
              8.0404e-01  -3.3257e+00 
        
              1.9650e+00   2.5873e+00 
        
              3.4365e+00   1.8219e+00 
        
              8.9074e-01   3.0885e+00 
        
              1.3192e+00   3.3725e+00 
        
              1.6933e+00  -2.9098e+00 
        
              2.9183e+00   2.6368e+00 
        
              2.1467e+00  -3.1101e+00 
        
              3.5106e+00   2.6571e-01 
        
              3.7817e+00  -5.1672e-01 
        
              3.9366e-01   3.6315e+00 
        
              3.7876e+00   9.9987e-01 
        
              4.6316e-01   3.8581e+00 
        
              2.8852e+00   1.8334e+00 
        
              3.6932e+00  -1.1497e+00 
        
              1.1934e+00  -2.8704e+00 
        
              2.9969e+00  -3.0917e-01 
        
              5.6266e-01   3.6611e+00 
        
              1.0397e+00   2.9162e+00 
        
              2.9713e+00  -2.5004e+00 
        
              3.0736e+00  -2.1952e+00 
        
              2.1191e+00  -3.0135e+00 
        
              3.2304e+00  -1.1271e+00 
        
              4.2136e-01   3.9310e+00 
        
              3.3863e+00   8.3441e-01 
        
              1.5903e-01  -3.7851e+00 
        
              8.6148e-01   3.3707e+00 
        
              2.4500e+00  -2.3769e+00 
        
              3.0894e+00  -4.7063e-01 
        
              2.9916e+00  -1.6288e+00 
        
              3.4119e+00  -1.2034e+00 
        
              2.9100e+00  -2.1200e+00 
        
              3.8424e-01  -3.8942e+00 
        
              3.7757e+00   4.7166e-01 
        
              8.7066e-01   3.2156e+00 
        
              2.7974e+00  -2.2442e+00 
        
              1.0020e+00   3.5294e+00 
        
              2.9193e+00   9.6362e-01 
        
              1.4394e+00  -3.2835e+00 
        
              1.8226e-01  -3.4827e+00 
        
              2.3539e+00   2.2029e+00 
        
              5.3804e-01  -3.8751e+00 
        
              1.9966e+00  -2.2741e+00 
        
              1.7686e+00  -2.9067e+00 
        
              3.9009e+00  -1.8150e-01 
        
              2.3992e+00   2.0427e+00 
        
              3.2850e+00   5.6488e-01 
        
              2.3042e+00  -2.5523e+00 
        
              2.3357e+00  -1.9213e+00 
        
              8.6909e-01  -3.0858e+00 
        
              1.1188e+00   3.8194e+00 
        
              2.9321e+00  -1.0424e+00 
        
              1.1160e+00   3.3635e+00 
        
              2.5811e+00   2.4498e+00 
        
              4.6646e-01   3.5346e+00 
        
              2.0842e+00   3.3019e+00 
        
              2.8357e+00  -1.8114e+00 
        
              1.5592e+00  -3.6037e+00 
        
              2.2012e+00   2.6138e+00 
        
              7.2512e-01   3.6860e+00 
        
              7.5530e-02  -3.9873e+00 
        
              8.2941e-01   3.3298e+00 
        
              1.0345e+00   3.4456e+00 
        
              2.5919e+00  -1.7574e+00 
        
              3.9705e+00  -2.2167e-01 
        
              3.7847e+00   1.0189e+00 
        
              2.1092e+00  -2.1520e+00 
        
              3.5445e+00  -7.6223e-01 
        
              3.7318e+00   1.3234e+00 
        
              3.8866e+00  -2.9519e-01 
        
              3.7695e+00  -2.8410e-01 
        
              1.7586e+00   2.9416e+00 
        
              3.4053e+00  -1.4835e+00 
        
              2.0183e+00   2.3202e+00 
        
              2.5091e+00  -2.5788e+00 
        
              2.1403e+00   2.9343e+00 
        
              2.0782e+00   2.6681e+00 
        
              3.4793e+00  -3.6475e-01 
        
              1.1006e+00   3.2158e+00 
        
              2.8225e+00  -1.0203e+00 
        
              5.5829e-01   3.7016e+00 
        
              2.0665e+00  -2.6634e+00 
        
              1.7900e+00  -3.2289e+00 
        
              2.5794e+00   1.8142e+00 
        
              1.8880e+00  -2.4712e+00 
        
              2.7883e-01   3.0388e+00 
        
              3.2136e+00   1.4619e+00 
        
              2.0484e+00  -2.3879e+00 
        
              2.4439e+00  -2.3100e+00 
        
              7.3675e-01  -3.7768e+00 
        
              5.3547e-01  -3.0894e+00 
        
              9.3167e-02  -3.2047e+00 
        
              1.5647e+00  -2.8067e+00 
        
              3.1307e+00   8.5950e-01 
        
              9.5680e-01  -3.4374e+00 
        
              7.7619e-01  -3.7667e+00 
        
              3.5295e+00   4.3717e-01 
        
              1.6903e+00   2.7718e+00 
        
              2.8744e+00  -1.4138e+00 
        
              3.5804e+00   2.3334e-01 
        
              3.5198e+00   3.5049e-01 
        
              1.3898e+00   2.8153e+00 
        
              1.6927e+00   2.5325e+00 
        
              7.1890e-01  -3.4086e+00 
        
              2.4239e+00  -2.7716e+00 
        
              3.0268e+00  -5.2831e-01 
        
              1.7130e+00  -3.5618e+00 
        
              3.2181e+00  -6.3201e-01 
        
              6.3490e-01   3.6426e+00 
        
              2.5348e+00  -1.8521e+00 
        
              1.1998e+00  -3.7480e+00 
        
              4.6597e-01  -3.6956e+00 
        
              3.6673e+00   7.7248e-01 
        
              3.3341e-01   3.3943e+00 
        
              2.6093e+00   1.9744e+00 
        
              7.8294e-01   3.1852e+00 
        
              3.6653e+00  -1.1211e+00 
        
              8.8452e-01  -3.0822e+00 
        
              2.2008e+00  -2.5728e+00 
        
              1.7767e+00   3.5169e+00 
        
              3.7772e+00   1.2704e+00 
        
              5.2182e-01  -3.3297e+00 
        
              2.2057e+00   2.0449e+00 
        
              3.0381e+00   7.5620e-01 
        
              3.4052e+00  -1.2178e+00]; 
        
           % Run dbscan clustering algorithm 
        
           epsilon = 0.3; 
        
           min_pts = 3; 
        
           [index] = mi.dbscan(X, epsilon, min_pts); 
        
           % Get number of clusters 
        
           clusters = max(index); 
        
           legends = {}; 
        
           for i = 0:clusters 
        
             % Get the clusters 
        
             Xi = X(find(index == i), :); 
        
             % Create legends 
        
             if i~=0 
        
               legends{end+1} = ['Cluster #' num2str(i)]; 
        
              else 
        
               if ~isempty(Xi) 
        
                 legends{end+1} = 'Noise'; 
        
               end 
        
             end 
        
             % Plot 
        
             if ~isempty(Xi) 
        
               scatter(Xi(:,1),Xi(:,2)); 
        
             end 
        
             hold on; 
        
           end 
        
           grid on; 
        
           legend(legends); 
        
           legend('Location', 'NorthEastOutside');

Results

ERA/DC - Eigensystem Realization Algorithm Data Correlations

ERA/DC was invented 1987 and is a successor from ERA, that was invented 1985 at NASA. The difference between ERA/DC and ERA is that ERA/DC can handle noise much better than ERA. But both algorihtm works as the same. ERA/DC want an impulse response. e.g called markov parameters. You will get a state space model from this algorithm. This algorithm can handle both SISO and MISO data.

Use this algorithm if you got impulse data from e.g structural mechanics.

[sysd, K] = mi.eradc(g, sampleTime, ktune, delay systemorder);

Example ERA/DC for MIMO systems

MataveID/examples/eradcExample.m

Lines 1 to 61 in 2014b74

    
           clc; clear; close all 
        
           %% Parameters 
        
           m1 = 2.3; 
        
           m2 = 3.1; 
        
           k1 = 8.5; 
        
           k2 = 5.1; 
        
           b1 = 3.3; 
        
           b2 = 5.1; 
        
           A=[0                 1   0                                              0 
        
             -(b1*b2)/(m1*m2)   0   ((b1/m1)*((b1/m1)+(b1/m2)+(b2/m2)))-(k1/m1)   -(b1/m1) 
        
              b2/m2             0  -((b1/m1)+(b1/m2)+(b2/m2))                      1 
        
              k2/m2             0  -((k1/m1)+(k1/m2)+(k2/m2))                      0]; 
        
           B=[0 0; 
        
              1/m1 0; 
        
              0 0 ; 
        
              (1/m1)+(1/m2) 1/m2]; 
        
           C=[0 0 1 0; 
        
              0 1 0 0]; 
        
           D=[0 0; 
        
              0 0]; 
        
           delay = 0; 
        
           %% Model 
        
           buss = mc.ss(delay,A,B,C,D); 
        
           %% Simulation 
        
           [g, t] = mc.impulse(buss, 10); 
        
           % Reconstruct the input impulse signal from impulse.m 
        
           u = zeros(size(g)); 
        
           u(1) = 1; 
        
           %% Add 15% noise 
        
           v = 2*randn(1, 1000); 
        
           for i = 1:length(g)-1 
        
             noiseSigma = 0.15*g(i); 
        
             noise = noiseSigma*v(i); % v = noise, 1000 samples -1 to 1 
        
             g(i) = g(i) + noise; 
        
           end 
        
           %% Identification 
        
           systemorder = 10; 
        
           ktune = 0.09; 
        
           sampleTime = t(2) - t(1); 
        
           delay = 0; 
        
           [sysd, K] = mi.eradc(g, sampleTime, ktune, delay, systemorder); 
        
           % Create the observer 
        
           observer = mc.ss(sysd.delay, sysd.A - K*sysd.C, [sysd.B K], sysd.C, [sysd.D sysd.D*0]); 
        
           observer.sampleTime = sysd.sampleTime; 
        
           %% Validation 
        
           [gf, tf] = mc.lsim(observer, [u; g], t); 
        
           close 
        
           %% Check 
        
           plot(t, g, tf, gf) 
        
           legend('Data 1', 'Data 2', 'Identified 1', 'Identified 2', 'location', 'northwest') 
        
           grid on

#Feature from accelerated segmentation test Use Feature from accelerated segmentation test(FAST) if you want to find corners inside an image. There is also an equivalent C-code FAST algorithm inside the CControl repository.

[corners, scores] = mi.fast(image, threshold, fast_method);

Example

MataveID/examples/fastExample.m

Lines 1 to 20 in cbd69ae

    
           % Close all 
        
           close all 
        
           % Read image 
        
           X = imread(fullfile('..', 'data', 'face.jpg')); 
        
           % Make image greyscale 
        
           if size(X, 3) > 1 
        
             X =  rgb2gray(X); 
        
           end 
        
           % Compute fast 
        
           threshold = 50; 
        
           fast_method = 9; 
        
           coordintes = mi.fast(X, threshold, fast_method); 
        
           % Show 
        
           imshow(uint8(X)); 
        
           hold on 
        
           plot(coordintes(:,1), coordintes(:,2), 'r.')

Result

Filtfilt - Zero Phase Filter

This filter away noise with a good old low pass filter that are being runned twice. Filtfilt is equal to the famous function filtfilt in MATLAB, but this is a regular .m file and not a C/C++ subroutine. Easy to use and recommended.

[y] = mi.filtfilt(y, t, K);

Filtfilt Example

MataveID/examples/filtfiltExample.m

Lines 1 to 31 in 2014b74

    
           clc; clear; close all; 
        
           %% Model of a mass spring damper system 
        
           M = 1; % Kg 
        
           K = 500; % Nm/m 
        
           b = 3; % Nm/s^2 
        
           G = mc.tf([1], [M b K]); 
        
           %% Input signal 
        
           t = linspace(0.0, 100, 3000); 
        
           u = 10*sin(t); 
        
           %% Simulation 
        
           y = mc.lsim(G, u, t); 
        
           close 
        
           %% Add 10% noise 
        
           v = 2*randn(1, length(y)); 
        
           for i = 1:length(y) 
        
             noiseSigma = 0.10*y(i); 
        
             noise = noiseSigma*v(i);  
        
             y(i) = y(i) + noise; 
        
           end 
        
           %% Filter away the noise 
        
           lowpass = 0.2; 
        
           [yf] = mi.filtfilt(y, t, lowpass); 
        
           %% Check 
        
           plot(t, yf, t, y); 
        
           legend("Filtered", "Noisy");

Hough Transform

Use this algorithm if you want to find lines inside an edge image. Important that the image needs to be an edge image.

[N, K, M] = mi.hough(X, p, epsilon, min_pts);

Hough Transform example

Assume that we have road that we want to track by writing two parallell lines that follows the road and we want to avoid everything else.

MataveID/examples/houghExample.m

Lines 1 to 30 in 0152b6b

    
           % Clear 
        
           clear all 
        
           close all 
        
           clc 
        
           % Read image 
        
           X = imread(fullfile('..', 'data', 'test_hough.png')); 
        
           % If the image is color 
        
           if(size(X, 3) > 1) 
        
             X = rgb2gray(X); 
        
           end 
        
           % Plot the image 
        
           imshow(X); 
        
           % Do hough transform 
        
           p = 0.3; % Percentage definition of a line e.g all lines shorter than p times longest line, should be classes as a non-line 
        
           epsilon = 10; % Minimum radius for hough cluster 
        
           min_pts = 2; % Minimum points for hough cluster 
        
           [N, K, M] = mi.hough(X, p, epsilon, min_pts); 
        
           % Plot the lines together with the image 
        
           [~, n] = size(X); 
        
           x = linspace(0, n); 
        
           hold on 
        
           for i = 1:N 
        
             y = K(i)*x + M(i); 
        
             plot(x, y) 
        
           end

Independent Component Analysis

Independent component analysis(ICA) is a tool if you want to separate independent signals from each other. This is not a filter algorithm, but instead of removing noise, it separate the disturbances from the signals. The disturbances are created from other signals. Assume that you have an engine and you are measuring vibration in X, Y and Z-axis. These axis will affect each other and therefore the signals will act like they are mixed. ICA separate the mixed signals into clean and independent signals.

[S] = mi.ica(X);

Independent Component Analysis example

MataveID/examples/icaExample.m

Lines 1 to 137 in 2014b74

    
           clear; close all; clc 
        
           % Tick clock 
        
           tic 
        
           %% Parameters 
        
           N = 6;                            %The number of observed mixtures 
        
           M = 1000;                         %Sample size, i.e.: number of observations 
        
           K = 0.1;                          %Slope of zigzag function 
        
           na = 8;                           %Number of zigzag oscillations within sample 
        
           ns = 5;                           %Number of alternating step function oscillations within sample 
        
           finalTime = 40*pi;                %Final sample time (s) 
        
           initialTime = 0;                  %Initial sample time (s) 
        
           %% Generating Data for ICA 
        
           % Create time vector data 
        
           timeVector  = initialTime:(finalTime-initialTime)/(M-1):finalTime; 
        
           % Create random, cos, sin and fast cos signal 
        
           source1     = rand(1, M);                    
        
           source2     = cos(0.25*timeVector);          
        
           source3     = sin(0.1*timeVector);           
        
           source4     = cos(0.7*timeVector);           
        
           % Ziggsack signal 
        
           source5 = zeros(1,M);    
        
           periodSource5 = (finalTime-initialTime)/na;  
        
           for i = 1:M 
        
               source5(i) = K*timeVector(i)-floor(timeVector(i)/periodSource5)*K*periodSource5; 
        
           end 
        
           source5 = source5 - mean(source5); 
        
           % PWM signal 
        
           source6 = zeros(1,M);   
        
           periodSource6 = (finalTime-initialTime)/ns/2; 
        
           for i = 1:M 
        
               if mod(floor(timeVector(i)/periodSource6),2) == 0 
        
                   source6(i) = 1; 
        
               else 
        
                   source6(i) = -1; 
        
               end     
        
           end 
        
           source6 = source6 - mean(source6); 
        
           % Create our source matrix. This matrix is what want to find 
        
           S = [source1;source2;source3;source4;source5;source6]; 
        
           % Create an matrix A that going to mix all signals in S, that we calling X 
        
           Amix = rand(N,N);                     
        
           X = Amix*S;                                  
        
           figure 
        
           plot(timeVector,source1)                     
        
           xlabel('time (s)') 
        
           ylabel('Signal Amplitude')  
        
           legend('source 1') 
        
           figure 
        
           plot(timeVector,source2)                     
        
           xlabel('time (s)') 
        
           ylabel('Signal Amplitude')  
        
           legend('source 2') 
        
           figure 
        
           plot(timeVector,source3)                     
        
           xlabel('time (s)') 
        
           ylabel('Signal Amplitude')  
        
           legend('source 3') 
        
           figure 
        
           plot(timeVector,source4)                     
        
           xlabel('time (s)') 
        
           ylabel('Signal Amplitude')  
        
           legend('source 4') 
        
           figure 
        
           plot(timeVector,source5)                    
        
           xlabel('time (s)') 
        
           ylabel('Signal Amplitude')  
        
           legend('source 5') 
        
           figure 
        
           plot(timeVector,source6)                     
        
           xlabel('time (s)') 
        
           ylabel('Signal Amplitude')  
        
           legend('source 6') 
        
           figure 
        
           plot(timeVector,X);                       
        
           xlabel('time (s)')  
        
           ylabel('Signal Amplitude') 
        
           legend('Observed Mixture 1', 'Observed Mixture 2', 'Observed Mixture 3', 'Observed Mixture 4', 'Observed Mixture 5', 'Observed Mixture 6') 
        
           % Use ICA to find S from X 
        
           S = mi.ica(X); 
        
           figure 
        
           plot(timeVector, S(1,:)) 
        
           xlabel('time (s)')  
        
           ylabel('Signal Amplitude')  
        
           legend('Source Estimation 1') 
        
           figure 
        
           plot(timeVector, S(2,:)) 
        
           xlabel('time (s)')  
        
           ylabel('Signal Amplitude')  
        
           legend('Source Estimation 2') 
        
           figure 
        
           plot(timeVector, S(3,:)) 
        
           xlabel('time (s)') 
        
           ylabel('Signal Amplitude')  
        
           legend('Source Estimation 3') 
        
           figure 
        
           plot(timeVector, S(4,:)) 
        
           xlabel('time (s)') 
        
           ylabel('Signal Amplitude')  
        
           legend('Source Estimation 4') 
        
           figure 
        
           plot(timeVector, S(5,:)) 
        
           xlabel('time (s)') 
        
           ylabel('Signal Amplitude')  
        
           legend('Source Estimation 5') 
        
           figure 
        
           plot(timeVector, S(6,:)) 
        
           xlabel('time (s)') 
        
           ylabel('Signal Amplitude')  
        
           legend('Source Estimation 6') 
        
           % End clock time and check the difference how long it took 
        
           toc

These signals are what we want to find

This is how the signals look when we are measuring them

This is how the signals are reconstructed as they were independent

IDBode - Identification Bode

This plots a bode diagram from measurement data. It can be very interesting to see how the amplitudes between input and output behaves over frequencies. This can be used to confirm if your estimated model is good or bad by using the bode command from Matavecontrol and compare it with idebode.

mi.idbode(u, y, w);

IDBode Example

MataveID/examples/idbodeExample.m

Lines 1 to 25 in 2014b74

    
           clc; clear; close all; 
        
           %% Model of a mass spring damper system 
        
           M = 5; % Kg 
        
           K = 100; % Nm/m 
        
           b = 52; % Nm/s^2 
        
           G = mc.tf([1], [M b K]); 
        
           %% Frequency response 
        
           t = linspace(0, 50, 3000); 
        
           [u, fs] = mc.chirp(t); 
        
           %% Simulation 
        
           y = mc.lsim(G, u, t); 
        
           close all 
        
           % Add noise 
        
           y = y + 0.0001*randn(1, length(y)); 
        
           %% Identify bode diagram 
        
           mi.idbode(u, y, fs); 
        
           %% Check 
        
           mc.bode(G);

Gaussian 2D filter

Use this filter if you want to blur an image

[Y] = mi.imgaussfilt(image, sigma);

Example

MataveID/examples/imgaussfiltExample.m

Lines 1 to 26 in 662c95d

    
           % Close all 
        
           close all 
        
           % Read image 
        
           image = imread(fullfile('..', 'data', 'campera.png')); 
        
           % Turn to grey scale 
        
           if(size(image, 3) == 3) 
        
             image = rgb2gray(image); 
        
           end 
        
           % Select sigma 
        
           sigma = 3; 
        
           % Compute gradients and orientations 
        
           Y = mi.imgaussfilt(image, sigma); 
        
           % Show original image 
        
           subplot(1, 2, 1); 
        
           imshow(image); 
        
           title('Original image'); 
        
           % Show gaussian image 
        
           subplot(1, 2, 2); 
        
           imshow(uint8(Y)); 
        
           title('Gaussian image');

Result

K-means clustering

K-means clustering is a tool that can identify the center of clusters. All you need to do is to specify how many cluster IDs you think there exist in your data. Use this algorithm if your data is gaussian and you know the numbers of clusters. All you want to know are the cetrums of the clusters.

[idx, C, success] = mi.kmeans(X, k);

K-means clustering example

MataveID/examples/kmeansExample.m

Lines 1 to 33 in cefeca3

    
           % Remove 
        
           clear  
        
           close all 
        
           clc 
        
           % Create data 
        
           t = linspace(0, 3*pi, 200)'; 
        
           data = [40 + 10*randn(200,3); 
        
                   50 + 5*sin(t) + 5*t, 10*randn(200, 1), 5*sqrt(t.^2); 
        
                   -20 + 23*rand(300, 3)]; 
        
           % Amount of clusters 
        
           K = 3; 
        
           % K-means clustering 
        
           [idx, C, success] = mi.kmeans(data, K); 
        
           % Check 
        
           if(success) 
        
             disp('K-means clustering success!'); 
        
           else 
        
             disp('You need to try with another K-value'); 
        
           end 
        
           % Plot cluster 
        
           figure; 
        
           scatter3(data(:,1), data(:,2), data(:,3), 10, idx, 'filled'); 
        
           hold on; 
        
           scatter3(C(:,1), C(:,2), C(:,3), 50, 'r', 'filled'); 
        
           title('K-means clustering', 'FontSize', 20); 
        
           xlabel('x', 'FontSize', 20); 
        
           ylabel('y', 'FontSize', 20); 
        
           zlabel('z', 'FontSize', 20);

Kernel Principal Component Analysis

Kernel Principal Component Analysis can be used for dimension reduction and projection on maximum variance between classes. Kernel methods make PCA suitable for nonlinear data. Kernels has proven very good results in nonlinear dimension reduction.

[P, W] = mi.kpca(X, c, kernel_type, kernel_parameters);

Kernel Principal Component Analysis example

MataveID/examples/kpcaExample.m

Lines 1 to 62 in e4274ae

    
           % First clear and close all figures 
        
           clear 
        
           close all 
        
           % Create data 
        
           X = [-3244 5324 1345; 
        
               10 30 100; 
        
               20 40 93; 
        
               30 60 163; 
        
               60 100 126; 
        
               55 13 134; 
        
               306 34 104; 
        
               316 56 120; 
        
               326 74 127; 
        
               337 80 128; 
        
               347 89 131; 
        
               358 103 139; 
        
               -31 -56 -120; 
        
               -32 -74 -127; 
        
               -33 -80 -128; 
        
               -34 -89 -131; 
        
               -35 -103 -139; 
        
               700 600 500; 
        
               1000 1000 1000; 
        
               -3231 4345 -4352; 
        
               -2342 4356 3453; 
        
               -2364 4326 3353; 
        
               658 143 1692]; 
        
           % Noisy data 
        
           X = [X; 10*randn(10, 3)]; 
        
           % Remove noise 
        
           y = mi.dbscan(X, 50, 2); 
        
           X = X(y ~= 0, :); 
        
           % Plot original data 
        
           cmap = jet(length(X)); 
        
           scatter3(X(:, 1), X(:, 2), X(:, 3), 50,cmap) 
        
           title('Original 3D data', 'FontSize', 20) 
        
           % Do KPCA 
        
           c = 2; 
        
           kernel_type = 'polynomial'; 
        
           kernel_parameters = [1, 2]; 
        
           [P, W] = mi.kpca(X, c, kernel_type, kernel_parameters); 
        
           % Reconstruct 
        
           X_reconstructed = W * P; 
        
           reconstruction_error = norm(X - X_reconstructed, 'fro') / norm(X, 'fro') 
        
           figure 
        
           switch(c) 
        
           case 3 
        
             scatter3(X_reconstructed(:, 1), X_reconstructed(:, 2), X_reconstructed(:, 3)); 
        
           case 2 
        
             scatter(X_reconstructed(:, 1), X_reconstructed(:, 2)); 
        
           case 1 
        
             scatter(X_reconstructed(:, 1), 0*X_reconstructed(:, 1)); 
        
           end 
        
           grid on 
        
           title(sprintf('Reconstruction %iD with error %f', c, reconstruction_error), 'FontSize', 20)

SRA - Stochastic Realization Algorithm

This is an algorithm that can identify a stochastic model from error measurement data.

[sysd, K] = mi.sra(e, k, sampleTime, ktune, delay, systemorder);

Example SRA 1

MataveID/examples/sraExample1.m

Lines 1 to 52 in 72c26ea

    
           % Clear all 
        
           clear all 
        
           % Create system model 
        
           G = mc.tf(1, [1 1.5 2]); 
        
           % Create disturbance model 
        
           H = mc.tf([2 3], [1 5 6]); 
        
           % Create input signal 
        
           [u, t] = mc.gensig('square', 10, 10, 100); 
        
           u = [u*5 u*2 -u 10*u -2*u]; 
        
           t = linspace(0, 30, length(u)); 
        
           % Create disturbance signal 
        
           e = randn(1, length(t)); 
        
           % Simulate with noise- 
        
           y = mc.lsim(G, u, t) + mc.lsim(H, e, t); 
        
           close 
        
           % Identify a system model 
        
           k = 50; 
        
           sampleTime = t(2) - t(1); 
        
           delay = 0; 
        
           systemorder = 2; 
        
           Ghat = mi.cca(u, y, k, sampleTime, delay, systemorder); 
        
           % Find the disturbance d = H*e 
        
           Ad = Ghat.A; 
        
           Bd = Ghat.B; 
        
           Cd = Ghat.C; 
        
           Dd = Ghat.D; 
        
           x = zeros(systemorder, 1); 
        
           for i = 1:size(t, 2) 
        
             yhat(:,i) = Cd*x + Dd*u(:,i); 
        
             x = Ad*x + Bd*u(:,i); % Update state vector 
        
           end 
        
           d = y - yhat; 
        
           % Identify the disturbance model 
        
           systemorder = 2; 
        
           ktune = 0.5; 
        
           [Hhat] = mi.sra(d, k, sampleTime, ktune, delay, systemorder); 
        
           % Simulate the disturbance model 
        
           [dy, dt] = mc.lsim(Hhat, e, t); 
        
           close 
        
           plot(dt, dy, t, d); 
        
           legend('d = Hhat*e(t)', 'd = y - yhat') 
        
           grid on

Example SRA 2

MataveID/examples/sraExample2.m

Lines 1 to 35 in 2014b74

    
           % Clear all 
        
           clear all 
        
           % Create disturbance signal 
        
           t = linspace(0, 100, 1000); 
        
           e = randn(1, length(t)); 
        
           % Create disturbance model 
        
           H = mc.tf([1], [1 3]); 
        
           % Simulate 
        
           y = mc.lsim(H, e, t); 
        
           close 
        
           % Identify a model 
        
           k = 100; 
        
           sampleTime = t(2) - t(1); 
        
           ktune = 0.01; 
        
           delay = 0; 
        
           systemorder = 2; 
        
           [H, K] = mi.sra(y, k, sampleTime, ktune, delay, systemorder); 
        
           % Observer 
        
           H.A = H.A - K*H.C; 
        
           % Create new signals 
        
           [y, t] = mc.gensig('square', 10, 10, 100); 
        
           y = [y*5 y*2 -y 10*y -2*y]; 
        
           t = linspace(0, 100, length(y)); 
        
           % Add some noise 
        
           y = y + 2*randn(1, length(y)); 
        
           % Simulate 
        
           mc.lsim(H, y, t);

Local Binary Pattern

Use LBP if you want to find a binary pattern inside of a matrix, or an image X around a pixel P = X(y, x)

[descriptor] = mi.lbp(X, x, y, radius, init_angle, lbp_bit);

Example

MataveID/examples/lbpExample.m

Lines 1 to 21 in 2e81144

    
           % Close all 
        
           close all 
        
           % Read image 
        
           X = imread(fullfile('..', 'data', 'lab.pgm')); 
        
           % Make image greyscale 
        
           if size(X, 3) > 1 
        
             X =  rgb2gray(X); 
        
           else 
        
             X = double(X); 
        
           end 
        
           % Compute Local Binary Pattern 
        
           radius = 100; 
        
           init_angle = 0; 
        
           lbp_bit = 32; 
        
           x = 110; 
        
           y = 110; 
        
           descriptor = mi.lbp(X, x, y, radius, deg2rad(init_angle), lbp_bit); 
        
           fprintf('%i = 0b%s\n', descriptor, dec2bin(descriptor))

Result

288882440 = 0b10001001101111111111100001000

Linear Discriminant Analysis

Linear Discriminant Analysis can be used for dimension reduction and projection on maximum distance between classes.

[P, W] = mi.lda(X, y, c);

Linear Discriminant Analysis example

MataveID/examples/ldaExample.m

Lines 1 to 70 in 74506ea

    
           % How much data 
        
           l = 50; 
        
           % Data for the first class 
        
           x1 = 2*randn(1, l); 
        
           y1 = 50 + 5*randn(1, l); 
        
           z1 = 1:l; 
        
           % Data for the second class 
        
           x2 = 5*randn(1, l); 
        
           y2 = -4 + 2*randn(1, l); 
        
           z2 = 2*l:-1:l+1; 
        
           % Data for the third class 
        
           x3 = 15 + 3*randn(1, l); 
        
           y3 = 50 + 2*randn(1, l); 
        
           z3 = -l:-1; 
        
           % Create the data matrix in this way 
        
           A = [x1 x2 x3]; 
        
           B = [y1 y2 y3]; 
        
           C = [z1 z2 z3]; 
        
           X = [A; B; C]; 
        
           % Create class ID, indexing from zero 
        
           y = [1*ones(1, l), 2*ones(1, l), 3*ones(1, l)]; 
        
           % How many dimensions 
        
           c = 2; 
        
           % Plot original data 
        
           close all 
        
           scatter3(x1, y1, z1, 'r') 
        
           hold on 
        
           scatter3(x2, y2, z2, 'b') 
        
           hold on 
        
           scatter3(x3, y3, z3, 'g') 
        
           title('Original 3D data', 'FontSize', 20) 
        
           legend('Class 1', 'Class 2', 'Class 3') 
        
           % Do LDA for 2D 
        
           [P, W] = mi.lda(X, y, c); 
        
           % Plot 2D were P is a c x l matrix 
        
           figure 
        
           scatter(P(1, 1:l), P(2,1:l), 'r') 
        
           hold on 
        
           scatter(P(1,l+1:2*l), P(2, l+1:2*l), 'b') 
        
           hold on 
        
           scatter(P(1, 2*l+1:3*l), P(2, 2*l+1:3*l), 'g') 
        
           grid on 
        
           title('Dimension reduction 2D data', 'FontSize', 20) 
        
           legend('Class 1', 'Class 2', 'Class 3') 
        
           % How many dimensions 
        
           c = 1; 
        
           % Do LDA for 1D 
        
           [P, W] = mi.lda(X, y, c); 
        
           % Plot 1D were P is a c x l matrix 
        
           figure 
        
           scatter(P(1, 1:l), 0*P(1, 1:l), 'r') 
        
           hold on 
        
           scatter(P(1, l+1:2*l), 0*P(1, l+1:2*l), 'b') 
        
           hold on 
        
           scatter(P(1, 2*l+1:3*l), 0*P(1, 2*l+1:3*l), 'g') 
        
           grid on 
        
           title('Dimension reduction 1D data', 'FontSize', 20) 
        
           legend('Class 1', 'Class 2', 'Class 3')

Logistic regression

Use this if you have a binary output (0, 1) or (-1, 1) and you want to to have a probabilistic output e.g 0 to 100%

[a, b, flag, iterations] = mi.logreg(x, y, function_type)

Example logreg

MataveID/examples/logregExample.m

Lines 1 to 29 in 351d612

    
           % Clear 
        
           close all 
        
           clear all 
        
           clc 
        
           % Logistic data 
        
           x = [-0.1, -0.2, -0.3, -0.311, -0.213, -0.133, -0.231, -0.4215, 0.13, 0.23, 0.19, 0.9, 1.2, 1.5, 0.423, 0.561]; 
        
           % Create the parameters a and b for tanh 
        
           y = [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1]; 
        
           [a, b, flag, iterations] = mi.logreg(x, y, 'tanh'); 
        
           % Plot the logistic function: tanh 
        
           t = linspace(-10, 10, 1000); 
        
           p = (exp(a*t + b) - exp(-a*t - b))./(exp(a*t + b) + exp(-a*t - b)); 
        
           plot(t, p) 
        
           title('tanh function', 'FontSize', 20) 
        
           grid on 
        
           % Create the parameters a and b for sigmoid 
        
           y = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1]; 
        
           [a, b, flag, iterations] = mi.logreg(x, y, 'sigmoid'); 
        
           % Plot the logistic function: sigmoid 
        
           t = linspace(min(x), max(x), 1000); 
        
           p = 1./(1 + exp(-(a*t + b))); 
        
           figure 
        
           plot(t, p) 
        
           title('sigmoid function', 'FontSize', 20)

Linear Support Vector Machine

This is the standard way to create a support vector machine. Even if it's only returning back a linear model, it's still very powerful and suits systems that need extreamly fast predictions such as embedded systems.

Notice that the Linear Support Vector Machine can only do two-class prediction only. But you can use multiple classes with the Linear Support Vector Machine by using multiple linear support vector machines. It's called One-VS-All method.

[w, b, accuracy, solution] = mi.lsvm(x, y, C, lambda)

Linear Support Vector Machine 2D example

MataveID/examples/lsvm2DExample.m

Lines 1 to 59 in 817443d

    
           % Data 
        
           X = [5 3; 
        
                2 1; 
        
                7 2; 
        
                8 3; 
        
                9 1; 
        
                15 23; 
        
                17 18; 
        
                18 13; 
        
                16 20; 
        
                19, 15]; 
        
           % Labels of the data for each class 
        
           y = [1; 
        
                1; 
        
                1; 
        
                1; 
        
                1; 
        
                -1; 
        
                -1; 
        
                -1; 
        
                -1; 
        
                -1]; 
        
           % Plot 2D 
        
           scatter(X(y == -1,1),X(y == -1,2), 'r'); 
        
           hold on 
        
           scatter(X(y == 1,1), X(y == 1,2), 'g'); 
        
           grid on 
        
           legend('Class A', 'Class B', 'location', 'northwest') 
        
           % Tuning parameters 
        
           C = 1; % For upper boundary limit 
        
           lambda = 1; % Regularization (Makes it faster to solve the quadratic programming) 
        
           % Compute weigths, bias and find accuracy 
        
           [w, b, accuracy, solution] = mi.lsvm(X, y, C, lambda); 
        
           % How long the line should be 
        
           min_value_column_1 = min(X(:,1)); 
        
           max_value_column_1 = max(X(:,1)); 
        
           % Create the separation line y = k*x + m 
        
           x1 = linspace(min_value_column_1, max_value_column_1); 
        
           x2 = (-w(1)*x1 - b) / w(2); 
        
           % Plot the separation line 
        
           plot(x1, x2, 'k', 'LineWidth', 2); 
        
           xlim([0 20]) % Max x-axis limit 
        
           ylim([0 20]) % Max y-axis limit 
        
           legend('Class A', 'Class B', 'Separation', 'location', 'northwest') 
        
           % Classify 
        
           x_unknown = [15; 5]; 
        
           class_ID = sign(w*x_unknown + b) 
        
           if(class_ID > 0) 
        
             disp('x_unknown class B') 
        
           else 
        
             disp('x_unknown is class A') 
        
           end

Linear Support Vector Machine 3D example

MataveID/examples/lsvm3DExample.m

Lines 1 to 57 in 817443d

    
           % Data 
        
           X = [5 3 2; 
        
                2 1 3; 
        
                7 2 4; 
        
                8 3 1; 
        
                9 1 2; 
        
                15 23 23; 
        
                17 18 13; 
        
                18 13 63; 
        
                16 20 24; 
        
                19, 15 52]; 
        
           % Labels of the data for each class 
        
           y = [1; 
        
                1; 
        
                1; 
        
                1; 
        
                1; 
        
                -1; 
        
                -1; 
        
                -1; 
        
                -1; 
        
                -1]; 
        
           % Plot 3D 
        
           scatter3(X(y == -1,1),X(y == -1,2), X(y == -1,3), 'r'); 
        
           hold on 
        
           scatter3(X(y == 1,1), X(y == 1,2), X(y == 1,3), 'g'); 
        
           grid on 
        
           legend('Class A', 'Class B', 'location', 'northwest') 
        
           % Tuning parameters 
        
           C = 1;              % For upper boundary limit 
        
           lambda = 1;         % Regularization (Makes it faster to solve the quadratic programming) 
        
           % Compute weigths, bias and find accuracy 
        
           [w, b, accuracy, solution] = mi.lsvm(X, y, C, lambda); 
        
           % Definiera området för 3D-plot 
        
           x1Range = linspace(min(X(:,1))-1, max(X(:,1))+1, 50); 
        
           x2Range = linspace(min(X(:,2))-1, max(X(:,2))+1, 50); 
        
           [x1Grid, x2Grid] = meshgrid(x1Range, x2Range); 
        
           x3Grid = (-w(1)*x1Grid - w(2)*x2Grid - b) / w(3); 
        
           % Plot the hyperplane 
        
           surf(x1Grid, x2Grid, x3Grid, 'FaceAlpha', 0.5); 
        
           colormap(gray); 
        
           legend('Class A', 'Class B', 'Separation', 'location', 'northwest') 
        
           % Classify 
        
           x_unknown = [15; 5; 7]; 
        
           class_ID = sign(w*x_unknown + b) 
        
           if(class_ID > 0) 
        
             disp('x_unknown class B') 
        
           else 
        
             disp('x_unknown is class A') 
        
           end

Linear Support Vector Machine 9D example

This example demonstrates how to use more than 3 columns in SVM. Notice that here we don't plot this 9D measurements

MataveID/examples/lsvmXDExample.m

Lines 1 to 45 in b35c7ca

    
           % This is an example for more than 3 columns of X. 
        
           % Here we cannot plot this 9D table because nobody knows how to plot 9D plots. 
        
           % Data - Each column inside the table is a sensor 
        
           X = [% Class B table 
        
                5 3 2 2 4 5 6 7 3; 
        
                2 1 3 3 3 4 3 2 5; 
        
                7 2 4 3 1 9 4 2 4; 
        
                8 3 1 1 5 7 8 2 9; 
        
                9 1 2 2 3 1 5 3 2; 
        
                % Class A table 
        
                15 23 23 32 43 52 13 64 34; 
        
                17 18 13 34 54 10 45 99 77; 
        
                18 13 63 56 33 95 35 65 55; 
        
                16 20 24 93 94 23 56 87 77; 
        
                19 15 52 36 20 45 44 22 32]; 
        
           % Labels of the data for each class 
        
           y = [1; % Class B 
        
                1; 
        
                1; 
        
                1; 
        
                1; 
        
                -1; % Class A 
        
                -1; 
        
                -1; 
        
                -1; 
        
                -1]; 
        
           % Tuning parameters 
        
           C = 1;              % For upper boundary limit 
        
           lambda = 1;         % Regularization (Makes it faster to solve the quadratic programming) 
        
           % Compute weigths, bias and find accuracy 
        
           [w, b, accuracy, solution] = mi.lsvm(X, y, C, lambda); 
        
           % Classify 
        
           x_unknown = [15; 5; 7; 2; 4; 6; 3; 5; 2]; 
        
           class_ID = sign(w*x_unknown + b) 
        
           if(class_ID > 0) 
        
             disp('x_unknown class B') 
        
           else 
        
             disp('x_unknown is class A') 
        
           end

Nonlinear Support Vector Machine with C code generation

This algorithm can do C code generation for nonlinear models. It's a very simple algorithm because the user set out the support points by using the mouse pointer. When all the supports are set ut, then the algorithm will generate C code for you so you can apply the SVM model in pure C code using CControl library.

All you need to have is two matrices, X and Y. Where the column length is the data and the row length is the amount of classes. The nlsvm.m file will plot your data and then when you have placed out your support points, then the svm.m will generate C code for you that contains all the support points.

If you have let's say more than two variables, e.g Z matrix or even more. Then you can create multiple models as well by just using diffrent data as arguments for the svm function below. The C code generation is very fast and it's very easy to build a model.

[X_point, Y_point, amount_of_supports_for_class] = mi.nlsvm(X, Y)

Nonlinear Support Vector Machine with C code generation example

MataveID/examples/nlsvmExample.m

Lines 1 to 54 in 817443d

    
           clc; clear; close all; 
        
           % How much data should we generate  
        
           N = 50; 
        
           % How many classes  
        
           c = 5; 
        
           % Create variance and average for X and Y data 
        
           X_variance = [2, 4, 3, 4, 5]; 
        
           Y_variance = [3, 5, 3, 4, 5]; 
        
           X_average = [50, 70, 10, 90, 20]; 
        
           Y_average = [20, 70, 60, 10, 20]; 
        
           % Create scatter data 
        
           X = zeros(c, N); 
        
           Y = zeros(c, N); 
        
           for i = 1:c 
        
             % Create data for X-axis  
        
             X(i, 1:N) = X_average(i) + X_variance(i)*randn(1, N); 
        
             % Create data for Y-axis 
        
             Y(i, 1:N) = Y_average(i) + Y_variance(i)*randn(1, N); 
        
           end 
        
           % Create SVM model - X_point and Y_point is coordinates for the Nonlinear SVM points. 
        
           % amount_of_supports_for_class is how many points there are in each row 
        
           [X_point, Y_point, amount_of_supports_for_class] = mi.nlsvm(X, Y); 
        
           % Do a quick re-sampling of random data again 
        
           for i = 1:c 
        
             % Create data for X-axis  
        
             X(i, 1:N) = X_average(i) + X_variance(i)*randn(1, N); 
        
             % Create data for Y-axis 
        
             Y(i, 1:N) = Y_average(i) + Y_variance(i)*randn(1, N); 
        
           end 
        
           % Check the SVM model 
        
           point_counter_list = zeros(1, c); 
        
           for i = 1:c 
        
             % Get the points  
        
             svm_points_X = X_point(i, 1:amount_of_supports_for_class(i)); 
        
             svm_points_Y = Y_point(i, 1:amount_of_supports_for_class(i)); 
        
             % Count how many data points this got - Use inpolygon function that return 1 or 0 back 
        
             point_counter_list(i) = sum(inpolygon(X(i,:) , Y(i, :), svm_points_X, svm_points_Y)); 
        
           end 
        
           % Plot how many each class got - Maximum N points per each class 
        
           figure  
        
           bar(point_counter_list); 
        
           xlabel('Class index'); 
        
           ylabel('Points');

Here is an application with SVM for a hydraulical system. This little box explains whats happening inside the hydraulical system if something happen e.g a motor or a valve is active. It can identify the state of the system.

NN - Neural Network

This generates a neural network back and an activation function. This Neural Network is tranied by Support Vector Machine.

[weight, bias, activation_function] = mi.nn(data, class_id, C, lambda);

Example

Here I'm using Fisher's Irish dataset to train a neural network.

MataveID/examples/nnExample.m

Lines 1 to 29 in 30c4904

    
           % Clear 
        
           close all 
        
           clear all 
        
           clc 
        
           % Data - Avoid the header and the last column 
        
           data_raw = double(csvread('..','data','iris.csv')(2:end, 1:end-1)); 
        
           % 3 class labels 
        
           class_id = [linspace(1, 1, 50)'; linspace(2, 2, 50)'; linspace(3, 3, 50)']; 
        
           % Create neural network 
        
           C = 2; 
        
           lambda = 0.5; 
        
           [weight, bias, activation_function] = mi.nn(data_raw, class_id, C, lambda); 
        
           % Check accuracy 
        
           X = weight*data_raw' + bias; 
        
           classes = length(class_id); 
        
           score = 0; 
        
           for i = 1:classes 
        
             class_id_predicted = activation_function(X(:, i)); 
        
             if(class_id_predicted == class_id(i)) 
        
               score = score + 1; 
        
             end 
        
           end 
        
           % Print status 
        
           fprintf('The accuracy of this model is: %i\n', score/classes*100);

Output:

Training: Neural Network success with accuracy: 1.000000 at class: 1
Training: Neural Network success with accuracy: 0.733333 at class: 2
Training: Neural Network success with accuracy: 0.986667 at class: 3
The accuracy of this model is: 96.6667

Robust Principal Component Analysis

Robust principal component analysis(RPCA) is a great tool if you want to separate noise from data X into a matrix S. RPCA is a better tool than PCA because it using optimization and not only reconstructing the image using SVD, which PCA only does.

[L, S] = rpca(X);

Robust Principal Component Analysis example

MataveID/examples/rpcaExample.m

Lines 1 to 13 in 2014b74

    
           clc; clear; close all; 
        
           file = fullfile('..','pictures','bob.png'); 
        
           X = imread(file); % Load Mr Bob 
        
           X = rgb2gray(X);       % Grayscale 8 bit 
        
           X = double(X);         % Must be double 40 => 40.0 
        
           [L, S] = mi.rpca(X);      % Start RPCA. Our goal is to get L matrix 
        
           figure(1) 
        
           imshow(uint8(X))       % Before RPCA 
        
           title('Before RPCA - Bob') 
        
           figure(2) 
        
           imshow(uint8(L))       % After RPCA 
        
           title('After RPCA - Bob')

Install

To install MataveID, download the folder "matave" and place it where you want it. Then the following code need to be written inside of the terminal of your MATLAB® or GNU Octave program.

path('path/to/the/folder/matave', path)
savepath

Example of a typical path.

path('C:\Users\dmn\Documents\Octave\matave\', path)
savepath

Package requriments:

MataveControl package.

Update

Write this inside the terminal. Then MataveID is going to download new .m files to MataveID from GitHub

mi.updatemataveid

License

DanielMartensson/MataveID

About

Topics

Resources

License

Stars

Watchers

Forks

Releases 3

Packages 0

Contributors 2

Languages

Packages