Richard Craig - Research Engineer
 
 

Archive for Matlab


Matlab Image Correlation

Visually Evoked Potentials (VEPs) are small changes in the brain signal cause by a visual stimulus. Responses are recorded using electroencephalogram (EEG) electrodes located near the occipital cortex, the area of the brain involved in receiving and interpreting visual signals. I used the following matlab code to check the correlation of images that were going to be used as stimuli for the EEG experiment.

I used the Corr2 function in Matlab. Corr2(A,B) computes the correlation coefficient between A and B, where A and B are matrices or vectors of the same size. Corr2 computes the correlation coefficient using the equation;

In case you want to try and type that into Latex
r = \frac{\sum_{m} \sum_{n}(A_{mn} - \bar{A}) (B_{mn} - \bar{B})}{\sqrt{\left ( \sum_{m} \sum_{n}(A_{mn} - \bar{A})^{2} \right ) \left ( \sum_{m} \sum_{n}(B_{mn} - \bar{B})^{2} \right )}}

Theres a handy Latex visual equation editor on the codecogs.com website.

Matlab Code

root ='C:\Users\Richard\Documents\My Dropbox\EXPERIMENTS\images\';
temp = (rgb2gray(imread([root 'FACE_v201.JPG'] ))); % Template Image to correlate against
for x=1:100
cc(x) = corr2(temp,(rgb2gray(imread([root 'FACE_v20' num2str(x) '.JPG'] )))); % correlate images
end

Matlab gives this example
I = imread(‘pout.tif’);
J = medfilt2(I);
R = corr2(I,J)
R =
0.9959

Create a Latin square in MATLAB

a Latin square is an n × n array filled with n different Latin letters, each occurring exactly once in each row and exactly once in each column.

Latin squares are used in statistics and in mathematics.

% Generate a Latin Square
% =============================
% for Matlab R13
% version 1.0 (feb 2010)
%
% Inputs
% ——————
% N – Integer value for the size of the square.
%
% Inputs
% ——————
% M – Your Latin Square
%
% Other functions
% ——————-
% RANDPERM(n) is a random permutation of the integers from 1 to n.
% For example, RANDPERM(6) might be [2 4 5 6 1 3].
%
% CUMSUM(X) is a vector containing the cumulative sum of the elements of X.
% Example: If X = [0 1 2
% 3 4 5]
% then cumsum(X,1) is [0 1 2 and cumsum(X,2) is [0 1 3
% 3 5 7] 3 7 12]
% email: richard [ a t ] richard-craig.co.uk

% Download if too lazy to copy and paste

function [M] = latsq(N)

M = [randperm(N); ones(N-1,N)] ; % Generate Square
M = rem(cumsum(M)-1,N) + 1 ; % Add one due to base zero

Update Simulink variables from MATLAB GUI

This example is designed to demonstrate how to control variables within a SIMULINK model before starting the simulation or in real time during the simulation. Although the same commands can be used from the MATLAB command window, I wanted to use a Graphical User Interface (GUIs) as it would co-ordinate with other applications.

Because the GUI sets parameters and runs the simulation, the 'EEG_BERT_FACE.mdl' model must be open when the GUI is displayed.
The purpose of the subfunction is to:

  • Determine if the model is open (find_system).
  • Open the block diagram for the model and the subsystem where the parameters are being set, if not open already (open_system).
  • Change the size of the controller Gain block so it can display the gain value (set_param).
  • Bring the GUI forward so it is displayed on top of the Simulink diagrams (figure).
  • Set the block parameters to match the current settings in the GUI.

Here is the code for the model_open subfunction:

function SIMULINK_model_open(handles)if isempty(find_system(‘Name’,'EEG_BERT_FACE’)),

% Open the model file ‘EEG_BERT_FACE.mdl’
open_system(‘EEG_BERT_FACE’);

% Variables updated from the MATLAB GUI
set_param(‘EEG_BERT_FACE/To File’,'Filename’,handles.sfilename); %Block Property Variable

%Gain Parameters
set_param(‘EEG_BERT_FACE/gStimulus’,'Gain’,'0′);
set_param(‘EEG_BERT_FACE/gFaceID’,'Gain’,'0′);
set_param(‘EEG_BERT_FACE/gFaceCount’,'Gain’,'0′);
set_param(‘EEG_BERT_FACE/gStepNumber’,'Gain’,'0′);
set_param(‘EEG_BERT_FACE/gNumberOfRuns’,'Gain’,'1′);
set_param(‘EEG_BERT_FACE/gVarDelay’,'Gain’,num2str(handles.time_addition));

end

function SIMULINK_model_update(handles)
%Gain Parameters
set_param(‘EEG_BERT_FACE/gStimulus’,'Gain’,'0′);
set_param(‘EEG_BERT_FACE/gFaceID’,'Gain’,num2str(handles.fID));
set_param(‘EEG_BERT_FACE/gFaceCount’,'Gain’,num2str(handles.iFace));
set_param(‘EEG_BERT_FACE/gStepNumber’,'Gain’,num2str(handles.iStep));
set_param(‘EEG_BERT_FACE/gNumberOfRuns’,'Gain’,num2str(handles.iRuns));
set_param(‘EEG_BERT_FACE/gVarDelay’,'Gain’,num2str(handles.time_addition));

function SIMULINK_model_close(handles)

set_param(‘EEG_BERT_FACE’,'SimulationCommand’,'pause’);
set_param(‘EEG_BERT_FACE’,'SimulationCommand’,’stop’);
save_system(‘EEG_BERT_FACE’);
close_system(‘EEG_BERT_FACE’);

Matlab Quadratic Programming Optimisation Problem

Problem

Find minimum of   x^2 + 4y^2 -3xy -3x -7y

Subject to   0<x<7    0<y<7

Change constraints till you hit both

Matlab Solution

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Programme  : EngD Maths for Systems
% Date       : 2011-06-01
% Session    : Optimisation
% Exercise   : Material Blending
% By         : Richard Craig
% % Description: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% quadprog4.m
%Find minimum of  x^2 + 4y^2 -3xy -3x -7y
%Subject to   0<x<7    0<y<7

clc; clear;
%% Set up function for optimization
% H will give expression of form:
%     h11.x1^2 +h12.x1.x2 +h21.x2.x1 +h22.x2^2
% but programme forces symmetric Hessian: h12=h21 etc.,
% so might as well make it symmetric now.

% from the eq we can see H = [x^2 + 4y^2 -3xy] f = [ -3x -7y]
H = [2 -3; -3 8];       %H = [x^2 + 4y^2 -3xy]

% f will give inner product f.x
f = [-3; -7];           % f = [ -3x -7y]
%% Constraints
%x = quadprog(H,f,A,b,Aeq,beq) solves the preceding problem while
%additionally satisfying the equality constraints Aeq*x = beq.
%If no inequalities exist, set A = [] and b = [].
A = [];b = [];
%lower bounds on x:
lb = zeros(2,1);
%upper bound on x:
ub = [16 7];    %[7 7];     %Subject to   0<x<7    0<y<7
%% View Function
% make a "meshgrid" to cover allowed ranges of the variables.
u=linspace(lb(1),ub(1),11); v=linspace(lb(2),ub(2),11);
[X, Y]=meshgrid(u,v);
% put a value to each point-pair of the grid:
Z=0.5*(H(1,1)*X.^2 + H(2,2)*Y.^2 + 2*H(1,2).*X.*Y) + f(1).*X + f(2).*Y;
%plot the surface with some colour and a contour plot underneath:
h = figure(1);
surfc(X,Y,Z)
saveas(h,'QuadSurface.jpg');
%% Minimize using 'quadprog'
% Call for quadratic programming method:
[x,fval] = quadprog(H,f,A,b,[],[],lb,ub)

MATLAB updates Access database via SQL

Matlab is able to store variables in files and can import external information from excel and access, but sometimes you want to keep the information in Access and only do some processing in MATLAB.  Here is some code to use Structured Query Language (SQL) to update an Access database from within MATLAB;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% File       : DatabaseConnection.m
% Date       : 2011-07-04
% Revision   : 1.0
% By         : Richard Craig
%
% Description:
% This file imports external information from an Access database into
% Matlab.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Clear all variables, close all windows
clear all;close all;clc;

%% % Taken from http://www.mathworks.com/help/toolbox/database/ug/database.html
% Specify the location of the database on disk.
dbpath = ['C:\Users\Richard\Desktop\Database.mdb'];

% Create the connection URL.
conurl = ['jdbc:odbc:Driver={Microsoft Access Driver (*.mdb)};' ...
   'DBQ=' dbpath];

% Connect to the database.
conn = database('','','','sun.jdbc.odbc.JdbcOdbcDriver', conurl);

if isconnection(conn)
    display('Connected');

    % Fetch data from the database.
    %Import table as-is
    %e = exec(conn,'SELECT * FROM tbl_table'); %Run an SQL statement
    % Apply a filter to results
    e = exec(conn,'SELECT * FROM tbl_table WHERE [Type] like ''Academic (Civil)'' '); %Run an SQL statement

    % Run SQL statement
    e = fetch(e);
    data = e.Data;          % Transfer data to workspace variable
    display('Information imported');

    [rows cols] = size(data);

    for i=1:rows
        % Informaton will be stored as a cell so will need converstions to
        % string (char(var)) and then numeric (str2num(char(var))
        % Assign to variables
        ID = cell2mat(data(i,1));%Numeric
        Name = char(data(i,2));%String
        Url = char(data(i,3));%String

        Url = strrep(Url,'#',''); %Remove errors in string
        data(i,3) = cellstr(Url);% Reassign string to cell

        % but that only updates the variable not the actual table in our
        % database.
        % We need to use an SQL statement;
        s =['UPDATE [URL] FROM tbl_table WHERE [ID] =' num2str(ID)];
        %sql = exec(conn,s);%execute statement
        %BE SURE the sql statement works as you want on the local data
        %variable and BACK UP your information before running ANY SQL
        %statement for the first time.

        %Display Result
        disp([num2str(ID) ' - ' Name ' ' Url ]);

    end%i=1:rows

end %isconnection(conn)

%% Close the connection.
close(conn);

Matlab Material Blending Optimisation Problem

Problem

Various source alloys, A,B,C,D,E are available with different proportions of elemental metals and costs. Need to make final blend with 30% Lead, 20% Zinc and 50% Tin at lowest cost.

A B C D E Req.
Lead 10 10 40 60 30 30
Zinc 10 30 50 30 30 20
Tin 80 60 10 10 40 50
Cost 8. 9. 11. 13 17

Adapted from Ref 2, Bhatti, p424

Matlab code

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Programme  : EngD Maths for Systems
% Date       : 2011-06-01
% Session    : Optimisation
% Exercise   : Material Blending
% By         : Richard Craig
% Description:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Clear previous variables
clc; clear all;

%% Blending Problem.
%         A        B        C        D        E        Req
% Lead    10        10        40        60        30        30
% Zinc    10        30        50        30        30        20
% Tin     80        60        10        10        40        50
% Cost    8.0        9.0        11.0    13.0    17.0
%
%
% Various source alloys, A,B,C,D,E are available with different proportions
% of elemental metals and costs.
M =  [10   10    40    60    30
10   30    50    30    30
80   60    10    10    40];

% Final solution material boundary
b = [30; 20; 50];

% Define a lower bound and an upper bound for each variable.
lb = zeros(5,1);    % Can't have minus material
ub = ones(5,1);     % Can't have more than 100% of one material 

% define the objective function for Minimize(f.x)
f = [8.0;  9.0; 11.0;  13.0; 17.0];

%Price per unit %% call the linear programming operation

[x,fval,exitflag,output,lambda] = linprog(f,[],[],M,b,lb,ub);
% lambda.ineqlin
% lambda.eqlin
% lambda.upper
% lambda.lower

%% ANSWER
x

% Final Blend of materials id given in the proportions;
% x = 0.5714
%     0.0000
%     0.0714
%     0.3571
%     0.0000

Matlab Open Box Optimisation Problem

Problem

Builder wants open boxes, carry 1200 cu ft of stuff;  walls, ends, base different construction, different (lifetime) costs. Wants to minimize (lifetime) cost.

Problem is Min{48bh+30bl+44hl} subject to bhl=1200.

Find: b=breadth, l=length, h=height.

Matlab Solution

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Programme  : EngD Maths for Systems
% Date       : 2011-06-01
% Session    : Optimisation
% Exercise   : OPEN BOX PROBLEM
% By         : Richard Craig
%
% Description:
% Builder wants open boxes, carry 1200 cu ft of stuff;
% walls, ends, base different construction, different (lifetime) costs.
% Wants to minimize (lifetime) cost.
% Problem is Min{48bh+30bl+44hl} subject to the constraint that bhl=1200.
% Find: b=breadth, l=length, h=height.
%
% Adapted from Ref 2, Bhatti, “Practical Optimization Methods”.
h = 12; % Height
b = 12; % Breadth
l = 12; % Length
x0 = [h; b; l]

%Cost function (Min{48bh+30bl+44hl}):
fun = @(x) (48*x(1)*x(2))+(30*x(2)*x(3))+(44*x(1)*x(3));

%Constraint
C = 1200;
lb =zeros(1,3); ub = [C,C,C];

%Constrain volume
%By using @nonlcon4 function in separate m file called nlcon4.m:
%       function [c,ceq] = nonlcon4(x)
%       %Input variable x = [H,B,L]
%           ceq = 1200-x(1)*x(2)*x(3);
%           c = [];
%       end

%Minimization variables A,b,Aeq,beq,

%% Find minimum of constrained nonlinear multivariable function
%x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,@nonlcon4)
x = fmincon(fun,x0,[],[],[],[],lb,ub,@nonlcon4)

% x =
%     7.9967
%    11.7285
%    12.7947

y = fun(x) % y = 1.3506e+004

%% Extra
% Vary capacity in nonlcon4 to 1210
C = 1210;
ub = [C,C,C]; % new upper limits
x = fmincon(fun,x0,[],[],[],[],lb,ub,@nonlcon4) %Re-eval
% x =
%
%     8.0188
%    11.7610
%    12.8301

y = fun(x) % y = 1.3581e+004

.