Декодирование DTFM

← →
Вадим © (2004-05-22 17:34) [0]

Здравствуйте. Задача стоит такая. На линейный вход приходит последовательность цифр в виде DTFM. Необходимо в реальном времени преобразовывать это снова в цифры. DTFM представляет собой смесь двух частот (одной нижней и одной верхней) из следующего ряда: 697, 770, 852, 941 (нижние), 1209, 1336, 1477 и 1633 Гц (верхние). Что позволяет закодировать 16 символов: 1, 2, 3, A, 4, 5, 6, B, 7, 8, 9, C, *, 0, #, D. Таким образом, проблема сводится к тому, чтобы обнаружить наличие или отсутствие данных частот во входном сигнале, но как это сделать не представляю… В обсуждении темы fft.dll (http://delphimaster.net/view/8-1082816087/) видел процедуру, делающую БПФ, и по этому (http://www.delphikingdom.com/asp/answer.asp?IDAnswer=5489) адресу было что-то подобное, но как этим воспользоваться до меня не доходит. Либо есть более простой выход?

← →
Jel © (2004-05-22 17:58) [1]

Посмотри готовые компоненты тут
http://www.torry.net/quicksearchd.php?SID=f8370ce5e349ab4b6be3a0a0bce1fd4b&String=DTMF&Title=No

← →
Вадим © (2004-05-22 21:07) [2]

To Jel. А ты сам-то видел, то, что советуешь?! Там либо генераторы DTFM, либо компоненты, которые "весят" мегабайты и стоят совсем нехило! Нужен дельный совет, расчитанный на мою пока не высокую квалификацию.

← →
Jel © (2004-05-22 22:16) [3]

> А ты сам-то видел, то, что советуешь?!

Вообще-то видел. Зря шумишь. Ты не говорил что тебе нужны бесплатные решения не весящие мегабайты. Ты хотел простой выход - я тебе его подсказал. Подходит - отлично, нет - ищи другой.
:)

← →
Вадим © (2004-05-22 23:16) [4]

> Ты хотел простой выход - я тебе его подсказал

Спасибо, но давай не устраивать полемику. Я же в вопросе даже дал ссылки на сравнительно несложные алгоритмы для преобразования БПФ. Мне нужна хотя бы подсказка как их применить для решения моей проблемы. Либо что-нибудь более простое. Ещё раз повторю, буду признатален за конкретную помощь или подсказу

← →
Вадим © (2004-05-23 21:54) [5]

Я, конечно, извиняюсь, но неужели не найдется человека, знающего как мне помочь?

← →
Evgeny V © (2004-05-24 06:21) [6]

В данном случае лучше дискретное преобразование фурье во скользящему окну(это точнее позволит определить момент смены сигнала), приведу пример для нахождения одной заданной частоты по N значениям, например частота дискретизации 8000, берем выборки из 64 значений( это цифра должна быть степенью 2 и достаточно большая, что бы обеспечить точность измерения частоты) и по ним делаем преобразование, шаг между частотами по спектру в этом случае 8000/64 125 герц (фактически это ширина полосы определяемой частоты, спектральной палки, точность), делаем массив выборок для мнимой и действительно части из 64 элеменотов, тогда делаем следующее (пример для одной частоты, посчитаешь сам )
EnX:array [0..63] double - сумма коэффициентов реальной части для заданной частоты
EnУ:[0..63] double - сумма коэффициентов мнимой части для заданной частоты
An:double - Энергия заданной частоты
buf:array [0..63] - твои выборки сигнала с саунда.
porog;double - величина, с которой ты сравниваешь на наличие частоты
j:integer - циклический счетчик, местонахождение фазы измерения, свяжем его с частотой дискретизации
Freq - заданная частота
Вначале все переменные, массивы раны нулю, кроме porog и Freq, это уже задается;

Предварительно заполни массив buf - очередной порцией выборок саунда - 64 выборки, можно и больше взять число, 128, 256 и т.д, влияет на точность определения и временной интевал, за который определяем
for i:=0 to 63 do
begin
j:=j mod 8000;
An:=An-EnX[i]*EnX[i]-EnY[i]*EnY[i];// - вычитаем из суммы коэффициентов самое старое значение, на его место запишем следом новое
EnX[i]:=buf[i]*Cos(2*pi*Freq*j/8000.0);// вычислим реальную часть
EnY[i]:=buf[i]*Sin(2*pi*Freq*j/8000.0);// вычислим мнимую часть
An:=An+EnX[i]*EnX[i]+EnY[i]*EnY[i];// значение энергии в заданной точке по результату 64 выборок
Inc(j);
if An >= porog then частота есть в сигнале
end;

Получим новые 64 выборки и повторим вычисления и так до .....

Вот весьма упрощенный алгоритм определения частоты в сигнале. Если ты не поленишься и построишь реальные графики(Tchart например весьма для этого удобен) сигнала и этих вычислений (можешь провести их для каждой частоты DTMF), то наглядно увидишь что и как, это лучше чем рассказывать, поиграйся с размером массива, можно и не 64, это дело твое, насколько помню, в DTMF частоты могут лежать довольно близко в строках или стобцах, и тебе лучше будет брать размер 128 точек или больше.

← →
Jeer © (2004-05-24 15:38) [7]

Алгоритм Герцеля (Goerzel)
http://factaee.elfak.ni.ac.yu/fu2k33/9mp.html

← →
Вадим © (2004-05-27 22:30) [8]

To Evgeny V. Извиняюсь, что долго думал, но так уж получилось. Вот что я делал, чтобы "поиграться":
{$N+} Uses CRT; var i,j: integer; EnX: array [0..63] of double; EnY: array [0..63] of double; An, Freq: double; buf: array [0..63] of double; begin ClrScr; j:=0; for i:=0 to 63 do begin Enx[i]:=0; Eny[i]:=0; buf [i]:= sin (2*pi*i*8/64);//- в ручную заполняю буфер из расчета - входящяя 1000 и дескридитация 8000 Герц write (i," ",buf[i]:5:3," ");//- посмотреть, чем заполнил end; begin writeln; An:=0; Freq:=1000; for i:=0 to 63 do begin j:=j mod 8000; An:=An-EnX[i]*EnX[i]-EnY[i]*EnY[i]; EnX[i]:=buf[i]*Cos(2*pi*Freq*j/8000.0); EnY[i]:=buf[i]*Sin(2*pi*Freq*j/8000.0); An:=An+EnX[i]*EnX[i]+EnY[i]*EnY[i]; write (i," ",EnX[i]:5:3," "); Inc(j); end; writeln; writeln (An:5:3); end; end.
Делалось это в Борланде. Так вот, какое бы значение Freq я не задавал, An остается одно и то же! Где здесь ошибка? В алгоритме ДПФ или я сделал что-то не так? И зачем былаа введена переменная j, ведь она в цикле принимает те же значения, что и i? А вообще спасибо, если это сработает, то будет то, что надо!

← →
KilkennyCat © (2004-05-28 05:10) [9]

Наиболее простой способ обойдется в 50 рублей, час работы паяльником.... микруха декодера DTMF -> LPT

← →
Evgeny V © (2004-05-28 05:46) [10]

Ошибки нет, дело в том что энергия не появится сразу, и в начале у тебя все инициализировано 0 и с первого буфера An будет расти медленно, ты сделай буфер на секунду, и в нем только на 100 миллисикунд заполни значения частотой, например в середине, тогда будет видно изменения. Про j когда будешь бежать по буферу более чем 64 выборки, 12000 например, то j не будет равно i

for i:=0 to N-1 do
begin
j:=j mod 8000;
An:=An-EnX[i mod 64]*EnX[i mod 64]-EnY[i mod 64]*EnY[i mod 64];
EnX[i mod 64]:=buf[i]*Cos(2*pi*Freq*j/8000.0);
EnY[i mod 64]:=buf[i]*Sin(2*pi*Freq*j/8000.0);
An:=An+EnX[i mod 64]*EnX[i mod 64]+EnY[i mod 64]*EnY[i mod 64];
write (i," ",EnX[i mod 64]:5:3," "); - это кстати не очень наглядный коэффициент, энергия должна сильно меняться
write (i," ",EnX[i mod 64]:5:3," "); - аналогично
Inc(j);
end;
writeln;
writeln (An:5:3);
end;
N -cделай побольше,ну пусть 4000, и тысячу значений в нем например заполни синусом нужной частоты. Посмотри результат тогда, если не получится, пришлю тебе прогу и исходник, где я делаю похожие вещи. Удачи, пиши, можешь на мыло в моей анкете.

← →
Вадим © (2004-05-28 12:33) [11]

To KilkennyCat: Не такой уж и проостой. Ведь внешние устройства нужно програмно обслуживать. Я с COM-портом кое-как разобрался, да и то не до конца, а тут ещё LPT! Так зачем городить огород, если можно всё решить только программным путём.
To Evgeny V. Сделал так:
{$N+} Uses CRT; var i,j: integer; EnX: array [0..63] of double; EnY: array [0..63] of double; An, Freq: double; buf: array [0..3999] of double; begin ClrScr; j:=0; An:=0; for i:=0 to 63 do begin EnX[i]:=0; EnY[i]:=0; end; for i:=0 to 3999 do buf [i]:= sin (2*pi*i*500/4000); //- опять же 1000 Гц for i:=0 to 999 do buf [i+3000]:= sin (2*pi*i*93.75/1000); //- как бы в конце появилось 750 Гц begin writeln; Freq:=750; for i:=0 to 3999 do begin j:=j mod 8000; An:=An-EnX[i mod 64]*EnX[i mod 64]-EnY[i mod 64]*EnY[i mod 64]; EnX[i mod 64]:=buf[i]*Cos(2*pi*Freq*j/8000.0); EnY[i mod 64]:=buf[i]*Sin(2*pi*Freq*j/8000.0); An:=An+EnX[i mod 64]*EnX[i mod 64]+EnY[i mod 64]*EnY[i mod 64]; Inc(j); end; writeln; writeln (An:5:3); end; end.
Всё равно An у меня принимает одно и то же значение (32.000) при любых раскладах. Оно меняется только если последние 64 элемента "буфера" заполнить нулями - тогда An=0, если последние 32, то An=16.000. "Нужную" частоту вставлял не только в конце, но и в середине (с 2000 элемента) - эффект тот-же. ???
А почитат исходники по нужной теме всегда поучительно, тем более разобраться в них! Ещё раз спасибо.

← →
Evgeny V © (2004-05-29 05:27) [12]

Извини Вадим, я и правда ввел тебя в заблуждение,
давно этим уже не занимался и подзабыл. Сорри. И так нам надо
посчитать энергию при выбранных частотах. En=Acos^2+ASin^2,
где Acos - косинусные и Asin -синусные коэффициенты (Re и Im).
Acos=Sum(buf[i])*Cos(2*pi*Freq*i/N) и Asin=Sum(buf[i])*Sin(2*pi*Freq*i/N) где i меняется oт 0 до N-1,
вот исправленный код с прорисовкой графиков сигнала и двух значений энергий для частот 1000 и 750
Графики рисую в Tchart, Series выбран Line или FastLine.

procedure TForm1.Button3Click(Sender: TObject);
var
i,j: integer;
EnX: array [0..63] of double;
EnY: array [0..63] of double;
EnX1: array [0..63] of double;
EnY1: array [0..63] of double;
Acos1000,Asin1000,Acos750,Asin750:double;
Freq1000: double;
Freq750,t1,t2,t3,t4:double;
En1000,En750:double;
buf: array [0..3999] of double;
begin
Freq1000:=1000;
Freq750:=750;
Chart1.Series[0].Clear;
Chart1.Series[1].Clear;
Chart1.Series[2].Clear;
j:=0;

for i:=0 to 63 do
begin
EnX[i]:=0;
EnY[i]:=0;
EnX1[i]:=0;
EnY1[i]:=0;
end;
for i:=0 to 999 do
buf [i]:= 0;// сигнала нет
for i:=0 to 2000 do
buf [i+1000]:= sin (2*pi*i*Freq1000/8000); //с 1000 точки сигнал 1000
for i:=0 to 999 do
buf [i+3000]:= sin (2*pi*i*Freq750/8000); //с 3000 тысячной точки сигнал 750
ACos1000:=0;
ASin1000:=0;
ACos750:=0;
ASin750:=0;

for i:=0 to 3999 do// внутри цикла у нас неявный подцикл по массивам En из 64 итераций
begin
j:=j mod 8000;
ACos1000:=ACos1000-EnX[i mod 64]; // вычитаем из значения коэффициента косинуса значения прошлого подцикла
ASin1000:=ASin1000-EnY[i mod 64];// // вычитаем из значения коэффициента синуса значения прошлого подцикла

EnX[i mod 64]:=buf[i]*Cos(2*pi*Freq1000*j/8000.0);// вычисляем значение в текущей точке
EnY[i mod 64]:=buf[i]*Sin(2*pi*Freq1000*j/8000.0);
ACos1000:=ACos1000+EnX[i mod 64];// cуммируем значения для получения значения коэффициента косинуса
ASin1000:=ASin1000+EnY[i mod 64];// cуммируем значения для получения значения коэффициента синуса
t1:=ACos1000*ACos1000;// квадрат
t2:=ASin1000*ASin1000;// квадрат
En1000:=Sqrt((t1*t1+t2*t2))/64;// находим значение энергии для заданной частоты
//корень квадратный и деление на 64 не нужно, сделано только для удобства отображения на графике,
//большие значения получаются

/ все то же для другой частоты
ACos750:=ACos750-EnX1[i mod 64];
ASin750:=ASin750-EnY1[i mod 64];
EnX1[i mod 64]:=buf[i]*Cos(2*pi*Freq750*j/8000.0);
EnY1[i mod 64]:=buf[i]*Sin(2*pi*Freq750*j/8000.0);
ACos750:=ACos750+EnX1[i mod 64];
ASin750:=ASin750+EnY1[i mod 64];
t3:=ACos750*ACos750;
t4:=ASin750*ASin750;
En750:=Sqrt((t3*t3+t4*t4))/64;
Inc(j);
Chart1.Series[0].AddY( buf[i],"",clRed);// показываем сигнал
Chart1.Series[1].AddY( En1000/50,"",clGreen);// значения энергии частоты 1000
// на 50 делю тоже только для смотрибельности на грфике
Chart1.Series[2].AddY( En750/50,"",clBlue);// значения энергии частоты 750
Edit1.Text:=IntToStr(i);// просто показываю сколько точек пробежал
if (i mod 1000) =0 then
AppLication.ProcessMessages;// даю возможность прорисоваться графику и эдиту
end;
end;
Тебе прийдется еще решать проблему с определением сигнала вообще, помнишь porog?
В линии все равно шумы есть и энергии сигнала будут присутствовать, но меньшими значениями,
поэтому их надо сравнивать с пороговым значением, его я подбирал экспериментально.
Для радио у меня получилась формула Sum(Abs(buf[i])*0.5 > Max(sqrt(En)/N) - i от 0 до N-1,
Может и тебе она подойдет, поиграйся с коэффициентом, не 0.5, а может другой.
Потом при смене частотной посылки возможны кратковременные появления энергий
"ложных" частот, так что суди о наличии частоты по нескольким опредлениям энергии,
например за 5 подциклов по 64 вычисления. Будут вопросы пиши, лучше наверное на почту.
evgeny_v@rambler.ru. Удачи, извини еще раз за ошибку в предыдущих постах.

← →
MegaVoltik (2004-05-31 10:26) [13]

А чем вас аппаратный декодер не устраивает то? Сигнал всё равно от куда то братся должен? Если качественно определять то жедаетельно АЦП ставить. А если АЦП не ставить а просто компаратор как в АОНах то я думаю для такого входного сигнала есть свои специальные методы. Сходи на http://www.telesys.ru/teleconf.shtml там помогут :)

← →
Вадим © (2004-06-01 01:34) [14]

To MegaVoltik
Спасибо, обязательно посмотрю, но у меня уже кое-что стало получаться из того, что предложил Евгений.

← →
Jeer © (2004-06-02 12:48) [15]

Еще раз - алгоритм Герцеля по приведенной ссылке.
http://factaee.elfak.ni.ac.yu/fu2k33/9mp.html
Это разновидность БПФ, ориентированная на распознавание
ограниченного числа частот в сигнале.
Имеет значительное преимущество в скорости (числе операций) перед обычным БПФ.

← →
Вадим © (2004-06-04 23:56) [16]

To Jeer: алгоритм Герцеля я посмотрел ещё когда Вы предложили его впервые (24.05.04 15:38). Но там ведь чистая математика, а разбираться в ней и, тем более, строить по этому алгоритму программу мне пока трудновато:-( Хотя, если он действительно даёт большой выигрыш по времени выполнения, то, возможно, придётся. Правда Вы сравниваете его с БПФ, а тут речь идёт о Дискретном Преобразовании, а оно тоже будет идти быстрее, чем БПФ.
В любом случае - всем спасибо за помощь.

← →
Evgeny V © (2004-06-05 06:48) [17]

> Jeer © (02.06.04 12:48) [15]

> Вадим © (04.06.04 23:56) [16]

То что я предложил,фактически и есть алгоритм Герцеля :-)) Быстрее оно только в данном случае, когда определяем весьма ограниченное количество частот из возможного спектра

← →
Anonym (2004-06-05 19:18) [18]

http://teworks.com/dtmf.htm

← →
Jeer © (2004-06-08 09:49) [19]

Вадим © (04.06.04 23:56) [16]
Рекомендую ознакомится с теорией по преобразованию Фурье, это не так сложно как кажется, заодно не будете "смешить тапочки"
по повожу отличий БПФ, ДПФ и тп.

Evgeny V © (05.06.04 06:48) [17]
>..фактически и есть алгоритм Герцеля
Не знаю, не вникал:)

>когда определяем весьма ограниченное количество частот из >возможного спектра
Это как раз "наш" случай для DTMF (DTFM)

← →
Damage (2004-06-08 20:39) [20]

2 Вадим:
Если интересует, могу на мыло кинуть прогу, которая спектр
WAV-файла рисует.
Не в реальном времени, конечно (там рукой двигать надо), но функцию полезную оттуда выцепить можно.

← →
Вадим © (2004-06-09 00:44) [21]

Jeer © (08.06.04 09:49) [19]
Ну, что ж, если я кого-нибудь насмешил, то это тоже не плохо. Хорошо, когда у кого-то поднимается настроение, даже если этот кто-то тапочки ;-)
А если серьёзно, то возможно Вы и правы, даже наверняка, и в преобразовании Фурье разобраться будет несложно, но мне сейчас не до этого (потребуется время, а есть более актуальные задачи). Тем более думаю, что теперь я смогу решить свою проблему и так (без обид!).

2 Damage:
А прога большая? И вообще это исходник или екзешник? Если екзешник, то я, конечно, не смогу "выцепить полезную функцию" оттуда. А так, посмотрел бы с удовольствием (и благодарностью).

← →
Damage (2004-06-09 13:39) [22]

Ясное дело, исходник :))

← →
Jeer © (2004-06-09 16:22) [23]

Вадим © (09.06.04 00:44) [21]

Это и плохо, т.к. вы полезли в цифровую обработку, а здесь надо иметь некоторый базис, чтобы в результате Ваша задача получила хотя бы решение.
Для занятий DSP безусловно необходимо знание математики, в том числе.
По той ссылке, что я Вам дал, вполне доступно описан алгоритм Герцеля. Если Вы не умеете переходить от формульной записи к алгоритмической, то очень странно, что Вы вообще взялись за это дело:)

Попробуем поправить вместе с Вами положение вещей, если Вы не против, конечно.

Итак, DTMF - способ тонального кодирования символов. (CCITT Recommendation Q.23 и Q.24)
Для представления одного символа используются две частоты из двух групп частот, в каждой из которых
по четыре частоты, итого: алфавит состоит из 16-ти символов.
Для декодирования символа необходимо определить мощности восьми частот и две преобладающие (упрощенно).
Могут быть использованы методы частотной фильтрации, основанные на полосовых рекурсивных или нерекурсивных фильтрах, а также методы основанные на спектральных преобразованиях, например - преобразовании Фурье (ПФ).
Жан Фурье высказал идею о том, что любой сигнал может быть представлен в виде суммы гармонических сигналов.
По сути ПФ - это способ перехода от описания сигнала во временном домене к описанию сигнала в частотном домене. Обратный способ перехода от частотного домена ко временному называется обратным ПФ (ОПФ).
Все сигналы в природе являются "аналоговыми", т.е. непрерывно изменяющимися.
Однако с появлением компьютерной техники, над естественными сигналами стала выполняться дискретизация - приведение сигнала к логической форме представления набором дискрет как во времени, так и по амплитуде.
Для таких оцифрованных (дискретных) сигналов тоже существует ПФ и называется оно ДИСКРЕТНЫМ (ДПФ).
В форме прямой реализации ДПФ требует большого числа операций (число комплексных умножений пропорционально квадрату длины выборки N*N).
Для снижения вычислительных затрат были разработаны ряд методов, получившие общее название быстрого ПФ (БПФ).
(в англо-язычной литературе Fast Fourie Transform - FFT).
Тем не менее, применение БПФ к задачам DTMF приводит к большим вычислительным затратам, т.к. вычисление ДПФ
для отдельных спектральных составляющих выполняется быстрее.
В общем случае, ДПФ работает с комплексными числами, но для случая DTMF восстановление фазовой информации неважно, поэтому можно обойтись только вещественным ДПФ.
Учитывая, что число анализируемых частот весьма мало (8) ряд ученых разработали ускоренные методы детектирования таких сигналов.
К одному из них относится метод Герцеля (Goertzel), основанный на ДПФ.
В основе лежит рекурсивное звено второго порядка.
Разностное уравнение:
Qk[n] = X[n] + Mk*Qk[n-1] - Qk[n-2] (1)
Yk[n] = Qk[n] - Wk*Qk[n-1] (2)

где
Mk = 2*Cos(2*PI*k/N)
Wk = exp(-(2*PI*k*/N)) - множители Фурье
X[n] - входной отсчет в дискретный момент времени n
Yk[n] - выходной отсчет (мощность k-ой гармоники)
Fs - частота дискретизации
N - число отсчетов.

Уравнение (1) вычисляется с частотой Fs,а уравнение (2) с частотой Fs/N, т.к. интересует только лишь
накопленное Yk[N-1] значение.
Величина Yk является комплексной, но т.к. фазовая информация не нужна, уравнение (2) для вычисления модуля трансформируется в (3)
|Yk[N-1]| = sqrt(Qk[N-1]^2 + Qk[N-2]^2 - Mk*Qk[N-1]*Qk[N-2]) (3)

Таким образом (1) и (3) являются искомым разностным уранением для вычисления мощности каждой из восьми частот.

Следующая важная задача - выбор N и Fs. Обычно Fs = 8 кГц.
Так как длительность одного сигнала должна быть не менее 40 ms, то мы имеем N = 0.04*8000 = 320 отсчетов.
Однако это приводит к выходу частот верхней группы за пределы разрешающей способности.
Поэтому необходимо решать задачу оптимального выбора N при заданном допуске различения частоты:
< 1.5% - есть частота
> 3.5 - нет частоты
Собственно, такие расчеты конечно же уже сделаны:)

Один вариантов:
N = 205

F k Mk
697 18 1.703275
770 20 1.635859
852 22 1.562297
941 24 1.482867
1209 31 1.163138
1336 34 1.008835
1477 38 0.790074
1633 42 0.559454

А вот и "сладкая пилюля" :)
Это несколько упрощенная формула.
Правда не сложно ?

function Goertzel(Buffer:array of Single; frequency, samplerate: single):single;
var
Qkn, Qkn1, Qkn2, Wkn, Mk : Single;
i : integer;
begin
Qkn:=0; Qkn1:=0; Qkn2:=0;
Wkn:=2*PI*frequency/samplerate;
Mk:=2*Cos(Wkn);
for i:=0 to Length(Buffer) do begin
Qkn2 = Qkn1; Qkn1 = Qkn;
Qkn =Buffer[i] + Mk*Qkn1 - Qkn2;
end;
Result:=(Qkn - exp(-Wkn)*Qkn1);
end;

Более точно - надо реализовать по ф-лам (1) и (3), но это домашнее задание:)

Jeer © (09.06.04 16:22) [23]
Это, конечно, здорово, но мне не совсем ясно следующее:

> Так как длительность одного сигнала должна быть не менее
> 40 ms, то мы имеем N = 0.04*8000 = 320 отсчетов.
> Однако это приводит к выходу частот верхней группы за пределы
> разрешающей способности.

Правильно ли я понял, что искомая "омега" вычисляется из отношения коэффициента k к кол-ву отсчётов N, и тогда полученная "омега" будет (или должна) отличаться от стандартной (из сетки DTFM) не более, чем на 1,5%? Но зачем вообще нужно подбирать эти k, если в формулу можно просто подставить искомую частоту и частоту дескритизации, как это сделано в примере? А если без k не обойтись, то из какого ряда их выбирать?

Jeer © (09.06.04 16:22) [23] (дополнительно)
Ну, хорошо, допустим понятно, что k показывает номер гармоники из спектральной сетки, разрешающая способность которой зависит от N. И здесь уже "насмешили тапочки" именно Вы, т.к. в описании алгоритма Герцеля, на которое Вы же дали ссылку, говорится примерно следующее: "значение N должно выбираться исходя из двух противоречивых условий: во-первых, оно должно быть достаточно малым, чтобы успеть различить кратковременные сигналы; во-вторых, достаточно большим, для обеспечения разрешающей способности, т.к. чем больше N, тем лучше эта самая способность". Таким образом, очевидно, что "за пределы разрешающей способности" выходят частоты нижней группы при малых N!
Теперь на счёт "сладкой пилюли". Во-первых даже я:-) заметил там грубую ошибку. Вместо:
for i:=0 to Length(Buffer) do begin
должно быть
for i:=0 to Length(Buffer)-1 do begin
думаю понятно почему. Во-вторых зависимость Qkn от времени носит синусоидальный, а не накопительный характер. В итоге Result имеет непредсказуемое значение, зависящее от фазы сигнала в последних ячейках буфера. Как же тогда выявлять сигнал? В-третьих мне до сих пор не ясно, почему мы обязаны искать частоты именно из частотной сетки, а не те, что нужны? Ведь в алгоритме, предложенном Евгением, искомые частоты так же "не попадают" в сетку и, тем не менее, прекрасно определяются (размер буфера (N) - 256). В то же время, если искать "по сетке", то действительно, получается не очень хорошо. Кстати ещё одно странно для меня - везде написано, что для преобразования Фурье размер буфера должен быть степенью двойки, а в алгоритме Герцеля совсем не обязательно...
Ну, что, теперь я прав, или уже вместе с тапочками и носки засмеялись? :-)

Ну, вот вы уже начали думать, это хороший признак.:)

>Ну, хорошо, допустим понятно, что k показывает номер гармоники
Без всяких "ну" и "допустим понятно" - это так и есть.

>Таким образом, очевидно, что "за пределы разрешающей способности"
выходят частоты нижней группы при малых N!

Вы не о том.
При 320 отсчетах ширина спектральной линии dF = 8000/320 = 25 Гц.
В соответствии с "Рекомендациями.." частота может "плавать" от среднего значения до 1.5%, а это дает на наивысшей частоте 1633 Гц возможное отклонение до +/- 25 Гц, что превышает ширину спектральной линии.
Т.е. искомая частота может быть идентифицирована как неверная.
Одним из вариантов, является снижение длины выборки, что расширяет спектральную линию.
В примере приведено N=205
Это дает ширину линии 8000/205 = 39 Гц. Наложение на соседнюю линию все равно есть, но это компромисс.

>Во-первых даже я:-) заметил там грубую ошибку
Это похвально:)
Там еще пара ошибок и тоже "грубых".

>Во-вторых зависимость Qkn от времени носит синусоидальный, а не накопительный характер.
Синусоидальный линейно нарастающий, а это и есть накопление.

>В итоге Result имеет непредсказуемое значение, зависящее от фазы сигнала в последних ячейках буфера
Это у Вас:)

>Как же тогда выявлять сигнал?
Если бы Вы просканировали по частоте от искомой и, например, выше, то увидели бы, что отклик функции Герцеля имеет единственный максимум в одной полуплоскости и много максимумов в другой, причем зависимость носит синусоидальный характер с экпоненциальным затуханием. Этот единственный максимум находится на искомой частоте.
Отсюда понятно, что и как искать.

>В-третьих мне до сих пор не ясно, почему мы обязаны искать частоты именно из частотной сетки, а не те, что нужны?
Потому, что Вы до сих пор не поняли главного: это ДИСКРЕТНОЕ преобразование Фурье и методически правильно
работать по сетке дискретных частот.
Можно представить себе ДПФ как набор из N узкополосных фильтров, однако частотная характеристика их не идеально прямоугольная, а с перекрытием. Поэтому частоты не попадающие точно в сетку тоже, естественно, будут "звенеть" на смежных гребенках.

>Кстати ещё одно странно для меня - везде написано, что для преобразования Фурье размер буфера должен быть степенью двойки, а в алгоритме Герцеля совсем не обязательно...

Это "странно" только для Вас, т.к. Вы еще только начинаете в данном вопросе "барахтаться".
ДПФ не привязано к числу точек никак, с другой стороны, действительно большинство РЕАЛИЗАЦИЙ БПФ чаще всего привязано к степени 2,а почему - читайте литературу или сами догадайтесь.
Могу еще добавить что есть быстрые преобразования основанные на числах Мерсенна, Фибоначчи.
Основная операция "свертка" весьма эффективна с их использованием.
Есть, наконец, другие ортогональные базисные функции (несинус), например Уолша, где тоже имеет место быть быстрое преобразование дискретных функций.

Так, что "тапочки отсмеялись, а носочки мои еще посмеиваются" :)

Jeer © (15.06.04 16:38) [26]

> >В итоге Result имеет непредсказуемое значение, зависящее от фазы сигнала в последних ячейках буфера
Это у Вас:)

Почему же только у меня?! Я делал по образцу, предложенному Вами. И, так как Result вычисляется из Qkn и Qkn1, которые в свою очередь являются расчитанными из двух последних ячеек массива Buffer, то, и у Вас, и у других, кто повторит этот алгоритм, результат будет зависеть от ДВУХ ПОСЛЕДНИХ ЯЧЕЕК БУФЕРА! Я, конечно, ценю Вашу попытку мне помочь и чему-то научить, но, в таком случае уж будьте добры не зазнавайтесь и дайте правильную информацию.

Jeer © (15.06.04 16:38) [26]

Вы ошибаетесь насчет двух последних ячеек Buffer.
Результат определяется по двух последним отсчетам полученным
в результате НАКОПЛЕНИЯ данных в ячейках (переменных) Qk и Qkn1.
Так, что здесь - все нормально.
Я уже упомянул, что стоит рассматривать более корректную реализацию по ф-ле (3). Там все однозначно.
Если интересно - могу выложить тестовую утилитку по DTMF.
Там все нагладно видно.
Там использованы оба варианта определения спектра - по формулам (2) и (3).
Кроме того, для улучшения качества спектрального анализа (устранения явления Гиббса) используется "оконное" взвешивание с различными окнами (Хеннинг, Каппелини, Ханн и тп)
В этой утилите есть возвожность предобработки серии с помощью
десятка различных окон - результат хорошо заметен.

Тестовый сигнал символа "1" ( частоты 697 и 1209 Гц )
http://algcom.nm.ru/snap_01.png
Спектр по Герцелю по ф-ле (3)
http://algcom.nm.ru/snap_02.png
Предварительная обработка окном Blakman
http://algcom.nm.ru/snap_03.png
Спектр обработанного сигнала
http://algcom.nm.ru/snap_03.png

Jeer © (30.06.04 19:56) [28]

> Если интересно - могу выложить тестовую утилитку по DTMF

Ещё бы, если это исходник конечно интересно!

To Jeer. Так как на счёт утилитки?

Завтра-послезавтра выложу.
Заскочил на минутку.
Насчет исходников - не обессудь, не будет. Не в моих правилах подавать ложечку в ротик, да еще и пережевывать.
Хочешь учиться - помогу, но не более.
Есть такой необходимый, но недостаточный параметр (да простят меня Дамы) - "жопочасы".
Ты их должен наработать:)

http://algcom.nm.ru/dtmf.zip

Тест DTMF-методов.
В этой версии оставлены только два (простейший и классический Герцель-методы).

Декодирование DTFM Найти похожие ветки