Автоматизированная система диагностики дефектов в конструкциях электронных средств на основе акустических сигналов

 

4. ПРОВЕРКА РАБОТЫ СИСТЕМЫ

Для проверки работы системы было записано при помощи физической модели конструкций «с шумами» и программы Adobe Audition 3, два wave-файла с сигналами дефектов «стука» и «треска» по 10 секунд каждый, информацию о которых имеется в банке дефектов. Частота дискретизации wave-файлов равна 44100Гц, что допустимо для программы. Проведем диагностику этих сигналов.

Загрузим wave-файл «Стуки.wav» содержащий звуковые сигналы стуков. Посмотрим график сигнала с помощью кнопки «График», чтобы убедиться о правильности импортирования wave-файла. График представлен на рис. 3

 

Рис. 3 График сигнала  из wave-файла «Стуки.wav»

 

Как видно из графика, сигнал загружен правильно  и содержит в себе 11 ярко выраженных всплеска амплитуды. Загрузим банк дефектов  Data.xlsx с помощью меню и проведем диагностику данного сигнала, нажав кнопку «Диагностика». Результат диагностики представлен на рис. 4.

 

Рис. 4 Результат  диагностики сигнала «Стук.wav»

 

Время процедуры  диагностирования 2 мин. 27 сек. Как видно из рисунка 4, программа выдала гораздо больше сообщений о наличие присутствия признаков дефекта «стук», чем самих звуков данного дефекта в сигнале. Это обусловлено тем, что программа сверяет последовательно отрезки сигнала по 0,2сек., что не предусматривает того, что сигнал дефекта может быть длиннее этого отрезка и находиться одновременно в нескольких таких рядом стоящих отрезках. В последующей модернизации сигнала необходимо учесть данный факт. Однако проанализировав данный сигнал на слух с помощью Adobe Auditio 3 убеждаемся, что все сообщения, выданные программой DiagSound соответствуют действительности.

Загрузим wave-файл «tre.wav» с набором сигналов дефекта «Треск». Построим его график. График представлен на рис. 5

Рис.5 График сигнала из wave-файла «tre.wav»

 

Как видно из графике, сигнал не имеет больших  пиков, как в случае с «Стук.wav», что говорит о низкой амплитуде сигналов дефектов. Проведем диагностик данного сигнала. Результат диагностики сигнала представлен на рис. 6. Время диагностики 3мин. 17сек.

 

Рис. 6 Результат диагностики  сигнала «tre.wav»

 

Как видно из рисунка 6 сообщений так же как  и в первом случае выдано больше, чем ожидалось и процент схожести ниже. При анализе на слух данного сигнала с помощью Adobe Audition обнаруживаем наличие в сигнале постороннего шума (шорканье ногтей по корпусу гитары, скрип струн, щелчки, возникшие из-за несовершенства АЦП и др.). Данная программа не учитывает погрешности вносимые внешним воздействием, что говорит о несовершенстве данной программы. При последующей модернизации программы необходимо будет учесть данный факт и добавить модуль по выявлению и очищению сигнала от постороннего шума.

Так же необходимо отметить ресурсоемкость данной программы. На рис. 7 видно, что программа из-за больших объемов данных занимает около 144мб оперативной памяти, что не приемлемо для практического применения. Так же видно, что большая часть времени тратиться на работу с банком дефектов реализованного в Excel-таблице. При модернизации программы необходимо создать подходящую базу данных, обеспечивающую более быструю работу с данными.

 

Рис. 7 Ресурсоемкость программы DiagSound

 

 

 

ЗАКЛЮЧЕНИЕ

В данной курсовой работе была создана автоматизированная система диагностирование дефектов конструкций электронных устройств и был составлен банк дефектов, реализованный в Excel-таблице, содержащий информацию о двух видов дефектов. Данная система демонстрирует работу модуля анализа Diag со звуковыми сигналами, работу с банком дефектов и алгоритм сравнения полученных данных с помощью модуля Diag. Таким образом были решены поставленные для курсовой работу задачи. Система не пригодна к практическому применению из-за следующих недостатков:

• Большие вычислительные затраты
• Отсутствия полноценной базы данных
• Отсутствие возможности записи сигнала
• Отсутствие модуля для выявления постороннего шума
• Низкая производительность
• Существенная погрешность

Для практического применения автоматизированную систему необходимо модернизировать и исключить перечисленные недостатки. А так же следует оптимизировать код для быстродействия достаточного для диагностики звукового сигнала в реальном времени.

 

СПИСОК ЛИТЕРАТУРЫ

1. В. В. Клюева «Неразрушающий контроль: Справочник в 7 томах» М.:Машиностроение, 2005.-829с
2. Справочник по технической акустике: Пер. с ием./ Под ред. Хекла и Х. А. Мюллера. – Л.: Судостроение, 1980.-440с., ил. 329.-ИСБН.
3. Борис Васильевич Павлов «Акустическая диагностика механизмов» Изд. «Машиностроение» Москва 1971г.
4. Увайсов Расул Исаевич «МЕТОД ДИАГНОСТИКИ ДЕФЕКТОВ БОРТОВЫХ РАДИОТЕХНИЧЕСКИХ УСТРОЙСТВ», Автореферат диссертации, Москва 2008г.
5. Кинтцель Т. «Программирование звука на ПК=A Programmer`s Guide to Sound»:Пер. С англ.-М.:ДМК Пресс, 2005.-432с., ил.
6. Ong Swee Eng “Use Matlab package to design software application for human voice synthesis”2011г.
7. Электронный ресурс: http://habrahabr.ru/post/144491/
8. Электронный ресурс: http://ru.wikipedia.org/wiki
9. Электронный ресурс: http://lnktd-opz.narod.ru/vd.html
10. Электронный ресурс: http://zetlab.ru/support/articles/programs.php
11. Электронный ресурс: http://matlab.exponenta.ru/
12. Электронный ресурс: http://www.steinberg.net
13. Электронный ресурс: http://gapeal.narod.ru/scince/01-scince.html
14. Электронный ресурс: http://habrahabr.ru/post/132487/
15. Электронный ресурс: http://matstats.ru/dover.html

 

ПРИЛОЖЕНИЕ

1. Модули MatLab
1. Листинг модуля Diag.m

 

function res=Diag(Name)

[y,fs]=wavread(Name);

speechIn=myVAD(y);

res=mfccf(13,speechIn,fs);

end

 

2. Листинг myVAD.m

 

function trimmedX = myVAD(x)

 

% Syntax:        trimmedSample = myVAD(samplex);

% This function accepts an audio sample 'samplex' as input and returns a

% trimmed down version with non-speech sections trimmed off. Also known as

% voice activity detection, it utilises the algorithm due to Rabiner &

% Sambur (1975)

 

Ini = 0.1;          % Initial silence duration in seconds

Ts = 0.01;          % Frame width in seconds

Tsh = 0.005;        % Frame shift in seconds

Fs = 16000;         % Sampling Frequency

counter1 = 0;   

counter2 = 0;

counter3 = 0;

counter4 = 0;

ZCRCountf = 0;      % Stores forward count of crossing rate > IZCT

ZCRCountb = 0;      % As above, for backward count

ZTh = 40;           % Zero crossing comparison rate for threshold

w_sam = fix(Ts*Fs);                   % No of Samples/window 

o_sam = fix(Tsh*Fs);                  % No of samples/overlap

lengthX = length(x);

segs = fix((lengthX-w_sam)/o_sam)+1;  % Number of segments in speech signal

sil = fix((Ini-Ts)/Tsh)+1;            % Number of segments in silent period

win = hamming(w_sam);

Limit = o_sam*(segs-1)+1;             % Start index of last segment

 

FrmIndex = 1:o_sam:Limit;             % Vector containing starting index for each segment

ZCR_Vector = zeros(1,segs);           % Vector to hold zero crossing rate for all segments

                                     

% Below code computes and returns zero crossing rates for all segments in

% speech sample  

for t = 1:segs

    ZCRCounter = 0; 

    nextIndex = (t-1)*o_sam+1;

    for r = nextIndex+1:(nextIndex+w_sam-1)

        if (x(r) >= 0) && (x(r-1) >= 0)

         

        elseif (x(r) >= 0) && (x(r-1) < 0)

         ZCRCounter = ZCRCounter + 1;

        elseif (x(r) < 0) && (x(r-1) < 0)

             elseif (x(r) < 0) && (x(r-1) >= 0)

         ZCRCounter = ZCRCounter + 1;

        end

    end

    ZCR_Vector(t) = ZCRCounter;

end

% Below code computes and returns frame energy for all segments in speech

% sample

Erg_Vector = zeros(1,segs);

for u = 1:segs

    nextIndex = (u-1)*o_sam+1;

    Energy = x(nextIndex:nextIndex+w_sam-1).*win;

    Erg_Vector(u) = sum(abs(Energy));

end

 

IMN = mean(Erg_Vector(1:sil));  % Mean silence energy (noise energy)

IMX = max(Erg_Vector);          % Maximum energy for entire utterance

I1 = 0.03 * (IMX-IMN) + IMN;    % I1 & I2 are Initial thresholds

I2 = 4 * IMN;

ITL = min(I1,I2);               % Lower energy threshold

ITU = 5 * ITL;                  % Upper energy threshold

IZC = mean(ZCR_Vector(1:sil));  % mean zero crossing rate for silence region

stdev = std(ZCR_Vector(1:sil)); % standard deviation of crossing rate for 

                                % silence region

IZCT = min(ZTh,IZC+2*stdev);    % Zero crossing rate threshold

indexi = zeros(1,lengthX);      % Four single-row vectors are created

indexj = indexi;                % in these lines to facilitate computation below

indexk = indexi;

indexl = indexi;

 

% Search forward for frame with energy greater than ITU

for i = 1:length(Erg_Vector)

    if (Erg_Vector(i) > ITU)

        counter1 = counter1 + 1;

        indexi(counter1) = i;

    end

end

ITUs = indexi(1);

% Search further forward for frame with energy greater than ITL

for j = ITUs:-1:1

    if (Erg_Vector(j) < ITL)

        counter2 = counter2 + 1;

         indexj(counter2) = j;

    end

end

start = indexj(1)+1;

 

Erg_Vectorf = fliplr(Erg_Vector);% Flips round the energy vector 

% Search forward for frame with energy greater than ITU

% This is equivalent to searching backward from last sample for energy > ITU

for k = 1:length(Erg_Vectorf)

    if (Erg_Vectorf(k) > ITU)

        counter3 = counter3 + 1;

        indexk(counter3) = k;

    end

end

ITUf = indexk(1);

% Search further forward for frame with energy greater than ITL

for l = ITUf:-1:1

    if (Erg_Vectorf(l) < ITL)

        counter4 = counter4 + 1;

        indexl(counter4) = l;

    end

end

 

finish = length(Erg_Vector)-indexl(1)+1;% Tentative finish index

    

% Search back from start index for crossing rates higher than IZCT

   

BackSearch = min(start,25);

for m = start:-1:start-BackSearch+1

    rate = ZCR_Vector(m);

    if rate > IZCT

        ZCRCountb = ZCRCountb + 1;

        realstart = m;

    end

end

if ZCRCountb > 3

    start = realstart;          % If IZCT is exceeded in more than 3 frames

                                % set start to last index where IZCT is exceeded

end

 

% Search forward from finish index for crossing rates higher than IZCT

FwdSearch = min(length(Erg_Vector)-finish,25);

for n = finish+1:finish+FwdSearch

    rate = ZCR_Vector(n);

     if rate > IZCT

        ZCRCountf = ZCRCountf + 1;

        realfinish = n;

    end

end

if ZCRCountf > 3

    finish = realfinish;        % If IZCT is exceeded in more than 3 frames

                                % set finish to last index where IZCT is exceeded

end

 

x_start = FrmIndex(start);      % actual sample index for frame 'start'

x_finish = FrmIndex(finish-1);  % actual sample index for frame 'finish'

trimmedX = x(x_start:x_finish); %T rim speech sample by start and finish indices

 

3. Листинг mfsscf.m

 

function FMatrix=mfccf(num,s,Fs)

 

% Syntax:        M=mfccf(num,s, Fs);

% Computes and returns the mfcc coefficients for a speech signal s

% where num is the required number of MFCC coefficients. It utilises the 

% function 'melbankm' from the toolbox 'Voicebox' by Mike Brooks 

 

n=512;              % Number of FFT points

Tf=0.025;           % Frame duration in seconds

N=floor(Fs*Tf);     % Number of samples per frame

fn=24;              % Number of mel filters

l=length(s);        % total number of samples in speech

Ts=0.01;            % Frame step in seconds

FrameStep=Fs*Ts;    % Frame step in samples

a=1;

b=[1, -0.97];       % a and b are high pass filter coefficients

noFrames=floor(l/FrameStep);    % Maximum no of frames in speech sample

FMatrix=zeros(noFrames-2, num); % Matrix to hold cepstral coefficients

lifter=1:num;                   % Lifter vector index

lifter=1+floor((num)/2)*(sin(lifter*pi/num));% raised sine lifter version

 

if mean(abs(s)) > 0.01

    s=s/max(s);                     % Normalises to compensate for mic vol differences

end

 

% Segment the signal into overlapping frames and compute MFCC coefficients

for i=1:noFrames-2

    frame=s((i-1)*FrameStep+1:(i-1)*FrameStep+N);  % Holds individual frames

    Ce1=sum(frame.^2);          % Frame energy

    Ce2=max(Ce1,2e-22);         % floors to 2 X 10 raised to power -22

    Ce=log(Ce2);

    framef=filter(b,a,frame);   % High pass pre-emphasis filter

    F=framef.*hamming(N);       % multiplies each frame with hamming window

    FFTo=fft(F,N);              % computes the fft

    melf=melbankm(fn,n,Fs);     % creates 24 filter, mel filter bank

    halfn=1+floor(n/2);    

    spectr1=log10(melf*abs(FFTo(1:halfn)).^2);% result is mel-scale filtered

    spectr=max(spectr1(:),1e-22);

    c=dct(spectr);              % obtains DCT, changes to cepstral domain

    c(1)=Ce;                    % replaces first coefficient

    coeffs=c(1:num);            % retains first num coefficients

     ncoeffs=coeffs.*lifter';    % Multiplies coefficients by lifter value

    FMatrix(i, :)=ncoeffs';     % assigns mfcc coeffs to succesive rows i

end

 

% Call the deltacoeff function to compute derivatives of MFCC

% coefficients; add all together to yield a matrix with 3*num columns 

d=(deltacoeff(FMatrix)).*0.6;   % Computes delta-mfcc

d1=(deltacoeff(d)).*0.4;        % as above for delta-delta-mfcc

FMatrix=[FMatrix,d,d1];         % concatenates all together

 

4. Листинг melbankm.m

 

function [x,mn,mx]=melbankm(p,n,fs,fl,fh,w)

% MELBANKM determine matrix for a mel-spaced filterbank [X,MN,MX]=(P,N,FS,FL,FH,W)

%

%  Inputs:   p   number of filters in filterbank

%       n   length of fft

%       fs  sample rate in Hz

%       fl  low end of the lowest filter as a fraction of fs (default = 0)

%       fh  high end of highest filter as a fraction of fs (default = 0.5)

%       w   any sensible combination of the following:

%             't'  triangular shaped filters in mel domain (default)

%             'n'  hanning shaped filters in mel domain

%             'm'  hamming shaped filters in mel domain

%

%             'z'  highest and lowest filters taper down to zero (default)

%             'y'  lowest filter remains at 1 down to 0 frequency and

%              highest filter remains at 1 up to nyquist freqency

%

%              If 'ty' or 'ny' is specified, the total power in the fft is preserved.

%

%  Outputs:  x     a sparse matrix containing the filterbank amplitudes

%             If x is the only output argument then size(x)=[p,1+floor(n/2)]

%             otherwise size(x)=[p,mx-mn+1]

%       mn   the lowest fft bin with a non-zero coefficient

%       mx   the highest fft bin with a non-zero coefficient

%

%  Usage:    f=fft(s);           f=fft(s);

%       x=melbankm(p,n,fs);     [x,na,nb]=melbankm(p,n,fs);

%       n2=1+floor(n/2);        z=log(x*(f(na:nb)).*conj(f(na:nb)));

%       z=log(x*abs(f(1:n2)).^2);

%       c=dct(z); c(1)=[];

%

%  To plot filterbanks e.g.  plot(melbankm(20,256,8000)')

%

%      %      Version: $Id: melbankm.m,v 1.3 2005/02/21 15:22:13 dmb Exp $

%

%   VOICEBOX is a MATLAB toolbox for speech processing.

%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html

%

if nargin < 6

  w='tz';

  if nargin < 5

    fh=0.5;

    if nargin < 4

      fl=0;

    end

  end

end

f0=700/fs;

fn2=floor(n/2);

lr=log((f0+fh)/(f0+fl))/(p+1);

%  convert to fft bin numbers with 0 for DC term

bl=n*((f0+fl)*exp([0 1 p p+1]*lr)-f0);

b2=ceil(bl(2));

b3=floor(bl(3));

if any(w=='y')

  pf=log((f0+(b2:b3)/n)/(f0+fl))/lr;

  fp=floor(pf);

  r=[ones(1,b2) fp fp+1 p*ones(1,fn2-b3)];

  c=[1:b3+1 b2+1:fn2+1];

  v=2*[0.5 ones(1,b2-1) 1-pf+fp pf-fp ones(1,fn2-b3-1) 0.5];

  mn=1;

  mx=fn2+1;

else

  b1=floor(bl(1))+1;

  b4=min(fn2,ceil(bl(4)))-1;

   pf=log((f0+(b1:b4)/n)/(f0+fl))/lr;

  fp=floor(pf);

  pm=pf-fp;

  k2=b2-b1+1;

  k3=b3-b1+1;

  k4=b4-b1+1;

  r=[fp(k2:k4) 1+fp(1:k3)];

  c=[k2:k4 1:k3];

  v=2*[1-pm(k2:k4) pm(1:k3)];

  mn=b1+1;

  mx=b4+1;

end

if any(w=='n')

  v=1-cos(v*pi/2);

elseif any(w=='m')

  v=1-0.92/1.08*cos(v*pi/2);

end