基于matlab MFCC+GMM的安全事件声学检测系统

一、安全事件声学检测简介（附lunwen）

1 选题背景

公共安全问题是社会安全稳定所聚焦的话题之一。近年来，检测技术与监控自动化正深刻地改变着人们的生活。尤其在安防领域，闭路电视CCTV（Closed Circuit Television）、视频流分析、智能监控等新技术得到了广泛应用，大大提高了安防监控的管理效率。然而值得注意的是，基于视频流的监控手段不可避免地也具有一定的先天性缺漏，例如存在视野盲区、易受光照影响等问题，对于事件检测，还可能存在语义不明的问题，监控手段不够全面。纯视频手段在枪击、爆炸、暴恐袭击、人群恐慌等具有较强语义性的突发公共安全事件中，往往不如声学检测分析手段敏感和有效。声学事件检测主要是使用一些声学处理方法，刻画现场音频流的声学特征，再结合适当的分类器进行检测分类，从而实现对音频流中出现的声学事件进行检测分析。基于声学的公共安全事件检测在反恐、维稳、社会治安等多个领域具有广泛的使用价值和应用前景。本课题重点针对枪击与爆炸两类突发公共安全事件，对相应的声学检测方法进行了研究。

2 研究现状

对于枪声的研究聚焦在对于膛口波和弹道波的研究分析上。吴松林等深基于弹丸的空气动力学模型，深入分析了弹道激波的成因和理论波形；蒋灏等分析了小口径武器发射的膛口波和弹道激波，并设计了基于膛口激波的DOA模型对弹丸弹道轨迹进行估计；卢慧洋分析了弹道波和膛口波在枪声检测与定位中的作用，并设计了一套基于正三角形麦克风阵列的枪声定位与测距软硬件系统。
对于声信号处理和声学特征的研究，赵力等给出了常用的信号加窗成帧、端点检测以及常用声学特征的计算方法；韩纪庆等对声学事件检测技术与常用模型做了综述性介绍。徐大为等对比了基于不同声信号特征的端点检测方法，并分析了他们对噪声的抵制能力和运算实时性。
针对枪声的信号处理与声学事件检测研究中，蒋小为和张文等[7]通过低通滤波和谱减法针对膛口波进行去噪处理，在实验中得到了与理论波形高度相近的膛口波信号波形，如图1.22所示，并提出可以使用相关分析进行枪声检测。张克刚等人研究了基于短时能量分析对枪声信号进行端点检测的方法，并提出使用持续时间处理来剔除瞬时大能量噪声。张涛、张文、朱强强等人的研究中指出，可以采用MFCC作为目标片段的特征，用于进一步给分类器进行分类检测。
对于声学事件检测的分类器，Clavel等讨论了监控环境中的枪声检测，并通过PCA选择13维特征作为GMM模型的输入特征；刘力维等提出使用10阶中值滤波处理端点检测中的能量序列，并用GMM对目标片段的按MFCC特征进行分类；朱强强分析了Logo、FFS、Adaboost三种特征选择算法，用特征选择算法对时域特征、频域特征、感知域特征、基于自相关函数的特征等共计9个特征组成的特征全集进行特征选择，并最后输入到GMM中进行分类；Pimentel等人提出了通过分析聚类过程中的WSS指标来确定聚类算法中聚类中心数目的方法。
关于声学事件数据库，Fonseca等人所在的庞培法布拉大学（Universitat Pompeu Fabra, Barcelona）音乐技术研究小组为了解决目前数据驱动型（data-driven）声学计算研究所遇到的瓶颈和困难，发起了Freesound Datasets项目，并建立了一个基于众包（crowdsourcing）、规模宏大、音频种类较齐全的大型公开数据库Freesound；坦佩雷理工大学(Tampere University of Technology，TUT)信号处理学系Mesaro等人发起了事件检测挑战TUT Sound Events Challenge与声学场景检测挑战Acoustic Scene Classification Challenge，加速了基于声学的事件检测和场景分析的相关研究。

3 算法流程设计

基于matlab MFCC+GMM的安全事件声学检测系统

二、部分源代码

for ii = 1:24 % 7 8 wav分别为背景声和背景+枪声，21-24是爆炸声
    if 8 < ii &&  ii < 21
        continue;
    end
    % 取信号
    file_name = strcat("gun",num2str(ii));
    file_name = strcat(file_name,".wav");
    fprintf("reading %s...
",file_name);
    [y,fs] = audioread(file_name);   
    sz = size(y);
    gun = (y(:,1))"; % 单声道
    % 原信号
    figure(ii);
    p2 = abs(fft(gun)/length(gun));
    % size(gun)
    % size(1:length(gun)/2+1)
    gun_fft = p2(1:length(gun)/2+1);
    gun_fft(2:end-1) = 2*gun_fft(2:end-1);
    f = fs*(0:(length(gun)/2))/length(gun);
    subplot(3,2,1);plot(gun);xlabel("t / s");title("signal");
    subplot(3,2,2);plot(f,gun_fft);title("spectrum");xlabel("frequency / Hz");
    % 短时能量分析
    N = 300; % 窗宽（张克刚）
    inc = 100; % 帧移（张克刚）
    win = hamming(N);
    % frameout: num x N 
    % t: num x 1, centers of frames
    % energy: 1 x num
    [frameout,t,energy]=enframe(y,win,inc);
    t = t";
    % 自适应短时能量阈值分割
    %size(energy)
    threshold = min(energy)+0.2*(max(energy)-min(energy));
    processed_energy = energy;
    for i = 1:length(energy)
        processed_energy(i) = 0;
        if energy(i) >= threshold 
            processed_energy(i) = 1;
        end
        %fprintf("%d: %f > %f = %d
",ii,energy(i),threshold,processed_energy(i));
    end
    subplot(3,2,3); plot(energy,"b");title("energy");
    hold on;plot(threshold*ones(size(energy)),"g");
    subplot(3,2,5); plot(processed_energy);title("binarized energy")
    % 持续时间分析
    thr = 30; % 持续采样点
    cnt = 0;
    for i = 1:length(processed_energy)
        if processed_energy(i) == 1
            if cnt > 0
                cnt = cnt+1; %计数器累加
            elseif cnt == 0
                cnt = 1; %初始化计数器
            end
            if i == length(processed_energy) && cnt < thr
                processed_energy((i-cnt):i) = 0;
            end
        elseif processed_energy(i) == 0
            if cnt > 0
                if cnt < thr
                    processed_energy((i-cnt):i) = 0;
                end
            end
            cnt = 0;
        end
        %fprintf("%f, %f
",i,processed_energy(i));
    end
    
    subplot(3,2,3);hold on;plot(processed_energy*max(energy),"r");hold off;
    subplot(3,2,5);hold on;plot(processed_energy,"r");hold off;
end
function [f,t,eng,zcr]=enframe(x,win,inc)
%ENFRAME split signal up into (overlapping) frames: one per row. [F,T]=(X,WIN,INC)
%
%	F = ENFRAME(X,LEN) splits the vector X(:) up into
%	frames. Each frame is of length LEN and occupies
%	one row of the output matrix. The last few frames of X
%	will be ignored if its length is not divisible by LEN.
%	It is an error if X is shorter than LEN.
%
%	F = ENFRAME(X,LEN,INC) has frames beginning at increments of INC
%	The centre of frame I is X((I-1)*INC+(LEN+1)/2) for I=1,2,...
%	The number of frames is fix((length(X)-LEN+INC)/INC)
%
%	F = ENFRAME(X,WINDOW) or ENFRAME(X,WINDOW,INC) multiplies
%	each frame by WINDOW(:)
%
%   The second output argument, T, gives the time in samples at the centre
%   of each frame. T=i corresponds to the time of sample X(i). 
%
nx=length(x);
nwin=length(win);
if (nwin == 1)
   len = win;
else
   len = nwin;
end
if (nargin < 3)
   inc = len;
end
len = nwin;
nf = fix((nx-len+inc)/inc);
f=zeros(nf,len);
indf= inc*(0:(nf-1)).";
inds = (1:len);
f(:) = x(indf(:,ones(1,len))+inds(ones(nf,1),:));
if (nwin > 1)
    w = win(:)";
    f = f .* w(ones(nf,1),:);
end
t = floor((1+len)/2)+indf;
%fprintf("size of f
");
szf = size(f);
% ff = f(:).*f(:);
for i = 1:szf(1)
    %ff = f(i,:).*f(i,:)
%     ff = abs(f(i,:));
%     eng(i) = sum(ff);
    eng(i) = 0;
    zcr(i) = 0;
    for j = 1:szf(2)
        eng(i) = eng(i)+abs(f(i,j));
        if j+1 <= szf(2)
            zcr(i) = zcr(i)+abs(sign(f(i,j+1))-sign(f(i,j)));
        end
    end
    zcr(i) = 0.5*zcr(i);
end