Идентификация пользователя по голосу [3]. Тестирование и исходники

Программирование
Продолжение темы «Идентификация пользователя по голосу»
Идентификация пользователя по голосу [1]. Фильтрация и разложение спектра звука
Идентификация пользователя по голосу [2]. Анализ спектра.

Тестирование алгоритма
Тестирование производилось с 8 пользователями. Каждый голос сначала сравнивался с эталонным, то есть голосом разработчика, а потом между собой, для того что бы выяснить как поведет себя система на однотипных голосах.
При тестировании использовались 6 мужских голосов и 2 женских. Схожесть голосов определяется в процентах, поэтому требовалось выяснить максимально возможный порог совпадения. Эталонный голос использовался мужской, поэтому для тестирования использовалось большое количество именно мужских голосов. При тестирование произносилась одна и та же кодовая фраза, которую я за много прошедших лет уже и не помню…

Графики спектральных характеристик визуально различаются достаточно сильно, но положение пиков у них абсолютно одинаковое. Именно поэтому на одинаковых фразах даже пользователь с похожим голосом не сможет добиться такой схожести. На его характеристике положение этих пиков совпадать не будет. Так же на спектрограммах видно что произносились фразы по разному, первый образец был самых отчетливый, второй был сказан с некоторым отдалении от микрофона, третий произнесен шепотом. Это должно было сильно усложнить задачу. Но как видно из графиков их спектральные характеристики оказались схожими.

Тестирование проводилось на очень слабой звуковой карте интегрированной в материнскую плату. Карточка с высоким уровнем шума и игнорированием высоких и низких частот. А также со слабым микрофоном, не обеспечивающим необходимый уровень записи. С хорошей звуковой подсистемой, можно добиться значительно лучших результатов.

Диалог добавления пользователя в систему


Анализатор голоса


Тестер работы алгоритма анализа спектра


Листинг модуля преобразования Фурье. Комплировал в dll, для вызова в дальнейшем из Delphi.
//---------------------------------------------------------------------------
//  VoiceSec module. FFT transform. 
//  10.2002
//
//  author email viman@pisem.net
//
//---------------------------------------------------------------------------

#include <vcl.h>
#include "fft.h"
#pragma hdrstop
#define WINDOW_RECTANGLE	0
#define WINDOW_HAMMING	1
#define WINDOW_BLACKMAN	2
#define WINDOW_BARTLETT	3
#define WINDOW_TRIANGLE	3
#define WINDOW_HANNING	4

#define MATH_PI	3.141592653589793


inline int	Round(double x) { return (int)(x+0.5); }
inline long double	sqr(long double x) { return x*x; }

int FFTSize;
#ifdef __cplusplus
extern "C" {
#endif
__declspec( dllexport )
void FFTC(Cmplx *X,int FFTSize)
{
	int	N, M, i, j, L, LE, LE2, ip, k, s;
	Cmplx t,z;
	RealData	UR, UI, SR, SI, TR, TI;

	N = FFTSize;
	M = Round(log(N)/log(2));
	// Bit-reverse
	i = 0;
	for (s=0;s<N-1;s++) {
		if (s<i) {
			t = *(X+i); *(X+i) = *(X+s); *(X+s) = t;
		}
		k = N >> 1;
		while (i&k) k >>= 1;
		i += k;
		k <<= 1;
		while (k<N) {
			i -= k;
			k <<= 1;
		}
	}
	// First pass
	for (i=0;i<N;i+=2) {
		t = *(X+i);
		(X+i)->Re = t.Re + (X+i+1)->Re;
		(X+i)->Im = t.Im + (X+i+1)->Im;
		(X+i+1)->Re = t.Re - (X+i+1)->Re;
		(X+i+1)->Im = t.Im - (X+i+1)->Im;
	}
	// Second pass
	for (i=0;i<N;i+=4) {
		t = *(X+i);
		(X+i)->Re = t.Re + (X+i+2)->Re;
		(X+i)->Im = t.Im + (X+i+2)->Im;
		(X+i+2)->Re = t.Re - (X+i+2)->Re;
		(X+i+2)->Im = t.Im - (X+i+2)->Im;
		t = *(X+i+1);
		z = *(X+i+3);
		(X+i+1)->Re = t.Re + z.Im;
		(X+i+1)->Im = t.Im - z.Re;
		(X+i+3)->Re = t.Re - z.Im;
		(X+i+3)->Im = t.Im + z.Re;
	}
	// Last passes
	for (L=3;L<=M;L++) {
		LE = 1 << L;
		LE2 = LE >> 1;
		UR = 1; UI = 0;
		SR = cos(MATH_PI/LE2);
		SI =	-sin(MATH_PI/LE2);
		for (j=0;j<LE2;j++) {
			for (i=j;i<N;i+=LE) {
				ip = i + LE2;
				TR = (X+ip)->Re*UR - (X+ip)->Im*UI;
				TI = (X+ip)->Re*UI + (X+ip)->Im*UR;
				(X+ip)->Re = (X+i)->Re - TR;
				(X+ip)->Im = (X+i)->Im - TI;
				(X+i)->Re = (X+i)->Re + TR;
				(X+i)->Im = (X+i)->Im + TI;
			}
			TR = UR;
			UR = TR*SR - UI*SI;
			UI = TR*SI + UI*SR;
		}
	}
}
}
#ifdef __cplusplus

#endif

#ifdef __cplusplus
extern "C" {
#endif
__declspec( dllexport )void
FFTR(Cmplx *X) {
	int	N, ND2, ND4;
	int	i, im, ip2, ipm, ip;
	RealData	UR, UI, SR, SI, TR, TI;
    FFTSize=256;
	// Separate even and odd points
	N = FFTSize;
	ND2 = N>>1;
	ND4 = ND2>>1;
	for (i=0;i<ND2;i++) {
		(X+i)->Re = (X+(i<<1))->Re;
		(X+i)->Im = (X+(i<<1)+1)->Re;
	}
	// Calculate N/2 point FFT
	FFTSize = ND2;
	FFTC(X,FFTSize);
	FFTSize = N;
	// Even/odd frequency domain decomposition
	for (i=1;i<ND4;i++) {
		im = ND2 - i;
		ip2 = i + ND2;
		ipm = im + ND2;
		(X+ipm)->Re = (X+ip2)->Re = ((X+i)->Im + (X+im)->Im)*0.5;
		(X+ip2)->Im = ((X+i)->Re - (X+im)->Re)*(-0.5);
		(X+ipm)->Im =  - (X+ip2)->Im;
		(X+im)->Re = (X+i)->Re = ((X+i)->Re + (X+im)->Re)*0.5;
		(X+i)->Im = ((X+i)->Im - (X+im)->Im)*0.5;
		(X+im)->Im = - (X+i)->Im;
	}
	(X+N*3/4)->Re = (X+ND4)->Im;
	(X+ND2)->Re = X->Im;
	(X+ND2+ND4)->Im = (X+ND2)->Im = (X+ND4)->Im = X->Im = 0;
	// Complete the last FFT stage
	// First step: calculate X[0] and X[N/2]
	TR = (X+ND2)->Re;
	TI = (X+ND2)->Im;
	(X+ND2)->Re = X->Re - TR;
	(X+ND2)->Im = X->Im - TI;
	X->Re = X->Re + TR;
	X->Im = X->Im + TI;
	// Other steps
	UR = SR = cos(MATH_PI/ND2);
	UI = SI = -sin(MATH_PI/ND2);
	ip = ND2+1;
	for (i=1;i<ND2;i++,ip++) {
		TR = (X+ip)->Re*UR - (X+ip)->Im*UI;
		TI = (X+ip)->Re*UI + (X+ip)->Im*UR;
		(X+ip)->Re = (X+i)->Re - TR;
		(X+ip)->Im = (X+i)->Im - TI;
		(X+i)->Re = (X+i)->Re + TR;
		(X+i)->Im = (X+i)->Im + TI;
		(X+i)->Re = (X+i)->Re + (X+ip)->Re*UR - (X+ip)->Im*UI;
		(X+i)->Im = (X+i)->Im + (X+ip)->Re*UI + (X+ip)->Im*UR;
		TR = UR;
		UR = TR*SR - UI*SI;
		UI = TR*SI + UI*SR;
	}
}}
#ifdef __cplusplus

#endif


//---------------------------------------------------------------------------
int WINAPI DllEntryPoint(HINSTANCE hinst, unsigned long reason, void*)
{
    return 1;
}
//---------------------------------------------------------------------------



Листинг модуля Spectr. В нем алгоритм анализа спектра, который в дальнейшем и вызывается из других модулей.
unit uSpectr;


interface
uses Sysutils,graphics,extctrls,windows;

procedure GetSpectr(filename:string;var spectr:array of integer);
function calc_fxy(energ11,energ21:array of integer):real;

implementation

type arr=array[0..10000]of byte;
     sampl=array[0..1024]of real;
     samplr=array[0..1024]of real;
     comp=record re,im:real;end;

procedure FFTR(buf:array of comp);stdcall; external 'fftc.dll';


const k=0.5;

type pikrec=record
  pos:integer;
  lpc:integer;
  l:integer;
end;

var  left:array[0..200000]of smallint;
     buf:array[0..200000]of byte;

procedure savevsd(s:string;energ:array of integer);
var f:file;
    x:integer;
begin
     assignfile(f,s);
     rewrite(f,1);
     blockwrite(f,energ,sizeof(energ),x);
     closefile(f);
end;

function calc_fxy(energ11,energ21:array of integer):real;
var mx,my,s1,s2,s3:real;
    i:integer;
begin
    mx:=0;my:=0;
    for i:=0 to 127 do
     begin
       mx:=mx+energ11[i];
       my:=my+energ21[i];
     end;
    mx:=mx/128;
    my:=my/128;
    s1:=0;s2:=0;s3:=0;
    for i:=0 to 127 do s1:=s1+((energ11[i]-mx)*(energ21[i]-my));
    for i:=0 to 127 do s2:=s2+(sqr(energ11[i]-mx));
    for i:=0 to 127 do s3:=s3+(sqr(energ21[i]-my));
    s2:=sqrt(s2);
    s3:=sqrt(s3)*s2;
    if s3<>0 then result:=abs(s1/s3) else result:=0;
end;

function I0(X:real):real;
var Y, T, E, DE, SDE:real;
     i:integer;
begin
	Y := X/2;
	T := 0.0000000001;
        de:= 1;
	E :=DE;
	for i:=1 to 50 do
         begin
	  DE := de*Y / i;
          SDE := DE * DE;
          E := e+SDE;
          if (E*T-SDE>0) then break;
	end;
	result:= E;
end;

procedure KaiserWindow(source:sampl; length,b:integer;var result:sampl);
var i, k:integer;
   iI0b:real;
   h:real;
begin
  iI0b:=1.0 / I0(b);
  k := -(Length shr 2);
  for i:=1 to Length shr 2 do
   begin
    inc(k);
    h := I0(b*sqrt(1-sqr(2.0*k/(Length-1)))) * iI0b;
    source[i]:=source[i] *h;
    source[Length-1-i] := source[Length-1-i]*h;
   end;
   for i:=1 to length do result[i]:=source[i];
end;

procedure Liftering(source:sampl; length,b:integer;var result:sampl);
var i:integer;
   w,xout:real;
begin
  for i:=1 to Length do
   begin
      xout:=source[i]-0.9*w;
      w:=0.54-0.46*cos(2*pi*(i-6)/179);
      source[i]:=xout*w;
   end;
   for i:=1 to length do result[i]:=source[i];
end;

procedure Fft(source:sampl; length,step:integer;var result:sampl);forward;
procedure Fft(source:sampl; length,step:integer;var result:sampl);
var  f,t:integer;
begin
  for f:=1 to Length do
    begin
     Result[f]:=0.0;
     for t:=1 to Length do
        Result[f]:=Result[f]+ Source[t]*(-cos(PI*f*t/Length/2));
    end;
end;

procedure GetSpectr(filename:string;var spectr:array of integer);

var step                        :integer;
   sImage                       :array[1..500,1..128]of smallint;
   x,y,n                     :integer;
   energ                        :array[0..128]of integer;
   bufc                         :array [0..512]of comp;
   f                            :file;
   chastota                     :sampl;
   l                            :integer;
   begin_crop,end_crop          :integer;
   max,min                      :integer;
   xout,dz                      :real;
   I                            :Integer;
   xx                           :integer;

begin
  assignfile(f,filename);
  reset(f,1);
  for n:=1 to 200000 do buf[n]:=0;
  for n:=1 to 200000 do left[n]:=0;
  seek(f,58);
  BlockRead(F, buf, filesize(f), x);
  end_crop:=x-100;
  closefile(f);
  for x:=1 to 128 do energ[x]:=0;
  x:=1;
  begin_crop:=200;
  max:=0;
  min:=255;

  for l:=begin_crop to end_crop do if buf[l]<>0 then break;
  begin_crop:=l;
  for l:=begin_crop to end_crop do if buf[l]<min then min:=buf[l];
  for l:=begin_crop to end_crop do
       left[l]:=round((buf[l]+buf[l-1]+buf[l+1])/3)-128;

  for l:=begin_crop to end_crop do
    begin
       xout:=LEFT[l]-0.90*dz;
       dz:=left[l];
       LEFT[l]:=round(xout*(0.54-0.46*cos((l-6)*6.2832/179)));
    end;

  repeat
   inc(begin_crop);
  until abs(left[begin_crop]-left[begin_crop+1])>4;

  repeat
   dec(end_crop);
  until abs(left[end_crop]-left[end_crop-1])>4;
  begin_crop:=begin_crop+512;
  end_crop:=end_crop-512;


  for l:=begin_crop to end_crop do
   begin
    if left[l]>max then max:=left[l];
   end;

  step:=round((end_crop-begin_crop)/500);
  for x:=1 to 500 do for y:=1 to 128 do
   begin
     sImage[x,y]:=0;
   end;
 step:=64;

  for y:=1 to 128 do  energ[y]:=0;

  //-=-=-=-=-=  Analize

  for x:=1 to 500 do
   begin

  if x*step+begin_crop< end_crop then begin

        for i:=0 to 511 do
         begin
         bufc[I].Re:= left[i+x*step+begin_crop];
         bufc[I].Im:=0;
         end;

  FFTR(bufc);//----

  for i:=128 to 256 do chastota[i-128]:=(sqrt(sqr(bufc[I].Im)+sqr(bufc[I].re)));

  Liftering(chastota,128,5,chastota);
  KaiserWindow(chastota,128,10,chastota);

//////////////////////
     for y:=1 to 125 do
      begin
      l:=abs(round(chastota[y]));
      energ[y]:=energ[y]+l;
      if l>255 then l:=255;
      if l<0 then l:=0;
      end;
   end;
end;

//-=-=-=-=-=-=- end analize

     energ[126]:=0;

    for x:=1 to 125 do
    begin
      if energ[x]>y then y:=energ[x];
    end;
    y:=round(y/128);
     xx:=1;
    if y>0 then
      begin
        for xx:=1 to 128 do ;//spectr[xx]:=round(energ[xx]/y);
        for xx:=0 to 127 do spectr[xx]:=round(energ[xx+1]/y);
        for xx:=1 to 128 do energ[xx]:=round(energ[xx]/y);
      end;

     SaveVSD(ChangeFileExt(filename,'.vsd'),energ);
end;

end.


Скачать текст диплома, чертежи, программу с исходниками

12 комментариев