Тестирование алгоритма
Тестирование производилось с 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 комментариев