下载数据

这个例子适用于声学场景和事件检测和分类(DCASE) 2022挑战赛的第二个任务[1]．该示例使用来自故障工业机器调查和检查声音数据集的公共数据集的子集[２]训练和评估自动编码器。它实现了自编码器基线系统的预处理步骤和网络设计思想[1]并提出网络在[２]并使用在[1]对测试结果进行分析。

下载数据集的子集[２]其中包含4种不同风扇类型的录音文件，以ID号标记。每种风扇类型都有正常和异常记录。这些文件包含一个以16 kHz采样的频道，长度为10秒。这些样本是有背景噪声的风扇运行记录，信噪比为6 dB。关于数据收集过程的详细说明见[２]．

dataFolder = tempdir;dataset = fullfile(数据文件夹，“fan6db”）;金宝appsupportFileLoc =“mimii / mono / fan6db.zip”；downloadFolder = matlab.internal.examples.download金宝appSupportFile(“音频”,金宝app supportFileLoc);解压缩(downloadFolder dataFolder)

调查数据

为了简单地检查数据集以及正常和异常录音之间的差异，从ID 00风扇数据集中选择每种类型的录音，并在扬声器上播放前两秒。

[normalSample,fs] = audioread(fullfile(dataset, fs)“id_00”，“normal_00”，“00000000. wav”));abnormalSample = audioread(fullfile(dataset，“id_00”，“abnormal_00”，“00000000. wav”));numSamples = 10*fs;pause(3) sound(abnormalSample(1:numSamples/5)，fs)

两个录音都是由一个音调主导的，这个音调在异常样本中明显更高。

数据进行预处理

您可以选择设置加速旗帜真正的以减少示例中使用的数据集的大小。如果你把这个设为真正的您可以快速地验证脚本是否按预期运行，但结果会有偏差。

加速=假；

将数据集分成两个audioDatastore(音频工具箱)对象，一个是正常样本，一个是异常样本。由于自动编码器只对正常样本进行训练，因此将异常样本包含在测试集中。

ads = audioDatastore(数据集，.．.IncludeSubfolders = true,.．.LabelSource =“foldernames”，.．.FileExtensions =“wav”）;normalLabels = categorical([“normal_00”，“normal_02”，“normal_04”，“normal_06”]);abnormalLabels = categorical([“abnormal_00”，“abnormal_02”，“abnormal_04”，“abnormal_06”]);isNormal = ismember(ads.Labels,normalLabels);isAbnormal = ~isNormal;adsNormal =子集(广告，isNormal);adsTestAbnormal =子集(ads,isAbnormal);提高(3);如果加速c = cvpartition(adsTestAbnormal.Labels,kFold=8,Stratify=true);adsTestAbnormal =子集(adsTestAbnormal,c.test(1));结束

将正常样本分为训练集、验证集和测试集，按ID号分层。然后将正常测试集与异常样本拼接成完整测试集。

c = cvpartition(adsNormal.Labels,kFold=8,Stratify=true);如果speedUp trainInd = c.test(3);其他的trainInd = ~boolean(c.test(1)+c.test(2));结束valInd = c.test(1);testInd = c.test(2);adsTrain =子集(adsNormal,trainInd);adsVal = adsNormal,valInd;adsTestNormal =子集(adsNormal,testInd);

通过应用帧长64 ms、跳长32 ms的STFT对每个数据存储进行转换，找到128个频段的log-mel能量，然后将这些帧连接到重叠的连续5组中，形成一个上下文窗口。使用log-mel能量作为音频深度学习任务的输入是很常见的，因为它们代表了音调的频谱，其规模类似于人类感知声音的方式。可视化log-mel谱图的两个片段播放之前使用plotLogMelSpect金宝app支持功能。

plotLogMelSpect (normalSample abnormalSample);

使用processData金宝app支持数据转换功能。

tdsTrain = transform(adsTrain，@processData);tdsVal = transform(adsVal，@processData);tdsTestNormal = transform(adsTestNormal，@processData);tdsTestAbnormal = transform(adsTestAbnormal，@processData);

将数据读入数组，其中每一列代表一个输入示例。如果已启用并行计算工具箱™，则并行执行此操作。然后将正常测试集和异常数据集合并为完整测试集，并对样本进行标记。

trainingData = readall(tdsTrain,UseParallel=canUseParallelPool);valData = readall(tdsVal,UseParallel=canUseParallelPool);normalTestData = readall(tdsTestNormal,UseParallel=canUseParallelPool);abnormalTestData = readall(tdsTestAbnormal,UseParallel=canUseParallelPool);testLabels = categorical([zero (length(adsTestNormal.Labels)，1);ones(length(adsTestNormal.Labels)，1)]，[0,1]，[“正常”，“不正常”]);testData = [normalTestData;abnormalTestData];

网络体系结构

编码器部分由2个完全连接的层组成，输出大小为128。瓶颈层将网络限制为原始640维输入的8维表示。解码器部分反映编码器体系结构，因为输入被重构并馈送到输出层。采用半均方误差作为损失函数对网络进行训练，量化重构误差。

layers = [featureInputLayer(640) fullyConnectedLayer(128,Name= .“Encoder1”) batchNormalizationLayer reluLayer fullyConnectedLayer(128,Name=“Encoder2”) batchNormalizationLayer reluLayer fullyConnectedLayer(8,Name=“瓶颈”) batchNormalizationLayer reluLayer fullyConnectedLayer(128,Name=“Decoder1”) batchNormalizationLayer reluLayer fullyConnectedLayer(128,Name=“Decoder2”batchNormalizationLayer reluLayer fullyConnectedLayer(640,Name=“输出”) regressionLayer);

列车网络的

使用ADAM优化器训练网络40个epoch。对每个纪元的小批进行洗牌，并设置ExecutionEnvironment字段“汽车”如果可用，则使用GPU而不是CPU。如果使用内存有限的GPU，则可能需要降低miniBatchSize字段。通过实验确定了训练参数的设置，以优化收敛速度。这可能需要10-15分钟，具体取决于您的硬件。

batchSize = length(trainingData)/2;如果2*batchSize;结束选项= trainingOptions(“亚当”，.．.MaxEpochs = 30,.．.InitialLearnRate = 1飞行,.．.LearnRateSchedule =“分段”，.．.LearnRateDropPeriod = 5,.．.LearnRateDropFactor = 7,.．.GradientDecayFactor = 8,.．.miniBatchSize = batchSize,.．.洗牌=“every-epoch”，.．.ExecutionEnvironment =“汽车”，.．.ValidationData = {valData, valData},.．.ValidationFrequency = 2,.．.Verbose = 0,.．.情节=“训练进步”）;

trainingData是网络试图用低维编码约束在自身上回归训练数据时的输入和目标输出。您的结果应该类似于下面的训练图。

[net,info] = trainNetwork(trainingData,trainingData,layers,options);

评估性能

对于网络的每个输入，自动编码器输出一个尝试重构。然而，每个网络输入只是一个更大的音频样本中的一个上下文窗口。对于这些网络输入中的每一个，误差被定义为原始输入与网络输出之间差值的L-2范数的平方。为了计算每个完整音频样本的决策度量，将与该音频样本相关的每个上下文窗口的误差加在一起，这个和除以网络输入维度和每个音频样本的上下文组的数量的乘积。对于音频样本 $\mathit{X}$，表示决策函数度量 $\mathit{一个}（\mathit{X}）$和定义 $\mathit{一个}（\mathit{X}）＝{\sum }_{\mathit{我}＝1}^{\mathit{n}}\frac{{为\mathit{f}（{\mathit{x}}_{\mathit{我}}）-{\mathit{x}}_{\mathit{我}}为}^2}{\mathit{n}＊\mathrm{昏暗的}（{\mathit{x}}_{\mathit{我}}）}$在哪里 $\mathit{n}$是每个样本的上下文组数， ${\mathit{x}}_{\mathit{我}}$是 ${\mathit{我}}^{\mathrm{th}}$从 $\mathit{X}$, $\mathit{f}（{\mathit{x}}_{\mathit{我}}）$是网络输出为 ${\mathit{x}}_{\mathit{我}}$． $\mathit{一个}（\mathit{X}）$表示与音频样本相关的所有上下文Windows的每个向量组件之间的均方重建误差 $\mathit{X}$．

对于每个输入 $\mathit{X}$， $\mathit{一个}（\mathit{X}）$也可以解释为相对衡量网络的信心 $\mathit{X}$异常，数值越高表示置信度越大。要部署此模型并对新数据进行预测，必须在的值上选择决策边界 $\mathit{一个}$区分积极和消极的预测。模型 $\mathit{一个}（\mathit{X}）$对于正态样本的伽马分布。伽玛分布通常用于自编码器重建误差的建模，因为误差通常是向右倾斜的，有一个沉重的尾部，这是伽玛分布的自然形状。在本例中，选择决策边界作为与预期假阳性率(FPR)对应的点。p= 0.1。该决策边界试图捕获所有真正异常的样本，同时容忍10%的正常样本将被错误地预测为异常。的特定值p以满足您的个别系统限制。

P = .1;

计算的值 $\mathit{一个}$然后把它们存储在变量中A_train使用helper函数getScore．然后求解伽玛分布参数的最大似然估计，从伽玛累积分布逆函数中选择截断点，用的直方图绘制拟合分布 $\mathit{一个}$使用getCutoffhelper函数。

trainRecons = predict(net,trainingData,MiniBatchSize=length(trainingData));A_train = getScore(trainingData,trainRecons);cutoff = getCutoff(A_train,p);

验证该截断点大致对应于训练集上的FPR为0.1:

sum(A_train > cutoff) / length(A_train)

Ans = 0.1086

测试该系统的分类精度与选定的截止点上的坚持测试集。

testRecons =预测(net,testData,MiniBatchSize=长度(testData));A_test = getScore(testData,testRecons);testPreds = categorical(A_test > cutoff，[false,true]，[“正常”，“不正常”])。';figure cm = confusichart (testLabels,testPreds);厘米。RowSummary =“row-normalized”；

使用这个截止点，该模型以0.125的FPR为代价获得了0.647的真阳性率(TPR)。

为了评估网络在决策边界范围内的准确性，通过接收器工作特征曲线(AUC)下的面积来测量测试集上的整体性能。同时使用完整AUC和部分AUC (pAUC)来分析网络性能。pAUC是子域上的AUC，其中FPR在区间上 $［0，\mathit{p}］$除以区间内可能的最大面积，也就是p．考虑pAUC是很重要的，因为异常检测系统需要能够在保持FPR最小化的同时实现高TPR，因为频繁出现假警报的系统是不值得信任和不可用的。函数计算AUCperfcurve功能从统计和机器学习工具箱™。

[X,Y,T,AUC] = perfcurve(testLabels,A_test,categorical(“不正常”));[~，cutoffIdx] = min(abs(T-cutoff));图绘制(X, Y);包含(“玻璃钢”）;ylabel (“TPR”）;标题(“测试集ROC曲线”）;持有在情节(X (cutoffIdx), Y (cutoffIdx),的r *）;持有从传奇(“ROC曲线”，“截止决策点”）;网格在

AUC

AUC =单0.8957

为了计算pAUC，在FPR域的前十分之一近似曲线下的面积使用trapz．作为参考，随机分类器的pAUC的期望值为0.05。

pX = X(X <= p);pY = Y(X <= p);pAUC = trapz(pX,pY)/p

pAUC =单0.5161

该网络很好地分离了正常和异常测试样本，并能够跨多个风扇id学习单一编码。通过直方图可见正常组和异常组重建误差的差异。

图保存在linspace(min(A_test)， 1100);直方图(A_test(testLabels == categorical(“正常”)),边缘,标准化=“概率”）;直方图(A_test(testLabels == categorical(“不正常”)),边缘,标准化=“概率”）;ylabel (“样本概率”）;包含(“重建误差(A)”）;传奇(“正常”，“不正常”）;

尽管存在一些重叠，但异常样本的重建误差分布进一步向右偏移，并且包含比正常样本上重建误差分布更重的尾部。

最后，评估模型在每个风扇ID上的性能，以揭示风扇类型之间的任何不平衡，并检查模型是否能够对所有ID进行普遍良好的预测。

id = [0;2;4;6];AUCs = 0 (4,1);pAUCs = AUCs;A_testNormal = A_test(1:sum(testInd));A_testAbnormal = A_test(sum(testInd)+1:end);为i = 1:4 normalMask = adsTestNormal。标签== normalLabels(i);abnormalMask = adsTestAbnormal。标签== abnormalLabels(i);A_testByID = [A_testNormal(normalMask) A_testAbnormal(abnormalMask)];testLabelsByID = [adsTestNormal.Labels(normalMask); adsTestNormal.Labels(abnormalMask)];[X_ID,Y_ID,T_ID,AUC_ID] = perfcurve(testLabelsByID,A_testByID,abnormalLabels(i));AUCs(i) = AUC_ID;pX_ID = X_ID(X_ID <= p);pY_ID = Y_ID(X_ID <= p);pAUCs(i) = trapz(pX_ID,pY_ID)/p;结束disp(表(id、auc pAUCs));

IDs AUCs pAUCs ___ _______ _______ 0 0.72147 0.19134 2 0.97739 0.73454 4 0.8761 0.42547 6 0.96017 0.6742

结果表明，不同风机类型的模型性能有显著差异。这一结果值得注意，因为与提交的最佳DCASE挑战相比，该网络相对较小且简单[3]．为了更好地概括风扇类型和不同领域，需要一个更复杂的模型。然而，如果您确切地知道要部署异常检测器的风扇类型，那么像本例中这样的非常轻量级的模型可能就足够了。

金宝app支持功能

函数plotLogMelSpect (normalSample abnormalSample)%PLOTLOGMELSPECT绘制正常和非正常的log-mel谱图% plotLogMelSpect(normalSample,abnormalSample)绘制log-mel两个输入并列的%谱图，参数一致与数据预处理转换用于准备信号%输入自动编码器。F =数字;f.位置(3)= 900;样本= {normalSample,abnormalSample};Fs = 16e3;winDur = 64e-3;winLen = windir * fs;numMelBands = 128;tiledlayout(1、2)为I = 1:2 nexttile x = samples{I};melSpectrogram (x, fs,窗口=汉明(winLen,“周期”)、FFTLength = winLen OverlapLength = winLen / 2, NumBands = numMelBands);xticks (1:10);xticklabels (string (1:10));colormap (“喷气机”）;如果I == 2 cbar = colorbar;cbar.Label.String =“权力(dB)”；标题(“异常Log-Mel谱图”）;ylabel ([]);其他的colorbar从标题(“正对数-梅尔谱图”）;结束结束结束函数features = processData(x)PROCESSDATA将音频文件输入x转换为自动编码器网络%输入格式% features = processData(x)取音频数据x的STFT，进行转换将STFT转换为log-mel谱图，然后构造上下文%组连续mel-谱图帧。函数返回%是一个numContextGroupsPerSample-by-contextGroupSize矩阵。为%此数据集，numContextGroupsPerSample = 309和contextGroupSize =% 640 = 128*5(因为每帧有128个MEL波段，5帧是为每个上下文组连接的%)Fs = 16e3;winDur = 64e-3;winLen = windir *fs;numMelBands = 128;afe = audioFeatureExtractor(.．.窗口=汉明(winLen,“周期”)，.．.FFTLength = winLen,.．.OverlapLength = winLen / 2,.．.SampleRate = fs,.．.melSpectrum = true);setExtractorParameters (afe“melSpectrum”numBands = numMelBands);%零垫numSamples =长度(x)numPad = winLen - mod(numSamples,winLen);numToPadFront =地板(numPad/2);numToPadBack = cell (numPad/2);xpadding = [0 (numToPadBack,1,like=x);x; 0 (numToPadBack,1,like=x)];%提取features = extract(afe, xpadding);Features = {log10(Features)};features = cellfun(@groupSTFT,features,UniformOutput=false);Features = vertcat(Features {:});结束函数groups = groupSTFT(x)%GROUPSTFT将STFT x转换为大小为5的上下文组% groups = groupSTFT(x)通过对每个STFT帧进行分组来转换STFT x%与以下4帧组成大小为5的上下文组。这%从每个音频样本中创建多个网络输入，每个输入大小相同% contextLen*numMelBands = 5*128 = 640。每个上下文组都是的目的，将%作为单独的640维向量% autoencoder。contextLen = 5;numMelBands = 128;X_flat =重塑(x'，1，[]);组= buffer(x_flat,contextLen*numMelBands,numMelBands*(contextLen-1)，“nodelay”)”;结束函数A = getScore(数据，preds)返回数据中每个样本的重构误差% A = getScore(data,preds)为样本集中的每个X返回A(X)%转换为网络输入数据。Err = sum((preds-data).^2,2);numSTFTFrames = 313;contexttwin = 5;numMelFilters = 128;numContextGroupsPerSample = numstftframes - contexttwin +1;numSamples = length(err)/numContextGroupsPerSample;A_total =重塑(err，[numContextGroupsPerSample,numSamples]);每一列包含一个样本中所有上下文组的重构误差A = sum(A_total)/(numMelFilters* contexttwin *numSTFTFrames);每个条目是每个样本的一个重构误差结束函数cutoff = getCutoff(A,p)%GETCUTOFF拟合gamma分布到a，并返回截断值为1-p的逆cdf% cutoff = getCutoff(A,p)拟合gamma分布的重建%误差数组A，求解截止点为逆伽马cdf%在1-p处评估，并将拟合分布与A的%直方图和计算的截断点。预期A为% one by numSamples数组，其中numSamples是音频样本的个数%用于计算A的重构误差值。gammaParams = gamfit(A);a = gammaParams(1);b = gammaParams(2);截止值= gaminv(1-p,a,b);图ax1 = subplot(4,1,1:3);柱状图(一个);xticks ([]);标题(A与拟合伽马区直方图。PDF）;ylTop = ylabel(“数”）;参照线(截止,”——“线宽= 2,标签=“截止”LabelOrientation =“水平”LabelVerticalAlignment =“中间”）;Ax2 = subplot(4,1,4);t = linspace(0, max(A)， 1000);Y = gampdf(t,a,b);情节(t、y);参照线(截止,”——“、线宽= 2);ylBottom = ylabel(“\伽马密度”）;yticks ([]);linkaxes ([ax₁ax2],“x”）;ylBottom.Position(1) = ylTop.Position(1);包含(“重建误差(A)”）;xlim([0。4]);ylBottom.Position(1) = ylTop.Position(1);结束

参考文献

[1]“应用领域泛化技术的机器状态监测的无监督异常声音检测”DCASE 2022．(在线)。可用:https://dcase.community/challenge2022/task-unsupervised-anomalous-sound-detection-for-machine-condition-monitoring。(访问:08 - 6月- 2022)。

[2]“Purohit, Harsh, Tanabe, Ryo, Ichige, Kenji, Endo, Takashi，二阶堂，Yuki, Suefusa, Kaori， &川口，Yohei。(2019)。MIMII数据集:工业机械故障调查与检测声音数据集(公开1.0)[数据集]。第四届声学场景与事件检测与分类研讨会(DCASE 2019研讨会)，美国纽约。Zenodo。https://doi.org/10.5281/zenodo.3384388．数据集采用知识共享署名- sharealike 4.0国际许可https://creativecommons.org/licenses/by-sa/4.0/

［3]“机器状态监测中异常声音的无监督检测，”DCASE 2020．(在线)。可用:https://dcase.community/challenge2020/task-unsupervised-detection-of-anomalous-sounds-results # Giri2020。(访问:08 - 6月- 2022)。