BP神经网络项目一、对声音的特征提取(代码分析)

BP神经网络项目一、对声音的特征提取(代码分析)

项目来源：b站链接 BP神经网络原理及编程实现_哔哩哔哩 _bilibili

1.数据读取+预处理

1.1加载文件+文件合成

%% 清空环境变量
clear,clc;

%% 训练数据预测数据提取及归一化

%下载四类语音信号
load data1 c1
load data2 c2
load data3 c3
load data4 c4

%四个特征信号矩阵合成一个矩阵
data(1:500,:)=c1(1:500,:);
data(501:1000,:)=c2(1:500,:);
data(1001:1500,:)=c3(1:500,:);
data(1501:2000,:)=c4(1:500,:);

将题目所给四个文件合成data矩阵，大小2000*25

1.2 划分标签与特征矩阵

%输入输出数据
input=data(:,2:25);%取出data中2到25列数据构成新矩阵input
output1 =data(:,1);%取出data中第一列数据构成矩阵output1

%把输出从1维变成4维
for i=1:2000
    switch output1(i)
        case 1
            output(i,:)=[1 0 0 0];
        case 2
            output(i,:)=[0 1 0 0];
        case 3
            output(i,:)=[0 0 1 0];
        case 4
            output(i,:)=[0 0 0 1];
    end
end

• 创建特征矩阵input，大小为2000*24

• 创建样本标签矩阵output1，大小为2000*1

• 拓展样本标签矩阵output，使其适应分四类的任务需求，大小为2000*4

1.3划分训练集和测试集

%从1到2000间随机排序
k=rand(1,2000);%生成0-1之中2000个随机数
[m,n]=sort(k);%对k进行从从小到大的排序，序结果放入m（1*2000）向量，n（1*2000）为原数据对应的索引号

%随机提取1500个样本为训练样本，500个样本为预测样本
input_train=input(n(1:1500),:)';   % M' 是求M矩阵的转置
output_train=output(n(1:1500),:)';
input_test=input(n(1501:2000),:)';
output_test=output(n(1501:2000),:)';

• 利用rand和sort函数随机分出训练集和测试集

• 训练集：input_train(大小为：24*1500)、output_train(大小为：4*1500)

• 测试集：input_test(大小为：24*500)、output_test(大小为：4*500)

1.4数据归一化+确定各层的维度

%输入数据归一化
[inputn,inputps]=mapminmax(input_train);%归一化到[-1,1]之间，inputps用来作下一次同样的归一化

%%网络结构初始化，设置节点输入层，隐藏层，输出层
innum=24;
midnum=25;%选择多少个节点才比较好？
outnum=4;

• 将数据归一化到[-1,1]区间内,inputn矩阵大小：24*1500

• Input_layer : 24D ; Hidden_layer : 25D ; output_layer : 4D

1.5 权值初始化+学习率的确定

%权值初始化 %rands在[-1,1]范围内生成随机数
w1=rands(midnum,innum);%随机给定隐藏层和是输入层间的初始神经元权重   W1= net. iw{1, 1}; 
b1=rands(midnum,1);%中间各层神经元阈值  B1 = net.b{1};
w2=rands(midnum,outnum);%中间层到输出层的权值  W2 = net.lw{2,1};
b2=rands(outnum,1);%输出层各神经元阈值  B2 = net. b{2}

%权重，偏差重新赋值
w2_1=w2;w2_2=w2_1;%把w2中的值分别配给w1w2
w1_1=w1;w1_2=w1_1;
b1_1=b1;b1_2=b1_1;
b2_1=b2;b2_2=b2_1;

%学习率
xite=0.1   %权值阈值更新

• 用rands函数随机生成随机权值和偏置矩阵w1,w2,b1,b1

• 生成矩阵大小：w1(大小：25*24) ; b1(大小：25*1) ; w2(大小：25*4) ; b2(大小：4*1)

• 学习率：xite = 0.1

2.模型训练

2.0总代码

for ii=1:10
    for i=1:1:1500
       %% （2）网络预测输出 
        x=inputn(:,i);
        % (3)隐含层输出
        for j=1:1:midnum
            I(j)=inputn(:,i)'*w1(j,:)'+b1(j); 
            Iout(j)=1/(1+exp(-I(j)));  % 激活函数:Sigmoid函数
        end
        % (4)输出层输出
        yn=w2'*Iout'+b2;
        
       %%(5) 权值阀值修正
        %计算误差
        e=output_train(:,i)-yn;     
        
        %计算权值变化率
        dw2=e*Iout;
        db2=e';
        
        for j=1:1:midnum
            S=1/(1+exp(-I(j)));
            FI(j)=S*(1-S);   %% FI得出的是Sigmoid函数的导数
        end      
        for k=1:1:innum
            for j=1:1:midnum
                dw1(k,j)=FI(j)*x(k)*(e(1)*w2(j,1)+e(2)*w2(j,2)+e(3)*w2(j,3)+e(4)*w2(j,4));
                db1(j)=FI(j)*(e(1)*w2(j,1)+e(2)*w2(j,2)+e(3)*w2(j,3)+e(4)*w2(j,4));
            end
        end
          %（6）权值阈值更新，学习率在这用 
        w1=w1_1+xite*dw1';
        b1=b1_1+xite*db1';
        w2=w2_1+xite*dw2';
        b2=b2_1+xite*db2';
        
        w1_2=w1_1;w1_1=w1;
        w2_2=w2_1;w2_1=w2;
        b1_2=b1_1;b1_1=b1;
        b2_2=b2_1;b2_1=b2;
    end
end

2.1最外层循环：循环轮数

for ii=1:10
    E(ii)=0;
%   (代码省略)
    E(ii)=E(ii)+sum(abs(e));
end

• 最外层循环为训练轮数，用ii记录，此处总共进行10轮训练

• E(ii)用于记录每一轮的训练的总误差，最后可以做误差图，研究最合适的训练次数

2.2 次外层循环：对所有样本进行遍历

for i=1:1:1500
%	(代码省略)
end

• i为样本标号，从1~1500

2.3 隐藏层输入和输出

 % 隐含层输出
for j=1:1:midnum
    I(j)=inputn(:,i)'*w1(j,:)'+b1(j); 
    Iout(j)=1/(1+exp(-I(j)));  % 激活函数:Sigmoid函数
end

% 输出层输出
yn=w2'*Iout'+b2;

• 矩阵相乘图解

• 数据大小：inputn的大小24*1500，w1大小25*24，b1大小25*1

• 循环参数j用来标记隐藏层结点的序号j

• 第j个结点的输入为：inputn的第i列的转置(1*24) × w1的第j行的转置(24*1)+b1的第j个元素

• 得到I(j)为第j个结点的输入值，大小为：1*1

• 结合激活函数Sigmoid函数得到Hidden_layer输出值,Iout(j)为第j个结点的输出值，大小为1*1

• 上述循环得到的结果为：隐藏层输入I矩阵，大小1*25，隐藏层输出矩阵，大小1*25

• 由公式得到输出层的输出（w2大小25*4，b2大小为4*1）yn

• 第i个样本的输出yn为：权值向量w2的转置(4*25) × 隐藏层输出向量Iout的转置(25*1),得到的输出向量yn，大小为4*1

2.4 计算误差+L关于w2,b2的偏导数(为梯度下降做准备)

 %% 权值阀值修正
%计算误差
e=output_train(:,i)-yn;     
        
%计算权值变化率
dw2=e*Iout;
db2=e';

• output_train为output经过随机分组后取转置的矩阵，大小4*1500，因而output_train( : , i )大小4*1

• 误差e计算：output_train( : , i )(大小4*1) - yn(大小4*1) 得到误差向量e大小4*1

$$1.损失函数对{\omega }_2的偏导:\frac{\partial L}{\partial {\omega }_2}=-E\ast Hidden\_Out$$

$$2.损失函数对b_2的偏导:\frac{\partial L}{\partial b_2}=-E$$

• 根据上述推导：可以得到权值w2和偏置b2的负变化率dw2，db2（由于是负变化率，因此做梯度下降法的之后可以直接w2 += η(学习率)*dw2）

• dw2 = e(大小4*1) × Iout(1*25) 得到dw2大小为4*25与w2对称√

• db2 = e转置（大小1*4）得到db2大小为4*1与b2对称√

• 最后在进行梯度下降时dw2和db2都需要求转置，因而转置之后dw2与db2分别与w2和b2的大小形状相同

2.5 计算损失函数L关于w1和b1的偏导数

for j=1:1:midnum
    S=1/(1+exp(-I(j)));
    FI(j)=S*(1-S);   %% FI得出的是Sigmoid函数的导数
end  

• 由激活函数得到每一个结点的S，其实就是上面的Iout(j):隐藏层第j个结点的输出

• Sigmoid函数求导公式如下↓：

$$S(x)=\frac{1}{1+e^{-x}}、S^{\prime }(x)=\frac{e^{-x}}{(1+e^{-x}{)}^2}=S(x)\ast (1-S(x))$$

• I的大小1*25，因此数据I(j)的大小为1*1

• 根据上述公式分析得：隐藏层第j个结点的输出的导数值FI(j) = S*(1-S)

• 最终得到的导数值矩阵为FI，大小为1*25

x=inputn(:,i);

for k=1:1:innum
    for j=1:1:midnum
       dw1(k,j)=FI(j)*x(k)*(e(1)*w2(j,1) + e(2)*w2(j,2) + e(3)*w2(j,3)        + e(4)*w2(j,4));
       db1(j)=FI(j)*(e(1)*w2(j,1) + e(2)*w2(j,2) + e(3)*w2(j,3) + e(4)        *w2(j,4));
    end
end

$$3.损失函数对{\omega }_1的偏导：\frac{\partial L}{\partial {\omega }_1}=-E\ast Hidden\_Out\ast x\ast S^{\prime }$$

$$4.损失函数对b_1的偏导：\frac{\partial L}{\partial b_1}=-E\ast Hidden\_Out\ast S^{\prime }$$

• 使用数据：x(输入层的数据为inputn的第i列)(大小24*1)，误差e(大小4*1)，w2(大小25*4)，w1(大小25*24)，b1（大小25*1）

• 最终求出的dw1矩阵和db1矩阵应分别与w1和b1对称，即dw1的大小为24*25，db1的大小为1*25

• 求w1矩阵第k行(第k行代表隐藏层第k个输入&特征)，第j列(第j列代表第j个结点)元素为：第j个结点的输出导数S’(j)（大小1*1） * 第k个输入x(k)（大小1*1） * 总的误差e的转置（大小14） * 第j个结点的w2转置，即w2‘( j , : )（大小41），最终得到dw1(k,j)的大小为1*1

• 求b1矩阵第j个(第j个代表第j个结点)元素为：第j个结点的输出导数S’(j)（大小1*1） * 总的误差e的转置（大小14） * 第j个结点的w2转置，即w2‘( j , : )（大小41），最终得到db1(k,j)的大小为1*1

2.6梯度下降法(Gradient Decent)更新权值w1，w2和偏置b1，b2

%权值阈值更新，学习率在这用 
w1=w1+xite*dw1';
b1=b1+xite*db1';
w2=w2+xite*dw2';
b2=b2+xite*db2';

• 由于前面求得的是负变化率（即负导数），因此此处的梯度下降法使用+而不是-

3.模型测试

3.1用同样的映射方法对测试集进行归一化

%% 语音特征信号分类
inputn_test=mapminmax('apply',input_test,inputps);

• 此处的映射方法来源于:[inputn,inputps]=mapminmax(input_train);

3.2 开始测试集测试

for ii=1:1
    for i=1:500
        %隐含层输出
        for j=1:1:midnum
            I(j)=inputn_test(:,i)'*w1(j,:)'+b1(j);
            Iout(j)=1/(1+exp(-I(j)));
        end
        
        fore(:,i)=w2'*Iout'+b2;
    end
end

• 使用的数据：w2大小：25*4，w1大小：25*24，b1大小：25*1，b2大小：，inputn_test大小：24*500

• ii用于标记测试的轮数：总共进行1轮测试

• i为样本标号：从1_{500，j为隐藏层结点标号：从1}25(midnum)

• 第k个隐藏层结点的输入I(j) = inputn_test的第j列的转置(大小1*24)×w1第j行的转置(大小24*1)+b1的j位置元素(大小1*1).最终得到隐藏层输入向量，大小为1*25

• 隐藏层第j个结点的输出值Iout(j)为Sigmoid(I(j)),最终得到隐藏层输出向量，大小为1*25

• 计算得预测值矩阵fore的第i列(第i个样本的预测值)为w2转置(大小4*25)×Iout转置(大小25*1)+b2(大小4*1)，最终得到全样本预测值矩阵大小为4*500

4.测试结果分析

4.1 找出每个预测值属于哪一类 + 获得预测误差

%% 结果分析
%（1）根据网络输出找出数据属于哪类
for i=1:500
    output_fore(i)=find(fore(:,i)==max(fore(:,i)));
end

%（2）BP网络预测误差
error=output_fore-output1(n(1501:2000))';

• find(条件A) 返回满足条件A的元素的下标

• output_fore向量记录预测的输出值，输出值为预测矩阵列的每一列的最大值

• 用预测值 - 实际值得到误差向量error，大小为1*500

4.2作出误差图+预测语音和实际语音种类分类图

%画出预测语音种类和实际语音种类的分类图
figure(2)
plot(output_fore,'r')
hold on
plot(output1(n(1501:2000))','b')
legend('预测语音类别','实际语音类别')

%画出误差图
figure(3)
plot(error)
title('BP网络分类误差','fontsize',12)
xlabel('语音信号','fontsize',12)
ylabel('分类误差','fontsize',12)

• fontsize是用来调整字号的

4.3统计每一类的错误

4.3.1 统计各类别错误的个数

k=zeros(1,4);  
%找出判断错误的分类属于哪一类
for i=1:500
    if error(i)~=0               %~表示非也就是error不等于0是
        [b,c]=max(output_test(:,i));
        switch c
            case 1 
                k(1)=k(1)+1;
            case 2 
                k(2)=k(2)+1;
            case 3 
                k(3)=k(3)+1;
            case 4 
                k(4)=k(4)+1;
        end
    end
end

• 创建全0矩阵k，用于承接四种类型的错误数

• i为测试样本的序号，判断error(i)~0的大小得到第i个样本分类有误

• 用max()得到分类有误的样本属于的类别，进而统计每个类别的正确率

• max(M)函数返回[M1,M2]:M1为最大值向量(默认取列的最大值)，M2为最大值的索引值向量

4.3.2 统计测试集中每一个的样本数目 + 正确率计算

%找出每类的个体和
kk=zeros(1,4);
for i=1:500
    [b,c]=max(output_test(:,i));
    switch c
        case 1
            kk(1)=kk(1)+1;
        case 2
            kk(2)=kk(2)+1;
        case 3
            kk(3)=kk(3)+1;
        case 4
            kk(4)=kk(4)+1;
    end
end

• 原理同上，kk向量用于统计测试集中属于各类的样本数

%正确率
rightridio=(kk-k)./kk

于的类别，进而统计每个类别的正确率

• max(M)函数返回[M1,M2]:M1为最大值向量(默认取列的最大值)，M2为最大值的索引值向量

4.3.2 统计测试集中每一个的样本数目 + 正确率计算

%找出每类的个体和
kk=zeros(1,4);
for i=1:500
    [b,c]=max(output_test(:,i));
    switch c
        case 1
            kk(1)=kk(1)+1;
        case 2
            kk(2)=kk(2)+1;
        case 3
            kk(3)=kk(3)+1;
        case 4
            kk(4)=kk(4)+1;
    end
end

• 原理同上，kk向量用于统计测试集中属于各类的样本数

%正确率
rightridio=(kk-k)./kk