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

最新推荐文章于 2023-12-26 21:37:13 发布

名字太简单就记不住

最新推荐文章于 2023-12-26 21:37:13 发布

阅读量1.1k

点赞数

分类专栏：毕设

本文链接：https://blog.csdn.net/weixin_42920780/article/details/104363733

版权

毕设专栏收录该内容

11 篇文章 12 订阅

订阅专栏

本文详细介绍了如何利用粒子群优化(PSO)算法结合极限学习机(ELM)和支援向量机(SVM)进行特征选择和模型优化，通过实例演示了PSO算法在ELM和SVM中的初始化、迭代过程，以及结果提取和评估的方法。

摘要由CSDN通过智能技术生成

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 annotation

fileName = 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 file

if 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;
end

age = 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 char
    rrStartTimeChar(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 grouping
    rrWhichGroup = floor(rrSecond/SEC_PER_EPOCH)*SEC_PER_EPOCH;

    % set epoch grouping
    if rrWhichGroup == 0
        rrStartTimeChar(end-5) = 48;
    elseif rrWhichGroup == 30
        rrStartTimeChar(end-5) = 51;
    end

    rrStartTimeChar(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)));
end

rrTime = 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:heaTotalEpoch
   anTimeGeneratedMat(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) >= 60
       anTimeGeneratedMat(i, 3) = 0;
       anTimeGeneratedMat(i, 2) = anTimeGeneratedMat(i-1, 2) + 1;

       if anTimeGeneratedMat(i, 2) >= 60
           anTimeGeneratedMat(i, 2) = 0;
           anTimeGeneratedMat(i, 1) = anTimeGeneratedMat(i-1, 1) + 1;
       end

   end

end

% convert anTimeGeneratedMat into anTimeGeneratedCell for easier comparison

anTimeGeneratedCell = cell(size(anTimeGeneratedMat, 1), 1);

for i = 1:size(anTimeGeneratedMat, 1)
    if anTimeGeneratedMat(i, 1) == 0 % when the 'hour' is 0
        anTimeGeneratedCell(i) = ...
            cellstr(strcat(sprintf('%d',anTimeGeneratedMat(i, 2)), ...
            sprintf(':%02d.000',anTimeGeneratedMat(i, 3))));
    else
        temp = 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));
else

    fprintf('[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) == heaTotalEpoch

    for i=1:heaTotalEpoch

        if ~strcmp(anTimeGeneratedCell{i}, anTime{i})
            fprintf('[FAILED ] anTimeGeneratedCell is NOT EQUAL to anTime\n');
            return
        end

    end

    fprintf('[SUCCESS] anTimeGeneratedCell is EQUAL to anTime\n');
else

    fprintf('[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');
        return
    end

end

fprintf('[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) ~= heaTotalEpoch
    fprintf('[WARNING] size(unique(rrConvertedTime), 1) (%d) != ', ...
        'heaTotalEpoch (%d)\n', size(unique(rrConvertedTime), 1), heaTotalEpoch);
else
    fprintf('[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 file

    if strcmp(rrConvertedTime(i), anTimeGeneratedCell(epochCounter))
        % when i-th RR time is equal to annotation time of current epoch
        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 not equal to annotation time of the next epoch

        while ~strcmp(rrConvertedTime(i), anTimeGeneratedCell(epochCounter+1))
            epochCounter = epochCounter + 1;
        end

        rrCounter=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 epoch

        rrCounter=1;
        epochCounter=epochCounter+1;
        rrCollection{epochCounter}(rrCounter) = rrNum(i);
        rrTimeCollection{epochCounter}(rrCounter) = rrTime(i);
        rrCounter=rrCounter+1;
    end

end

% 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:heaTotalEpoch

   if strcmp(anTimeGeneratedCell(i), anTime(je))

	   % if the time is the same, copy the annotation
       anClassGeneratedCell(i) = anClass(je);

       if je < size(anClass, 1)
           je = je + 1;
       end
   else
	   % if the time is different, fill with 'none'
       anClassGeneratedCell(i) = {'none'};
   end

end

fprintf('Removing invalid epoch:\n');
isExists = 0;

for i=heaTotalEpoch:-1:1
    flag = 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 epoch
        flag = 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 epoch
        flag = 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 data
    if flag == 1
        anTimeGeneratedCell{i} = [];
        rrCollection{i} = [];
        anClassGeneratedCell{i} = [];
        isExists = 1;
    end

end

% print message if no invalid epoch
if ~isExists
   fprintf('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 returned

    Data = 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:nSamples
        rr_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 annotation
        switch SlpdbData(i).annotation
            case '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;
            otherwise
                fprintf('Invalid Annotation');
                return
        end
    end

    hrv( :, ~any(hrv,1) ) = [];

    % create a new dir if not exists
    dirList = dir;
    isDirExists = 0;

    for i=1:length(dir)
        if dirList(i).isdir && strcmp(dirList(i).name, destination)
            isDirExists = 1;
        end
    end

    if ~isDirExists
        mkdir(destination);
    end

    % save the data into destination
    hrv_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);
    end

    if 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');
    end
    
end

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)) >= 1
        disp('Index limit exceeded');
        return
    end

    range = [];

    for i=1:length(index)
        if index(i) == 1
            startNum = 1;
            endNum = nSamplesEachData(index(i));
        else
            startNum = sum(nSamplesEachData(1:index(i)-1)) + 1;
            endNum = startNum + nSamplesEachData(index(i)) - 1;
        end
        range = [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 ends

startTime = 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.nParticles
    while bintodec(populationPosition(i, nFeatures+1:end)) < nFeatures || ...
          bintodec(populationPosition(i, nFeatures+1:end)) > size(trainingData, 1) || ...
          sum(populationPosition(i, 1:nFeatures)) == 0
        populationPosition(i, :) = rand(1, nFeatures+nHiddenBits) > 0.5;   
    end
end

populationVelocity = 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.nParticles
    pBest(i).position = false(1, nFeatures+nHiddenBits);   
    % max fitness value
    pBest(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 Velocity
    r1 = 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 position
        newPosDec = abs(int64(positionDec + populationVelocity(i, 1)));
        newPosBin = dectobin(newPosDec);

        % if the total bits is lower than nFeatures + nHiddenBits,add zeros in front
        if 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 samples
        if bintodec(newPosBin(1, nFeatures+1:end)) > size(trainingData, 1) ...
                || size(newPosBin(1, nFeatures+1:end), 2) > nHiddenBits
            newPosBin = ...
                [newPosBin(1, 1:nFeatures) dectobin(size(trainingData, 1))];
        end

        % if the number of selected features is 0 如果选择的特征数量为0
        while sum(newPosBin(1, 1:nFeatures)) == 0
            newPosBin(1, 1:nFeatures) = rand(1, nFeatures) > 0.5;
        end

        % set the new value of position
        populationPosition(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 data
    result(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 ends

startTime = clock;

%% PSO PARAMETER PREPARATION
populationPosition = rand(PSOSettings.nParticles, nFeatures) > 0.5;

for i=1:PSOSettings.nParticles
    while sum(populationPosition(i, 1:nFeatures)) == 0
        populationPosition(i, :) = rand(1, nFeatures) > 0.5;
    end
end

populationVelocity = 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.nParticles
    pBest(i).position = false(1, nFeatures);

    % max fitness value
    pBest(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 ITERATION

for iteration=1:PSOSettings.MAX_ITERATION
    % Update Velocity
    r1 = 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 position
        newPosDec = abs(int64(positionDec + populationVelocity(i, 1)));
        newPosBin = dectobin(newPosDec);

        % if the total bits is lower than nFeatures, add zeros in front
        if size(newPosBin, 2) < nFeatures
            newPosBin = [zeros(1, (nFeatures)-size(newPosBin, 2)) newPosBin];
        end

        % if the total bits is higher than nFeatures, get first nFeatures
        if size(newPosBin, 2) > nFeatures
            newPosBin = newPosBin(1, 1:nFeatures);
        end

        % if the number of selected features is 0
        while sum(newPosBin(1, 1:nFeatures)) == 0
            newPosBin(1, 1:nFeatures) = rand(1, nFeatures) > 0.5;
        end

        % set the new value of position
        populationPosition(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 data
    result(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 | PSOELM
    nClassClassifiers = length(classNum);   % 4
    headerEachExp = [];

    switch method
        case 'PSOELM'
            headerEachExp = {'Experiment', 'gBestFitness', 'TrainAcc', ...
                'TestAcc', 'ProcessTime(Sec)', 'HiddenNodes', 'SelectedFeatures'};
        case 'PSOSVM'
            headerEachExp = {'Experiment', 'gBestFitness', 'TrainAcc', ...
                'TestAcc', 'ProcessTime(Sec)', 'SelectedFeatures'};
    end

    headerBestExp = 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)));
    end

    for iFile=1:length(fileNames) % loop for each file
        % eachFileFolder example: PSOELM_raw_result/PSOELM_slp01a_raw_result
        eachFileFolder = sprintf('%s/%s_%s_raw_result', resultRootFolder, ...
            method, fileNames{iFile});

        for iClass=1:nClassClassifiers % loop for each class number
            matFileName = 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 experiment
            for iExp=1:nExperiments
                lastResult = 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));
                end

                tempCell(iExp, 1) = ...
                    {bintostringorder(lastResult.gBest.position(1, 1:nFeatures))};
                gBestParticles(iExp, :) = lastResult.gBest.position;
            end

            eachFileExcelPath = 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 classification
            bestExpIdx = find(temp(:, 2) == max(temp(:, 2)));

            if length(bestExpIdx) > 1
                % if have the same gBest fitness value, get the max of testAcc
                bestExpIdx = ...
                    bestExpIdx(temp(bestExpIdx, 4) == max(temp(bestExpIdx, 4)));
                if length(bestExpIdx) > 1
                    % if have the same testAcc, get the max of trainAcc
                    bestExpIdx = bestExpIdx(temp(bestExpIdx, 3) == ...
                        max(temp(bestExpIdx, 3)));
                    if length(bestExpIdx) > 1
                        % if have the same trainAcc,
                        % get the min of selected features
                        bestExpIdx = 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 used
                            switch method
                                case 'PSOELM'
                                    bestExpIdx = ...
                                        bestExpIdx(temp(bestExpIdx, 6) == ...
                                        min(temp(bestExpIdx, 6)));
                                    if length(bestExpIdx) > 1
                                        % if have the same selected feature,
                                        % get the first
                                        bestExpIdx = bestExpIdx(1);
                                    end
                                case 'PSOSVM'
                                    bestExpIdx = bestExpIdx(1);
                            end
                        end
                    end
                end
            end

            % mark the best index
            xlswrite(eachFileExcelPath, {'BEST EXPERIMENT'}, ...
                sprintf('%d classes', classNum(iClass)), ...
                sprintf('%s%d', ...
                getexcelcolumncode(length(headerEachExp)+1), bestExpIdx+1));
            
            %{
            % gather gBest fitness of the best experiment
            gBest = zeros(nIterations, 1);
            for iItr=1:nIterations
                gBest(iItr) = ...
                    ExperimentResult(bestExpIdx). ...
                    iterationResult(iItr+1).gBest.fitness;
            end

            % save graphics
            f = 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 excel
            switch method
                case '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));
            end
        end
    end

end

evaluatefitness、gbestupdate函数

function [modelArr, trainAccArr, testAccArr, timeArr, populationFitness, ...
    pBest] = evaluatefitness(PSOSettings, nFeatures, trainingData, ...
    testingData, populationPosition, pBest)

    modelArr(PSOSettings.nParticles).inputWeight = [];   % 定义modelArr结构体字段inputWeight、outputWeight
    modelArr(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.nParticles
        tic;
        % TRAINING
        maskedTrainingFeature = 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)));

        % TESTING
        maskedTestingFeature = 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 update
        ischanged = 0;
        % 满足以下任意一条件，更新变量
        if populationFitness(i, 1) > pBest(i).fitness
            ischanged = 1;
        elseif populationFitness(i, 1) == pBest(i).fitness
            if pBest(i).testingAccuracy < testAcc
                ischanged = 1;
            elseif pBest(i).trainingAccuracy < trainAcc
                ischanged = 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;
            end
        end

        if ischanged
            pBest(i).fitness = populationFitness(i, 1);
            pBest(i).position = populationPosition(i, :);
            pBest(i).trainingAccuracy = trainAcc;
            pBest(i).testingAccuracy = testAcc;
        end

        % end of pBest update
        modelArr(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.fitness
        found = 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 node
                        hn = zeros(length(found), 1);
                        for i=1:length(found)
                            hn(i, 1) = bintodec(populationPosition(found(i), ...
                                nFeatures+1:end));
                        end

                        found = found(hn == min(hn));

                        if length(found) > 1
                            found = found(1);
                        end
                    end
                end
            end
        end

        gBest.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.nParticles
        tic;

        % TRAINING
        maskedTrainingFeature = featuremasking(trainingData, ...
            populationPosition(i, 1:nFeatures)); % remove unselected features
        Model = trainSVM(maskedTrainingFeature, trainingData(:,end), 'RBF');
        trainAcc = testSVM(maskedTrainingFeature, trainingData(:,end), Model);

        % TESTING
        maskedTestingFeature = featuremasking(testingData, ...
            populationPosition(i, 1:nFeatures));     % remove unselected features
        testAcc = testSVM(maskedTestingFeature, testingData(:,end), Model);
        populationFitness(i, 1) = fitness(PSOSettings.Wa, PSOSettings.Wf, ...
            testAcc, populationPosition(i, 1:nFeatures));

        % pBest update
        ischanged = 0;

        if populationFitness(i, 1) > pBest(i).fitness
            ischanged = 1;
        elseif populationFitness(i, 1) == pBest(i).fitness
            if pBest(i).testingAccuracy < testAcc
                ischanged = 1;
            elseif pBest(i).trainingAccuracy < trainAcc
                ischanged = 1;
            elseif sum(pBest(i).position(1, 1:nFeatures)) > ...
                    sum(populationPosition(i, 1:nFeatures))
                ischanged = 1;
            end
        end

        if ischanged
            pBest(i).fitness = populationFitness(i, 1);
            pBest(i).position = populationPosition(i, :);
            pBest(i).trainingAccuracy = trainAcc;
            pBest(i).testingAccuracy = testAcc;
        end

        % end of pBest update

        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.fitness
        found = find(populationFitness == max(populationFitness));
        if length(found) > 1
            % if have the same gBest fitness value, get the max of testAcc
            found = found(testAccArr(found) == max(testAccArr(found)));
            if length(found) > 1
                % if have the same testAcc, get the max of trainAcc
                found = found(trainAccArr(found) == max(trainAccArr(found)));
                if length(found) > 1
                    % if have the same trainAcc, get the min of selected features
                    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 first
                        found = found(1);
                    end
                end
            end
        end

        gBest.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 features

    maskedFeature = zeros(size(feature, 1), sum(mask));
    j = 1;

    for i=1:sum(mask)
        if mask(1, i) == 1
           maskedFeature(:, j) = feature(:, i);
           j = j + 1;
        end
    end
end


function [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 accuracy

    if size(feature, 2) == 0
        fprintf('Someting went wrong, no feature selected.');
        return
    end

    % 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 weight
    outputWeight = 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);
    end

    trainAcc = sum(predictedClass == vec2ind(target')')/ ... % vec2ind：向量到索引，向量中每一列中的所有元素有且只能有一个为1，向量的列数即为索引矩阵的列数
        size(predictedOutput, 2) * 100;
    ELMModel.inputWeight = inputWeight;
    ELMModel.outputWeight = outputWeight;
end


function 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 output
    hiddenOutput = (ELMModel.inputWeight(:, 1:end-1) * feature')+ ...
        repmat(ELMModel.inputWeight(:, end), 1, size(feature, 1)); 
    % apply activation function on hidden output
    hiddenOutput = sigmoid(hiddenOutput);

    % linear combination of predicted output
    predictedOutput = ELMModel.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);
    end

    testAcc = sum(predictedClass == vec2ind(target')')/ ...
        size(predictedOutput, 2) * 100;    
end