睡眠阶段分期——SVM和ELM分别与粒子群算法结合（function）

importslpdb函数

loadmatobject函数

extractfeatures函数

getindexrange函数

PSOforELM函数

PSOforSVM函数

extractresults函数

importslpdb函数

——将头文件、RR文件、注释文件构成结构体：OutputData = struct('filename', fileName, 'time', anTimeGeneratedCell, ...'rr', rrCollection, 'annotation', anClassGeneratedCell, 'age', age, ...'gender', gender, 'weight', weight);

function OutputData = importslpdb(fileName)
%Import and synchronize a slpdb recording
%   Syntax:
%   OutputData = importslpdb(fileName)
%
%   Input:
%   *) fileName   - slpdb file name to be imported. Example: 'slp01a'.
%                   file must be located in 'slpdb' folder,
%                   three file formats needed: .hea, .rr, and .an
%
%   Output:
%   *) OutputData - struct contains synchronized RR and annotationfileName = strcat('slpdb/', cell2mat(fileName));
SEC_PER_EPOCH = 30; % amount of seconds in one epoch (value for slpdb is 30)
OutputData = [];
fprintf('\n%s DATA IMPORT...\n', fileName);

%% IMPORT HEADER DATA   导入头文件数据fprintf('Importing header file...\n');
fid = fopen(strcat(fileName, '.hea'), 'r');
headerFile = textscan(fid, '%s%s%s%s%s%s%s%*[^\n]');
fclose(fid);
heaSamplingFreq = strsplit(char(headerFile{3}(1)), '/');
heaSamplingFreq = str2double(cell2mat(heaSamplingFreq(1)));
heaTotalSamples = str2double(cell2mat(headerFile{4}(1)));
heaRecLengthInSec = ceil(heaTotalSamples/heaSamplingFreq);
heaTotalEpoch = ceil(heaRecLengthInSec/SEC_PER_EPOCH);
heaIdx = size(headerFile{1}, 1); % get last line index of fileif cell2mat(headerFile{1}(end-1)) == '#'% decease index by 1 for 'slp37.hea',% because the last line is not age, gender, and weight information    heaIdx = heaIdx - 1;
endage = headerFile{2}(heaIdx);
gender = headerFile{3}(heaIdx);
weight = headerFile{4}(heaIdx);% output of "IMPORT HEADER DATA" section:
% *) heaTotalEpoch - total epoch according to header data
% *) age           - age of the subject
% *) gender        - gender of the subject
% *) weight        - weight of the subject
% END OF IMPORT HEADER DATA

%% IMPORT ANNOTATION DATA     导入注释数据fprintf('Importing annotation file...\n');
fid = fopen(strcat(fileName, '.an'), 'r');
anFile = textscan(fid, '%s%s%s%s%s%s%s%*[^\n]');
fclose(fid);% remove header of annotation data (first line)
% remove this string: 'Elapsed time  Sample #  Type  Sub Chan  Num Aux'for i=1:size(anFile, 2)anFile{i}(1) = [];
end% change first epoch's 'start time' into 0:00.000
anTemp = cell2mat(anFile{1}(1));
anTemp(end-2:end) = 48;
anFile{1}(1) = cellstr(anTemp);
anTime = anFile{1};
anClass = anFile{7};% output of "IMPORT ANNOTATION DATA" section:
% *) anTime  - time from annotation file (cell array)
% *) anClass - annotation (cell array)
% END OF IMPORT ANNOTATION DATA

%% IMPORT RR DATA      导入RR数据fprintf('Importing RR file...\n');
fid = fopen(strcat(fileName, '.rr'), 'r');
rrFile = textscan(fid, '%s%s%s%s%s%*[^\n]');
fclose(fid);rrConvertedTime = rrFile{1};for i=1:size(rrConvertedTime, 1);rrStartTimeChar = cell2mat(rrConvertedTime(i));     % convert cell into charrrStartTimeChar(end-2:end) = 48;    % xx:xx:xx.aaa -> change 'aaa' part to '000'% split start time by ":" into matrix    求时长rrStartTimeMat = strsplit(char(rrConvertedTime(i)), ':')';    % get seconds from the last element      rrSecond = str2double(cell2mat(rrStartTimeMat(end)));% epoch groupingrrWhichGroup = floor(rrSecond/SEC_PER_EPOCH)*SEC_PER_EPOCH;% set epoch groupingif rrWhichGroup == 0rrStartTimeChar(end-5) = 48;elseif rrWhichGroup == 30rrStartTimeChar(end-5) = 51;endrrStartTimeChar(end-4) = 48;rrConvertedTime(i) = mat2cell(rrStartTimeChar, 1);end% change RR value from 'array of cell' into 'array of double'  改变RR值属性rrNum = zeros(size(rrFile{3}, 1), 1);for i=1:size(rrFile{3}, 1)rrNum(i) = str2double(cell2mat(rrFile{3}(i)));
endrrTime = rrFile{1};% output of "IMPORT RR DATA" section:
% *) rrConvertedTime - rounded RR start time according to the epoch
%                      example: 1:34:31.328 -> 1:34:30.000
%                               1:50:13.616 -> 1:50:00.000
% *) rrNum           - RR value (array of double)
% *) rrTime          - start time of un-rounded RR
% END OF IMPORT RR DATA

%% VALIDITY CHECK     数据有效性检验（A.注释、B.RR数据）fprintf('Data Validity Check:\n');% A. Annotation File Check
% *) generate annotation time acocrding to total epoch from header file
% *) result: anTimeGeneratedMat -> Matrix size: number of epoch X 3 (h,m,s)anTimeGeneratedMat = zeros(heaTotalEpoch, 3);for i=2:heaTotalEpochanTimeGeneratedMat(i, 3) =  anTimeGeneratedMat(i-1, 3) + SEC_PER_EPOCH;anTimeGeneratedMat(i, 2) = anTimeGeneratedMat(i-1, 2);anTimeGeneratedMat(i, 1) = anTimeGeneratedMat(i-1, 1);if anTimeGeneratedMat(i, 3) >= 60anTimeGeneratedMat(i, 3) = 0;anTimeGeneratedMat(i, 2) = anTimeGeneratedMat(i-1, 2) + 1;if anTimeGeneratedMat(i, 2) >= 60anTimeGeneratedMat(i, 2) = 0;anTimeGeneratedMat(i, 1) = anTimeGeneratedMat(i-1, 1) + 1;endendend% convert anTimeGeneratedMat into anTimeGeneratedCell for easier comparisonanTimeGeneratedCell = cell(size(anTimeGeneratedMat, 1), 1);for i = 1:size(anTimeGeneratedMat, 1)if anTimeGeneratedMat(i, 1) == 0 % when the 'hour' is 0anTimeGeneratedCell(i) = ...cellstr(strcat(sprintf('%d',anTimeGeneratedMat(i, 2)), ...sprintf(':%02d.000',anTimeGeneratedMat(i, 3))));elsetemp = strcat(sprintf('%d', anTimeGeneratedMat(i, 1)), ...sprintf(':%02d',anTimeGeneratedMat(i, 2)));anTimeGeneratedCell(i) = ...cellstr(strcat(temp, sprintf(':%02d.000',anTimeGeneratedMat(i, 3))));end
end% *) ANNOTATION FILE CHECK 1 (Total epoch of each data): fprintf('  CHECK 1: ');if heaTotalEpoch == size(anTime, 1) && ...size(unique(rrConvertedTime), 1) == size(anTimeGeneratedMat, 1) && ...heaTotalEpoch == size(unique(rrConvertedTime), 1)fprintf('[SUCCESS] heaTotalEpoch (%d) == size(anTime, 1) (%d) == ', ...'size(unique(rrConvertedTime), 1) (%d) == ', ...'size(anTimeGeneratedMat, 1) (%d)\n', heaTotalEpoch, size(anTime, 1), ...size(unique(rrConvertedTime), 1), size(anTimeGeneratedMat, 1));
elsefprintf('[WARNING] heaTotalEpoch (%d) != size(anTime, 1) (%d) != ', ...'size(unique(rrConvertedTime), 1) (%d) != ', ...'size(anTimeGeneratedMat, 1) (%d)\n', heaTotalEpoch, size(anTime, 1), ...size(unique(rrConvertedTime), 1), size(anTimeGeneratedMat, 1));
end% *) ANNOTATION FILE CHECK 2 (Check equality of anTimeGeneratedCell and anTime): fprintf('  CHECK 2: ');
if size(anTime, 1) == heaTotalEpochfor i=1:heaTotalEpochif ~strcmp(anTimeGeneratedCell{i}, anTime{i})fprintf('[FAILED ] anTimeGeneratedCell is NOT EQUAL to anTime\n');returnendendfprintf('[SUCCESS] anTimeGeneratedCell is EQUAL to anTime\n');
elsefprintf('[WARNING] size(anTime, 1) (%d) != heaTotalEpoch (%d), ', ...'anTimeGeneratedCell will be used\n', size(anTime, 1), heaTotalEpoch);end% *) ANNOTATION FILE CHECK 3 (Check annotation value must be '1', '2', '3','4', 'W', 'R', or {'MT', 'M' -> these two will be removed later}): fprintf('  CHECK 3: ');
distinctClass = char(unique(anClass));for i=1:size(distinctClass, 1)if distinctClass(i) ~= '1' && distinctClass(i) ~= '2' ...&& distinctClass(i) ~= '3' && distinctClass(i) ~= '4' ...&& distinctClass(i) ~= 'W' && distinctClass(i) ~= 'R' ...&& distinctClass(i) ~= 'M'fprintf('[WARNING ] Annotation values is NOT OK\n');returnendendfprintf('[SUCCESS] Annotation values is OK\n');% B. RR File Check
% *) RR FILE CHECK 1 (Check equality of size(unique(rrConvertedTime), 1)
% and heaTotalEpoch):fprintf('  CHECK 4: ');if size(unique(rrConvertedTime), 1) ~= heaTotalEpochfprintf('[WARNING] size(unique(rrConvertedTime), 1) (%d) != ', ...'heaTotalEpoch (%d)\n', size(unique(rrConvertedTime), 1), heaTotalEpoch);
elsefprintf('[SUCCESS] size(unique(rrConvertedTime), 1) (%d) == ', ...'heaTotalEpoch (%d)\n', size(unique(rrConvertedTime), 1), heaTotalEpoch);
end% END OF VALIDITY CHECK

%% SYNCHRONIZE RR AND ANNOTATION DATA   同步RR和注释数据epochCounter = 1;
rrCounter = 1;% rrCollection = each row contains RRs of associated epoch
rrCollection = cell(heaTotalEpoch, 1);% rrTimeCollection = each row contains RR time of associated epoch
rrTimeCollection = cell(heaTotalEpoch, 1);for i=1:size(rrConvertedTime, 1) % looping for each rrConvertedTime in that fileif strcmp(rrConvertedTime(i), anTimeGeneratedCell(epochCounter))% when i-th RR time is equal to annotation time of current epochrrCollection{epochCounter}(rrCounter) = rrNum(i);rrTimeCollection{epochCounter}(rrCounter) = rrTime(i);rrCounter=rrCounter+1;elseif ~strcmp(rrConvertedTime(i), anTimeGeneratedCell(epochCounter)) ...&& ~strcmp(rrConvertedTime(i), anTimeGeneratedCell(epochCounter+1))% when i-th RR time is not equal to annotation time of current epoch% and i-th RR time is not equal to annotation time of the next epochwhile ~strcmp(rrConvertedTime(i), anTimeGeneratedCell(epochCounter+1))epochCounter = epochCounter + 1;endrrCounter=1;epochCounter=epochCounter+1;rrCollection{epochCounter}(rrCounter) = rrNum(i);rrTimeCollection{epochCounter}(rrCounter) = rrTime(i);rrCounter=rrCounter+1;elseif ~strcmp(rrConvertedTime(i), anTimeGeneratedCell(epochCounter)) ...&& strcmp(rrConvertedTime(i), anTimeGeneratedCell(epochCounter+1))% when i-th RR time is not equal to annotation time of current epoch% and i-th RR time is equal to annotation time of the next epochrrCounter=1;epochCounter=epochCounter+1;rrCollection{epochCounter}(rrCounter) = rrNum(i);rrTimeCollection{epochCounter}(rrCounter) = rrTime(i);rrCounter=rrCounter+1;endend% END OF SYNCHRONIZE RR AND ANNOTATION DATA

%% SYNCHRONIZED DATA VALIDITY CHECK    同步数据有效性检验% generate new annotation matrix, fill the time without annotation 生成新的注释矩阵,在没有注释的情况下填充时间
anClassGeneratedCell = cell(heaTotalEpoch, 1);je = 1;for i=1:heaTotalEpochif strcmp(anTimeGeneratedCell(i), anTime(je))% if the time is the same, copy the annotationanClassGeneratedCell(i) = anClass(je);if je < size(anClass, 1)je = je + 1;endelse% if the time is different, fill with 'none'anClassGeneratedCell(i) = {'none'};endendfprintf('Removing invalid epoch:\n');
isExists = 0;for i=heaTotalEpoch:-1:1flag = 0;if sum(rrCollection{i}) < 28 || sum(rrCollection{i}) > 32% set flag to remove incomplete RR data of that epoch by:% check the sum of RR interval from each epoch,% can't be below 28 or higher than 32% (according to slp04 data, min sum is 29 and max is 30)flag = 1;fprintf('  Epoch %d (time: %s) of %s data is removed because ', ...'incomplete RR data\n', i, anTimeGeneratedCell{i}, fileName);elseif strcmp(anClassGeneratedCell(i), {'none'})% set flag to remove no annotation epochflag = 1;fprintf('  Epoch %d (time: %s) of %s data is removed because ', ...'no annotation\n', i, anTimeGeneratedCell{i}, fileName);elseif strcmp(anClassGeneratedCell(i), {'MT'}) || ...strcmp(anClassGeneratedCell(i), {'M'})% set flag to remove 'MT' or 'M' annotation epochflag = 1;fprintf('  Epoch %d (time: %s) of %s data is removed because ', ...'the annotation is %s\n', i, anTimeGeneratedCell{i}, fileName, ...anClassGeneratedCell{i});end% when the flag is 1, remove the dataif flag == 1anTimeGeneratedCell{i} = [];rrCollection{i} = [];anClassGeneratedCell{i} = [];isExists = 1;endend% print message if no invalid epoch
if ~isExistsfprintf('No invalid epoch\n');
end% delete empty row
anTimeGeneratedCell = ...anTimeGeneratedCell(~cellfun(@isempty, anTimeGeneratedCell));
rrCollection = rrCollection(~cellfun(@isempty, rrCollection));
anClassGeneratedCell = ...anClassGeneratedCell(~cellfun(@isempty, anClassGeneratedCell));%END OF SYNCHRONIZED DATA VALIDITY CHECK

%% PREPARE THE OUTPUT   输出信息OutputData = struct('filename', fileName, 'time', anTimeGeneratedCell, ...'rr', rrCollection, 'annotation', anClassGeneratedCell, 'age', age, ...'gender', gender, 'weight', weight);% END OF PREPARE THE OUTPUT
end

loadmatobject函数

——导入第index例样本的字段及对应数据

function Data = loadmatobject(fileName, index)%Load .mat object by index
%   Syntax:
%   Data = loadmatobject(dir, index)
%
%   Input:
%   *) fileName - file name
%   *) index    - index of .mat's variable to be returned
%
%   Output:
%   *) Data  - index-th variable returnedData = load(fileName);fieldName = fieldnames(Data);  %返回结构体的字段名称Data = Data.(fieldName{index});  %返回数据end

extractfeatures函数

——没有输出变量，输出非归一化、归一化的特征文件（'xlsx', 'mat'）

function extractfeatures(SlpdbData, destination, outputFormat)%Extract HRV Features
%   Syntax:
%   extractfeatures(SlpdbData, destination, outputFormat)
%
%   Input:
%   *) SlpdbData    - struct generated from importslpdb() function  % 从importslpdb()函数生成的结构
%   *) destination  - directory of the result
%   *) outputFormat - output format: 'xlsx', 'mat', 'all'
%
%   Output:
%   No output variables, but there are two files output:
%   hrv_features_unorm - unnormalized features
%   hrv_features_norm  - normalized features
%   target         - matrix total samples X 6 (1 - 6 classes target)nSamples = size(SlpdbData, 1);nClasses = length(unique([SlpdbData.annotation]));     % unique：找出数据矩阵中所有不重复数，确定分类类数hrv = zeros(nSamples, 25);target = zeros(nSamples, nClasses);target(:, [1 5]) = NaN;for i=1:nSamplesrr_diff = diff(SlpdbData(i).rr);                % diff函数式用于求导数和差分hrv(i, 1) = HRVFeature.AVNN(SlpdbData(i).rr);hrv(i, 2) = HRVFeature.SDNN(SlpdbData(i).rr);hrv(i, 3) = HRVFeature.RMSSD(rr_diff);hrv(i, 4) = HRVFeature.SDSD(rr_diff);hrv(i, 5) = HRVFeature.NNx(50, rr_diff);hrv(i, 6) = HRVFeature.PNNx(hrv(i, 5), size(SlpdbData(i).rr, 2));hrv(i, 7) = HRVFeature.HRV_TRIANGULAR_IDX(SlpdbData(i).rr);hrv(i, 8) = HRVFeature.SD1(hrv(i, 4));hrv(i, 9) = HRVFeature.SD2(hrv(i, 2), hrv(i, 4));hrv(i, 10) = HRVFeature.SD1_SD2_RATIO(hrv(i, 8), hrv(i, 9));hrv(i, 11) = HRVFeature.S(hrv(i, 8), hrv(i, 9));[TP,pLF,pHF,LFHFratio,VLF,LF,HF,f,Y,NFFT] = ...HRVFeature.fft_val_fun(SlpdbData(i).rr,2);hrv(i, 12) = TP;hrv(i, 13) = pLF;hrv(i, 14) = pHF;hrv(i, 15) = LFHFratio;hrv(i, 16) = VLF;hrv(i, 17) = LF;hrv(i, 18) = HF;% set class annotationswitch SlpdbData(i).annotationcase '1'target(i,6) = 1;        % 标记类别target(i,4) = 1;target(i,3) = 1;target(i,2) = 1;case '2'target(i,6) = 2;target(i,4) = 1;target(i,3) = 1;target(i,2) = 1;case '3'target(i,6) = 3;target(i,4) = 2;target(i,3) = 1;target(i,2) = 1;case '4'target(i,6) = 4;target(i,4) = 2;target(i,3) = 1;target(i,2) = 1;case 'R'target(i,6) = 5;target(i,4) = 3;target(i,3) = 2;target(i,2) = 1;case 'W'target(i,6) = 6;target(i,4) = 4;target(i,3) = 3;target(i,2) = 2;otherwisefprintf('Invalid Annotation');returnendendhrv( :, ~any(hrv,1) ) = [];% create a new dir if not existsdirList = dir;isDirExists = 0;for i=1:length(dir)if dirList(i).isdir && strcmp(dirList(i).name, destination)isDirExists = 1;endendif ~isDirExistsmkdir(destination);end% save the data into destinationhrv_features_unorm = hrv;hrv_features_norm = normalizedata(hrv, -1, 1);if strcmp(outputFormat, 'xlsx') || strcmp(outputFormat, 'all')    % strcmp是用于做字符串比较，保存成相应的特征文件xlswrite(strcat(destination, 'hrv_features_unorm.xlsx'), ...hrv_features_unorm);xlswrite(strcat(destination, 'hrv_features_norm.xlsx'), ...hrv_features_norm);xlswrite(strcat(destination, 'target.xlsx'), target);endif strcmp(outputFormat, 'mat') || strcmp(outputFormat, 'all')save(strcat(destination, 'hrv_features_unorm.mat'), 'hrv_features_unorm');save(strcat(destination, 'hrv_features_norm.mat'), 'hrv_features_norm');save(strcat(destination, 'target.mat'), 'target');endend

getindexrange函数

——通过index得到输入向量的范围

function range = getindexrange(nSamplesEachData, index)%Get range of inputted vector by index. For example [2 3 4 5] is the
%inputted nSamplesEachData and index is 2. Then, the output is [3 4 5].
%Explanation:
%The sum of [2 3 4 5] is 14 (there are 14 items).
%If index = 1, so the output is [1 2] -> total elements are 2
%If index = 2, so the output is [3 4 5] -> total elements are 3
%If index = 3, so the output is [6 7 8 9] -> total elements are 4
%If index = 4, so the output is [10 11 12 13 14] -> total elements are 5
%   Syntax:
%   range = getindexrange(nSamplesEachData, index)
%
%   Input:
%   *) nSamplesEachData - total number of data in each index   % 每个index中的数据总数
%   *) index            - index to be retrieved
%
%   Output:
%   *) range - a vector contains ordered number of associated index  % 包含相关index的有序数量向量if sum(index > length(nSamplesEachData)) >= 1disp('Index limit exceeded');returnendrange = [];for i=1:length(index)if index(i) == 1startNum = 1;endNum = nSamplesEachData(index(i));elsestartNum = sum(nSamplesEachData(1:index(i)-1)) + 1;endNum = startNum + nSamplesEachData(index(i)) - 1;endrange = [range startNum:endNum];end
end

PSOforELM函数

——PSO的初始化、结合ELM的PSO迭代

function [result, startTime, endTime] = PSOforELM(nFeatures, trainingData, ...testingData, PSOSettings)%Running PSO with ELM for feature selection and number of hidden nodes  % 使用ELM运行PSO来选择特征和隐藏节点的数量
%optimization
%   Syntax:
%   [result, startTime, endTime] = PSOforELM(nFeatures, trainingData, ...
%       testingData, PSOSettings)
%
%   Input:
%   *) nFeatures    - total number of features to be selected
%   *) trainingData - training data (Matrix size: total samples X nFeatures)
%   *) testingData  - testing data (Matrix size: total samples X nFeatures)
%   *) PSOSettings  - struct contains PSO parameters, examples:
%                       PSOSettings.MAX_ITERATION = MAX_ITERATION;
%                       PSOSettings.nParticles = 20;
%                       PSOSettings.W = 0.6;
%                       PSOSettings.c1 = 1.2;
%                       PSOSettings.c2 = 1.2;
%                       PSOSettings.Wa = 0.95;
%                       PSOSettings.Wf = 0.05;
%
%   Output:
%   *) result    - struct contains records of PSO ELM result   % 此结构包含PSO ELM结果的记录
%   *) startTime - time when the experiment starts
%   *) endTime   - time when the experiment endsstartTime = clock;%% PSO PARAMETER PREPARATION  PSO各参数定义及初始化% max total bits for hidden nodes
nHiddenBits = length(dectobin(size(trainingData, 1)));  % dectobin十-二进制转换
populationPosition = rand(PSOSettings.nParticles, nFeatures+nHiddenBits) > 0.5; % 随机矩阵中的每一个数与0.5比较，若值小于0.5，populationPosition矩阵中相对应的值返回1，否则返回0% 不符合条件的重新更新
for i=1:PSOSettings.nParticleswhile bintodec(populationPosition(i, nFeatures+1:end)) < nFeatures || ...bintodec(populationPosition(i, nFeatures+1:end)) > size(trainingData, 1) || ...sum(populationPosition(i, 1:nFeatures)) == 0populationPosition(i, :) = rand(1, nFeatures+nHiddenBits) > 0.5;   end
endpopulationVelocity = int64(zeros(PSOSettings.nParticles, 1)); % 初始化为0（十进制）% 定义pBest结构体的字段：position、fitness、trainingAccuracy、testingAccuracy
pBest(PSOSettings.nParticles).position = [];
pBest(PSOSettings.nParticles).fitness = [];
pBest(PSOSettings.nParticles).trainingAccuracy = [];
pBest(PSOSettings.nParticles).testingAccuracy = [];% 各字段初始化
for i=1:PSOSettings.nParticlespBest(i).position = false(1, nFeatures+nHiddenBits);   % max fitness valuepBest(i).fitness = repmat(-1000000, PSOSettings.nParticles, 1);pBest(i).trainingAccuracy = 0;pBest(i).testingAccuracy = 0;
end% 定义gBest结构体的字段：position、fitness、trainingAccuracy、testingAccuracy、fromIteration、fromParticle
gBest.position = false(1, nFeatures+nHiddenBits);
gBest.fitness = -1000000; % max fitness value all particle all iteration
gBest.trainingAccuracy = [];
gBest.testingAccuracy = [];
gBest.fromIteration = [];
gBest.fromParticle = [];% 定义result结构体的字段
result(PSOSettings.MAX_ITERATION+1).iteration = [];
result(PSOSettings.MAX_ITERATION+1).populationPosition = [];
result(PSOSettings.MAX_ITERATION+1).pBest = [];
result(PSOSettings.MAX_ITERATION+1).time = [];
result(PSOSettings.MAX_ITERATION+1).trainingAccuracy = [];
result(PSOSettings.MAX_ITERATION+1).testingAccuracy = [];
result(PSOSettings.MAX_ITERATION+1).model = [];
result(PSOSettings.MAX_ITERATION+1).gBest = [];% END OF PSO PARAMETER PREPARATION

%% INITIALIZATION STEP 初始化步骤
%fitness function evaluation 适应度函数的评价[modelArr, trainAccArr, testAccArr, timeArr, populationFitness, pBest] = ...evaluatefitness(PSOSettings, nFeatures, trainingData, testingData, ...populationPosition, pBest);         % evaluatefitness函数具体见本节最后gBest = gbestupdate(nFeatures, trainAccArr, testAccArr, populationFitness, ...populationPosition, gBest, 0);      % gbestupdate函数具体见本节最后% save initial data  保存原始数据
result(1).iteration = 0;
result(1).populationPosition = populationPosition;
result(1).pBest = pBest;
result(1).time = timeArr;
result(1).trainingAccuracy = trainAccArr;
result(1).testingAccuracy = testAccArr;%result(1).model = modelArr;
result(1).gBest = gBest;% END OF INITIALIZATION STEP

%% PSO ITERATION  PSO的迭代for iteration=1:PSOSettings.MAX_ITERATION% Update Velocityr1 = rand();r2 = rand();for i=1:PSOSettings.nParticles% calculate velocity value 计算速度值positionDec = int64(bintodec(populationPosition(i, :)));populationVelocity(i, 1) = PSOSettings.W * populationVelocity(i, 1) + ...PSOSettings.c1 * r1 * (bintodec(pBest(i).position) - positionDec) ...+ PSOSettings.c2 * r2 * (bintodec(gBest.position) - positionDec);% update particle positionnewPosDec = abs(int64(positionDec + populationVelocity(i, 1)));newPosBin = dectobin(newPosDec);% if the total bits is lower than nFeatures + nHiddenBits,add zeros in frontif size(newPosBin, 2) < (nFeatures + nHiddenBits)newPosBin = ...[zeros(1, (nFeatures + nHiddenBits) - size(newPosBin, 2)) ...newPosBin];end% if the number of hidden node is more than the number of samplesif bintodec(newPosBin(1, nFeatures+1:end)) > size(trainingData, 1) ...|| size(newPosBin(1, nFeatures+1:end), 2) > nHiddenBitsnewPosBin = ...[newPosBin(1, 1:nFeatures) dectobin(size(trainingData, 1))];end% if the number of selected features is 0 如果选择的特征数量为0while sum(newPosBin(1, 1:nFeatures)) == 0newPosBin(1, 1:nFeatures) = rand(1, nFeatures) > 0.5;end% set the new value of positionpopulationPosition(i, :) = newPosBin;end% fitness function evaluation[modelArr, trainAccArr, testAccArr, timeArr, populationFitness, pBest] = ...evaluatefitness(PSOSettings, nFeatures, trainingData, testingData, ...populationPosition, pBest);gBest = gbestupdate(nFeatures, trainAccArr, testAccArr, ...populationFitness, populationPosition, gBest, iteration+1);% save dataresult(iteration+1).iteration = iteration;result(iteration+1).populationPosition = populationPosition;result(iteration+1).pBest = pBest;result(iteration+1).time = timeArr;result(iteration+1).trainingAccuracy = trainAccArr;result(iteration+1).testingAccuracy = testAccArr;%result(iteration+1).model = modelArr;result(iteration+1).gBest = gBest;
end% END OF PSO ITERATION
endTime = clock;
end

PSOforSVM函数

function [result, startTime, endTime] = PSOforSVM(nFeatures, trainingData, ...testingData, PSOSettings)%Running PSO with SVM for feature selection
%   Syntax:
%   [result, startTime, endTime] = PSOforSVM(nFeatures, trainingData, ...
%       testingData, PSOSettings)
%
%   Input:
%   *) nFeatures    - total number of features to be selected
%   *) trainingData - training data (Matrix size: total samples X nFeatures)
%   *) testingData  - testing data (Matrix size: total samples X nFeatures)
%   *) PSOSettings  - struct contains PSO parameters, examples:
%                       PSOSettings.MAX_ITERATION = MAX_ITERATION;
%                       PSOSettings.nParticles = 20;
%                       PSOSettings.W = 0.6;
%                       PSOSettings.c1 = 1.2;
%                       PSOSettings.c2 = 1.2;
%                       PSOSettings.Wa = 0.95;
%                       PSOSettings.Wf = 0.05;
%
%   Output:
%   *) result    - struct contains records of PSO SVM result
%   *) startTime - time when the experiment starts
%   *) endTime   - time when the experiment endsstartTime = clock;%% PSO PARAMETER PREPARATION
populationPosition = rand(PSOSettings.nParticles, nFeatures) > 0.5;for i=1:PSOSettings.nParticleswhile sum(populationPosition(i, 1:nFeatures)) == 0populationPosition(i, :) = rand(1, nFeatures) > 0.5;end
endpopulationVelocity = int64(zeros(PSOSettings.nParticles, 1)); % in decimal value% pBest
pBest(PSOSettings.nParticles).position = [];
pBest(PSOSettings.nParticles).fitness = [];
pBest(PSOSettings.nParticles).trainingAccuracy = [];
pBest(PSOSettings.nParticles).testingAccuracy = [];for i=1:PSOSettings.nParticlespBest(i).position = false(1, nFeatures);% max fitness valuepBest(i).fitness = repmat(-1000000, PSOSettings.nParticles, 1);pBest(i).trainingAccuracy = 0;pBest(i).testingAccuracy = 0;
end% gBest
gBest.position = false(1, nFeatures);
gBest.fitness = -1000000; % max fitness value all particle all iteration
gBest.trainingAccuracy = [];
gBest.testingAccuracy = [];
gBest.fromIteration = [];
gBest.fromParticle = [];% initialize struct data
result(PSOSettings.MAX_ITERATION+1).iteration = [];
result(PSOSettings.MAX_ITERATION+1).populationPosition = [];
result(PSOSettings.MAX_ITERATION+1).pBest = [];
result(PSOSettings.MAX_ITERATION+1).time = [];
result(PSOSettings.MAX_ITERATION+1).trainingAccuracy = [];
result(PSOSettings.MAX_ITERATION+1).testingAccuracy = [];
result(PSOSettings.MAX_ITERATION+1).gBest = [];% END OF PSO PARAMETER PREPARATION

%% INITIALIZATION STEP%fitness function evaluation
[trainAccArr, testAccArr, timeArr, populationFitness, pBest] = ...evaluatefitness(PSOSettings, nFeatures, trainingData, testingData, ...populationPosition, pBest);gBest = gbestupdate(nFeatures, trainAccArr, testAccArr, populationFitness, ...populationPosition, gBest, 0);% save initial data
result(1).iteration = 0;
result(1).populationPosition = populationPosition;
result(1).pBest = pBest;
result(1).time = timeArr;
result(1).trainingAccuracy = trainAccArr;
result(1).testingAccuracy = testAccArr;
result(1).gBest = gBest;% END OF INITIALIZATION STEP

%% PSO ITERATIONfor iteration=1:PSOSettings.MAX_ITERATION% Update Velocityr1 = rand();r2 = rand();for i=1:PSOSettings.nParticles% calculate velocity valuepositionDec = int64(bintodec(populationPosition(i, :)));populationVelocity(i, 1) = PSOSettings.W * populationVelocity(i, 1) + ...PSOSettings.c1 * r1 * (bintodec(pBest(i).position) - positionDec) ...+ PSOSettings.c2 * r2 * (bintodec(gBest.position) - positionDec);% update particle positionnewPosDec = abs(int64(positionDec + populationVelocity(i, 1)));newPosBin = dectobin(newPosDec);% if the total bits is lower than nFeatures, add zeros in frontif size(newPosBin, 2) < nFeaturesnewPosBin = [zeros(1, (nFeatures)-size(newPosBin, 2)) newPosBin];end% if the total bits is higher than nFeatures, get first nFeaturesif size(newPosBin, 2) > nFeaturesnewPosBin = newPosBin(1, 1:nFeatures);end% if the number of selected features is 0while sum(newPosBin(1, 1:nFeatures)) == 0newPosBin(1, 1:nFeatures) = rand(1, nFeatures) > 0.5;end% set the new value of positionpopulationPosition(i, :) = newPosBin;end% fitness function evaluation[trainAccArr, testAccArr, timeArr, populationFitness, pBest] = ...evaluatefitness(PSOSettings, nFeatures, trainingData, testingData, ...populationPosition, pBest);gBest = gbestupdate(nFeatures, trainAccArr, testAccArr, ...populationFitness, populationPosition, gBest, iteration+1);% save dataresult(iteration+1).iteration = iteration;result(iteration+1).populationPosition = populationPosition;result(iteration+1).pBest = pBest;result(iteration+1).time = timeArr;result(iteration+1).trainingAccuracy = trainAccArr;result(iteration+1).testingAccuracy = testAccArr;result(iteration+1).gBest = gBest;
end% END OF PSO ITERATION
endTime = clock;end

extractresults函数

——没有输出变量，有Excel输出文件（每次实验结果、所有实验的最好结果）、保存gBest随迭代变化图片

function extractresults(resultRootFolder, nFeatures, classNum, nExperiments, ...nIterations)%Extract raw results of experiment using PSOELM or PSOSVM method
%   Syntax:
%   extractresults(resultRootFolder, nFeatures, classNum, nExperiments, ...
%       nIterations)
%
%   Input:
%   *) resultRootFolder - root directory of the results
%   *) nFeatures        - total features
%   *) classNum         - number of class in vector -> [2 3 4 6]
%   *) nExperiments     - total experiments of PSO
%   *) nIterations      - total iterations of each PSO
%
%   Output:
%   No output variables, there are excel output files:
%   [method]_[filename]_extracted_result - results of each experiment
%   [method]_result                      - best of all experiments
%
%   Parameter Example:
%   resultRootFolder = 'PSOELM_raw_result';
%   nFeatures = 18;
%   classNum = [2 3 4 6];
%   nExperiments = 25;
%   nIterations = 100;fileNames = {'slp01a' 'slp01b' 'slp02a' 'slp02b' 'slp03' 'slp04' ...'slp14' 'slp16' 'slp32' 'slp37' 'slp41' 'slp45' 'slp48' ...'slp59' 'slp60' 'slp61' 'slp66' 'slp67x'};    method = strsplit(resultRootFolder, '_');method = method{1};   % PSOSVM | PSOELMnClassClassifiers = length(classNum);   % 4headerEachExp = [];switch methodcase 'PSOELM'headerEachExp = {'Experiment', 'gBestFitness', 'TrainAcc', ...'TestAcc', 'ProcessTime(Sec)', 'HiddenNodes', 'SelectedFeatures'};case 'PSOSVM'headerEachExp = {'Experiment', 'gBestFitness', 'TrainAcc', ...'TestAcc', 'ProcessTime(Sec)', 'SelectedFeatures'};endheaderBestExp = headerEachExp;headerBestExp{1} = 'RecordingName';% write header for main excel result (only 1 excel)for i=1:length(classNum)xlswrite(sprintf('%s/%s_result.xlsx', resultRootFolder, method), ...headerBestExp, sprintf('%d classes', classNum(i)));endfor iFile=1:length(fileNames) % loop for each file% eachFileFolder example: PSOELM_raw_result/PSOELM_slp01a_raw_resulteachFileFolder = sprintf('%s/%s_%s_raw_result', resultRootFolder, ...method, fileNames{iFile});for iClass=1:nClassClassifiers % loop for each class numbermatFileName = sprintf('%s/%s_%s_%dclasses_raw_result.mat', ...eachFileFolder, method, fileNames{iFile}, classNum(iClass));ExperimentResult = loadmatobject(matFileName, 1);temp = zeros(nExperiments, length(headerEachExp)-1);nBits = ...length(ExperimentResult(4). ...iterationResult(end).gBest.position) - nFeatures;gBestParticles = false(nExperiments, nFeatures+nBits);tempCell = cell(nExperiments, 1);% get the last gBest result of each experimentfor iExp=1:nExperimentslastResult = ExperimentResult(iExp).iterationResult(end);temp(iExp, 1) = iExp;temp(iExp, 2) = lastResult.gBest.fitness;temp(iExp, 3) = lastResult.gBest.trainingAccuracy;temp(iExp, 4) = lastResult.gBest.testingAccuracy;temp(iExp, 5) = ...etime(ExperimentResult(iExp).endTime, ...ExperimentResult(iExp).startTime);if strcmp(method, 'PSOELM')temp(iExp, 6) = ...bintodec(lastResult.gBest.position(nFeatures+1:end));endtempCell(iExp, 1) = ...{bintostringorder(lastResult.gBest.position(1, 1:nFeatures))};gBestParticles(iExp, :) = lastResult.gBest.position;endeachFileExcelPath = sprintf('%s/%s_%s_extracted_result.xlsx', ...eachFileFolder, method, fileNames{iFile});xlswrite(eachFileExcelPath, headerEachExp, ...sprintf('%d classes', classNum(iClass)), 'A1');xlswrite(eachFileExcelPath, temp, ...sprintf('%d classes', classNum(iClass)), 'A2');xlswrite(eachFileExcelPath, tempCell, ...sprintf('%d classes', classNum(iClass)), ...sprintf('%s2', getexcelcolumncode(length(headerEachExp))));% get the best experiment of each classificationbestExpIdx = find(temp(:, 2) == max(temp(:, 2)));if length(bestExpIdx) > 1% if have the same gBest fitness value, get the max of testAccbestExpIdx = ...bestExpIdx(temp(bestExpIdx, 4) == max(temp(bestExpIdx, 4)));if length(bestExpIdx) > 1% if have the same testAcc, get the max of trainAccbestExpIdx = bestExpIdx(temp(bestExpIdx, 3) == ...max(temp(bestExpIdx, 3)));if length(bestExpIdx) > 1% if have the same trainAcc,% get the min of selected featuresbestExpIdx = bestExpIdx( ...sum(gBestParticles(bestExpIdx, 1:nFeatures), 2) == ...min(sum(gBestParticles(bestExpIdx, 1:nFeatures), 2)));if length(bestExpIdx) > 1% if have the same selected feature,% check the method usedswitch methodcase 'PSOELM'bestExpIdx = ...bestExpIdx(temp(bestExpIdx, 6) == ...min(temp(bestExpIdx, 6)));if length(bestExpIdx) > 1% if have the same selected feature,% get the firstbestExpIdx = bestExpIdx(1);endcase 'PSOSVM'bestExpIdx = bestExpIdx(1);endendendendend% mark the best indexxlswrite(eachFileExcelPath, {'BEST EXPERIMENT'}, ...sprintf('%d classes', classNum(iClass)), ...sprintf('%s%d', ...getexcelcolumncode(length(headerEachExp)+1), bestExpIdx+1));%{% gather gBest fitness of the best experimentgBest = zeros(nIterations, 1);for iItr=1:nIterationsgBest(iItr) = ...ExperimentResult(bestExpIdx). ...iterationResult(iItr+1).gBest.fitness;end% save graphicsf = figure;plot(1:nIterations, gBest);ylabel('gBest Fitness'); xlabel('Iteration');title(sprintf('[%s] Best Experiment of %s (%d classes)', ...method, fileNames{iFile}, classNum(iClass)));saveas(f, ...sprintf('%s/[%s] Best Experiment of %s (%d classes).png', ...eachFileFolder, method, fileNames{iFile}, classNum(iClass)));close all;%}% save result to main excelswitch methodcase 'PSOELM'xlswrite( ...sprintf('%s/%s_result.xlsx', resultRootFolder, ...method), ...[fileNames(iFile) temp(bestExpIdx, 2) ...temp(bestExpIdx, 3) temp(bestExpIdx, 4) ...temp(bestExpIdx, 5) temp(bestExpIdx, 6) ...tempCell(bestExpIdx)], ...sprintf('%d classes', classNum(iClass)), ...sprintf('A%d', iFile+1));case 'PSOSVM'xlswrite( ...sprintf('%s/%s_result.xlsx', resultRootFolder, ...method), ...[fileNames(iFile) temp(bestExpIdx, 2) ...temp(bestExpIdx, 3) temp(bestExpIdx, 4) ...temp(bestExpIdx, 5) tempCell(bestExpIdx)], ...sprintf('%d classes', classNum(iClass)), ...sprintf('A%d', iFile+1));endendendend

evaluatefitness、gbestupdate函数

function [modelArr, trainAccArr, testAccArr, timeArr, populationFitness, ...pBest] = evaluatefitness(PSOSettings, nFeatures, trainingData, ...testingData, populationPosition, pBest)modelArr(PSOSettings.nParticles).inputWeight = [];   % 定义modelArr结构体字段inputWeight、outputWeightmodelArr(PSOSettings.nParticles).outputWeight = [];  trainAccArr = zeros(PSOSettings.nParticles, 1);testAccArr = zeros(PSOSettings.nParticles, 1);timeArr = zeros(PSOSettings.nParticles, 1);populationFitness = zeros(PSOSettings.nParticles, 1);for i=1:PSOSettings.nParticlestic;% TRAININGmaskedTrainingFeature = featuremasking(trainingData, ...  % featuremasking函数具体见最后populationPosition(i, 1:nFeatures));    % remove unselected features% prepare the target data% (example: transformation from 4 into [0 0 0 1 0 0])trainingTarget = full(ind2vec(trainingData(:,end)'))'; [Model, trainAcc] = trainELM(maskedTrainingFeature, trainingTarget, ...  % trainELM函数具体见最后bintodec(populationPosition(i, nFeatures+1:end)));% TESTINGmaskedTestingFeature = featuremasking(testingData, ...    % featuremasking函数具体见最后populationPosition(i, 1:nFeatures)); % remove unselected features% prepare the target data% (example: transformation from 4 into [0 0 0 1 0 0])testingTarget = full(ind2vec(testingData(:,end)'))';testAcc = testELM(maskedTestingFeature, testingTarget, Model);    % testELM函数具体见最后populationFitness(i, 1) = fitness(PSOSettings.Wa, PSOSettings.Wf, ... testAcc, populationPosition(i, 1:nFeatures));     % function fitnessValue = fitness(Wa, Wf, acc, featureMask)% fitnessValue = Wa * acc + Wf * (1 - (sum(featureMask)/length(featureMask)));% end% pBest updateischanged = 0;% 满足以下任意一条件，更新变量if populationFitness(i, 1) > pBest(i).fitnessischanged = 1;elseif populationFitness(i, 1) == pBest(i).fitnessif pBest(i).testingAccuracy < testAccischanged = 1;elseif pBest(i).trainingAccuracy < trainAccischanged = 1;elseif sum(pBest(i).position(1, 1:nFeatures)) > ...sum(populationPosition(i, 1:nFeatures))ischanged = 1;elseif bintodec(pBest(i).position(1, nFeatures+1:end)) > ...bintodec(populationPosition(i, nFeatures+1:end))ischanged = 1;endendif ischangedpBest(i).fitness = populationFitness(i, 1);pBest(i).position = populationPosition(i, :);pBest(i).trainingAccuracy = trainAcc;pBest(i).testingAccuracy = testAcc;end% end of pBest updatemodelArr(i) = Model;timeArr(i) = toc;trainAccArr(i) = trainAcc;testAccArr(i) = testAcc;end
end

function gBest = gbestupdate(nFeatures, trainAccArr, testAccArr, ...populationFitness, populationPosition, gBest, iteration)if max(populationFitness) >= gBest.fitnessfound = find(populationFitness == max(populationFitness));if length(found) > 1% if have the same gBest fitness value, get the max of testAcc 如果具有相同的gBest适应度值，则获取testAcc的最大值found = found(testAccArr(found) == max(testAccArr(found)));if length(found) > 1% if have the same testAcc, get the max of trainAcc  如果有相同的testAcc，得到trainAcc的最大值found = found(trainAccArr(found) == max(trainAccArr(found)));if length(found) > 1% if have the same trainAcc, get the min of selected features 如果具有相同的trainAcc，则获取所选特征的最小值found = ...found(sum(populationPosition(found, 1:nFeatures), 2) ...== min(sum(populationPosition(found, 1:nFeatures), 2)));if length(found) > 1% if have the same selected feature,get the min of hidden nodehn = zeros(length(found), 1);for i=1:length(found)hn(i, 1) = bintodec(populationPosition(found(i), ...nFeatures+1:end));endfound = found(hn == min(hn));if length(found) > 1found = found(1);endendendendendgBest.fitness = populationFitness(found);gBest.position = populationPosition(found, :);gBest.trainingAccuracy = trainAccArr(found);gBest.testingAccuracy = testAccArr(found);gBest.fromIteration = iteration;gBest.fromParticle = found;end
end

function [trainAccArr, testAccArr, timeArr, populationFitness, pBest] = ...evaluatefitness(PSOSettings, nFeatures, trainingData, testingData, ...populationPosition, pBest)trainAccArr = zeros(PSOSettings.nParticles, 1);testAccArr = zeros(PSOSettings.nParticles, 1);timeArr = zeros(PSOSettings.nParticles, 1);populationFitness = zeros(PSOSettings.nParticles, 1);for i=1:PSOSettings.nParticlestic;% TRAININGmaskedTrainingFeature = featuremasking(trainingData, ...populationPosition(i, 1:nFeatures)); % remove unselected featuresModel = trainSVM(maskedTrainingFeature, trainingData(:,end), 'RBF');trainAcc = testSVM(maskedTrainingFeature, trainingData(:,end), Model);% TESTINGmaskedTestingFeature = featuremasking(testingData, ...populationPosition(i, 1:nFeatures));     % remove unselected featurestestAcc = testSVM(maskedTestingFeature, testingData(:,end), Model);populationFitness(i, 1) = fitness(PSOSettings.Wa, PSOSettings.Wf, ...testAcc, populationPosition(i, 1:nFeatures));% pBest updateischanged = 0;if populationFitness(i, 1) > pBest(i).fitnessischanged = 1;elseif populationFitness(i, 1) == pBest(i).fitnessif pBest(i).testingAccuracy < testAccischanged = 1;elseif pBest(i).trainingAccuracy < trainAccischanged = 1;elseif sum(pBest(i).position(1, 1:nFeatures)) > ...sum(populationPosition(i, 1:nFeatures))ischanged = 1;endendif ischangedpBest(i).fitness = populationFitness(i, 1);pBest(i).position = populationPosition(i, :);pBest(i).trainingAccuracy = trainAcc;pBest(i).testingAccuracy = testAcc;end% end of pBest updatetimeArr(i) = toc;trainAccArr(i) = trainAcc;testAccArr(i) = testAcc;end
endfunction gBest = gbestupdate(nFeatures, trainAccArr, testAccArr, ...populationFitness, populationPosition, gBest, iteration)if max(populationFitness) >= gBest.fitnessfound = find(populationFitness == max(populationFitness));if length(found) > 1% if have the same gBest fitness value, get the max of testAccfound = found(testAccArr(found) == max(testAccArr(found)));if length(found) > 1% if have the same testAcc, get the max of trainAccfound = found(trainAccArr(found) == max(trainAccArr(found)));if length(found) > 1% if have the same trainAcc, get the min of selected featuresfound = ...found(sum(populationPosition(found, 1:nFeatures), 2) ...== min(sum(populationPosition(found, 1:nFeatures), 2)));if length(found) > 1% if have the same selected feature, get the firstfound = found(1);endendendendgBest.fitness = populationFitness(found);gBest.position = populationPosition(found, :);gBest.trainingAccuracy = trainAccArr(found);gBest.testingAccuracy = testAccArr(found);gBest.fromIteration = iteration;gBest.fromParticle = found;end
end

featuremasking、trainELM、testELM函数

function maskedFeature = featuremasking(feature, mask)
%Retrieve masked features
%   Syntax:
%   maskedFeature = featuremasking(feature, mask)
%
%   Input:
%   *) feature       - feature collection
%      (Matrix size: total samples X total features)
%   *) mask          - logical matrix, 1 means selected, 0 is not selected
%      (Matrix size: 1 X total features)
%   Output:
%   *) maskedFeature - matrix with selected features only
%      (Matrix size: total samples X total selected featuresmaskedFeature = zeros(size(feature, 1), sum(mask));j = 1;for i=1:sum(mask)if mask(1, i) == 1maskedFeature(:, j) = feature(:, i);j = j + 1;endend
endfunction [ELMModel, trainAcc] = trainELM(feature, target, nHiddenNode)
%Train Extreme Learning Machine (ELM) model
%   Syntax:
%   [ELMModel, trainAcc] = trainELM(feature, target, nHiddenNode)
%
%   Input:
%   *) feature     - feature collection
%      (Matrix size: total samples X total features)
%   *) target      - target of each sample
%      (Matrix size: total samples X total classes)
%       Example: class 4 -> target is [0 0 0 1]
%   *) nHiddenNode - total hidden nodes of ELM
%
%   Output:
%   *) ELMModel.inputWeight  - input weight of ELM
%      (Matrix size: nHiddenNode (+1 for bias) X total features)
%   *) ELMModel.outputWeight - output weight of ELM
%      (Matrix size: total classes X nHiddenNode)
%   *) trainAcc              - training accuracyif size(feature, 2) == 0fprintf('Someting went wrong, no feature selected.');returnend% STEP 1: RANDOMLY ASSIGN INPUT WEIGHT AND BIAS 随机分配输入权重和偏差minWeight = -1;maxWeight = 1;inputWeight = (maxWeight-minWeight) .* ...rand(nHiddenNode, size(feature, 2)+1) + minWeight;% STEP 2: CALCULATE THE HIDDEN LAYER OUTPUT MATRIX H  计算隐含层输出矩阵H% linear combination of hidden output 线性组合的隐藏输出hiddenOutput = (inputWeight(:, 1:end-1) * feature')+ ...repmat(inputWeight(:, end), 1, size(feature, 1));% apply activation function on hidden output 对隐藏输出应用激活函数hiddenOutput = sigmoid(hiddenOutput);  % function y = sigmoid(x)% y = 1./(1 + exp(-1.*x));% end% STEP 3: CALCULATE THE OUTPUT WEIGHT B   计算输出权值B% estimate output weightoutputWeight = target' * pinv(hiddenOutput);   % pinv函数-求矩阵的伪逆矩阵% STEP 4: APPLY MODEL TO TRAINING DATA   将模型应用训练数据% linear combination of predicted output  预测输出的线性组合predictedOutput = outputWeight * hiddenOutput;% apply activation function on predicted output  对预测输出应用激活函数predictedOutput = sigmoid(predictedOutput);maxPred = max(predictedOutput); predictedClass = zeros(size(predictedOutput, 2), 1);for i=1:size(predictedOutput, 2)class = find(predictedOutput(:, i) == maxPred(i));predictedClass(i) = class(1, 1);endtrainAcc = sum(predictedClass == vec2ind(target')')/ ... % vec2ind：向量到索引，向量中每一列中的所有元素有且只能有一个为1，向量的列数即为索引矩阵的列数size(predictedOutput, 2) * 100;ELMModel.inputWeight = inputWeight;ELMModel.outputWeight = outputWeight;
endfunction testAcc = testELM(feature, target, ELMModel)
%Test Extreme Learning Machine (ELM) model
%   Syntax:
%   testAcc = testELM(feature, target, ELMModel)
%
%   Input:
%   *) feature  - feature collection
%      (Matrix size: total samples X total features)
%   *) target   - target of each sample
%      (Matrix size: total samples X total classes)
%      Example: class 4 -> target is [0 0 0 1]
%   *) ELMModel - ELMModel generated from trainELM() function
%
%   Output:
%   *) testAcc  - testing accuracy% linear combination of hidden outputhiddenOutput = (ELMModel.inputWeight(:, 1:end-1) * feature')+ ...repmat(ELMModel.inputWeight(:, end), 1, size(feature, 1)); % apply activation function on hidden outputhiddenOutput = sigmoid(hiddenOutput);% linear combination of predicted outputpredictedOutput = ELMModel.outputWeight * hiddenOutput;% apply activation function on predicted outputpredictedOutput = sigmoid(predictedOutput);maxPred = max(predictedOutput);predictedClass = zeros(size(predictedOutput, 2), 1);for i=1:size(predictedOutput, 2)class = find(predictedOutput(:, i) == maxPred(i));predictedClass(i) = class(1, 1);endtestAcc = sum(predictedClass == vec2ind(target')')/ ...size(predictedOutput, 2) * 100;
end