Мне как-то для подобной проблемы (АОH+dialtone detection+voice detection на
зюхе) понадобилось то же самое.
Так я тоже сначала начал Фурье делать. Для АОH прошло классно (там всего 6
частот, причем строго фиксированных), для диалтона хуже (частота плавает в
достаточно широких пределах), а вот с войсом проблемы.
Так мне недавно рассказали алгоритм, от которого я просто приторчал.
Суть: резонансный фильтр. Т.е. просто "маятник" с собственной резонансной
частотой, на который действуем нашим сигналом и смотрим получившуюся амплитуду.
Он выигрывает в том, что покрывает сразу полосу частот, причем ее ширина
конфигурабельна. Для диалтона просто незаменимо, для голоса их, конечно,
несколько надо делать.
Реализация совершенно фанатна.
x[i+1]=x[i]-x[i-1]. Это синусоида с шестью точками на период (легко проверяется
через тригонометрические тождества). К ней добавляется исходный сигнал и
затухание, и через некоторое время смотрим получившуюся амплитуду.
Это делалось очень давно на какой-то машине, на которой не было плавающей точки
и не было умножения. Так это все делается сдвигами:
#define k 4
#define T 200
int x0=0, x1=0, x2, t;
for (t=0;t<T;t++)
{
x2=x1-x0+F();
x2-=(x2>>k);
x0=x1;
x1=x2;
}
где F() - исходный сигнал,
T - время опроса в элементарных единицах,
1/2**k - коэффициент затухания, он же (точнее, 2**k) - добротность фильтра (т.е.
отношение амплитуды в резонансе к амплитуде при шуме). Если, конечно, я ничего
не перепутал, облом проверять.
Мне это рассказал Э.М.Куссуль, (c) его, судя по всему.
По сравнению с "честным" фильтром получается единственная неприятность: всплески
амплитуды на частотах F*13 и F*27. Для большинства применений это несущественно.
Еще момент: нужно выбирать F() с нужной частотой (в шесть раз большей частоты
фильтра), из-за чего теряется полезная информация. Hа самом деле, это влияет
только на то, что помехи будут играть большую роль. От этого можно избавиться,
беря за F() среднее из близких по времени результатов опроса, чтобы не просто
выкидывать "ненужные". Если же частота опроса слишком низка (например, 9600 при
том, что надо определять 3000), то промежуточные значения, которых не хватает,
можно интерполировать.
Lucky carrier,
Гуля
Algorithm
Begin WINDOW
allocate array space in VIFF structure for window data
(1-D array of WINDOW_NPTS points)
write WINDOW_TYPE and WINDOW_NPTS to VIFF header
(from values in WTYPE and NPTS)
write PLOT1_XTITLE, PLOT1_ATITLE to VIFF header
set XMIN=1.0, XMAX=NPTS
write PLOT1_XMIN, PLOT1_XMAX to VIFF header
Case switch on window type (WINDOW_TYPE)
user spec: set WINDOW_IFILE = I_FILE
read window values from WINDOW_IFILE
transfer values into IMG1 array area
write WINDOW_IFILE to VIFF header
Rectangle: compute rectangle window values in IMG1
write WINDOW_SIDX, WINDOW_EIDX to VIFF header
(from values in SIDX, EIDX)
calculate the window as:
for i=1 to SIDX-1
IMG1(i) = (0.0, 0,0)
end for
for i=SIDX to EIDX
IMG1(i) = (1.0, 0.0)
end for
for i=EIDX+1 to NPTS
IMG1(i) = (0.0,0.0)
end for
Triangle: compute triangle window in IMG1
calculate the window as:
slope=2./FLOAT(NPTS-1)
for i=1 to (NPTS+1)/2
j=NPTS+1-i
IMG1(i) = ((i-1)*slope,0.0)
IMG1(j) = IMG1(i)
end for
Hamming: compute Hamming window in IMG1
calculate the window as:
for i=1 to NPTS
arg = (2*PI)*FLOAT(I-1)/FLOAT(NPTS-1)
W = .54-.46*COS(arg)
IMG1(i) = (W, 0.0)
end for
Hanning: compute Hanning window in IMG1
calculate the window as:
for i=1 to NPTS
arg = (2*PI)*FLOAT(I-1)/FLOAT(NPTS-1)
W = .5(1-COS(arg))
IMG1(i) = (W, 0.0)
end for
Kaiser: set WINDOW_KMAX TO 11 (sets precision)
compute Kaiser window in img1
write ALPHA as WINDOW_ALPHA to VIFF header
write WINDOW_KMAX to VIFF header
calculate the window as:
for i=1 to NPTS
arg = 1.0 -
(FLOAT(NPTS/2+1)/FLOAT(NPTS/2))**2
W = IO(PI*ALPHA*SQRT(arg))/IO(PI*ALPHA)
IMG1 = (W, 0.0)
end for
where IO is the zero-order modified Bessel
function of the first kind
Taylor: compute Taylor window in IMG1
write WINDOW_SLB, WINDOW_NBAR to VIFF header
(from values in SLB and NBAR)
calculate the window as:
for i=1 to NPTS
arg = 2*PI*FLOAT(I-(NPTS/2 + 1))/FLOAT(NPTS)
W(i) = F(1)
for n=2,NBAR
W(i) = W(i) + F(n)*cos((n-1)*arg)
end for
end for
where F(n) are the Taylor coefficients
calculated as:
NBAR
[(NBAR-1)!]**2 * PI (1 - j**2/Zm**2)
m = 1
F(n) = ---------------------------------
(NBAR-1+n)! * (NBAR-1-n)!
and
Zm = sigma*sqrt{(A**2 + (m -.5)**2}
A = (1/PI) * inv_cosh(R)
R = 10**(PSLL_DB/20)
Dolph Chebyshev:
compute Dolph Chebyshev window in IMG1
write PSLL to VIFF header
calculate the window as:
N = NPTS
for i=1 to NPTS
W(i) =
L
2/N *[S + 2* SUM{T_N-1(arg1(n,N))*cos(arg2)}]
n = 1
end for
where S is the mainlobe to sidelobe ratio given as
S = 10 ** PSLL_DB/20
T_m is the Chebyshev Polynomial given as
T_m(z) = cos(m*inv_cos(z)) for |z| le 1
= cosh(m*inv_cosh(z)) for |z| gt 1
arg1(n,N) = Zo*cos(n*PI/N)
arg2 = 2*PI*n/N * (k - (N+1)/2)
and
Zo = cosh[inv_cosh(S)/(N-1)]
L = (N-1)/2 for N odd
= N/2 - 1 for N even
End Case switch on window type
write appropriate string as PLOT1_TITLE to VIFF header
write IMG1 to output file O_FILE
End WINDOW