Определение передаваемых частот Опции

[Kseniia]

Всем доброго времени суток!

Задача: есть DTMF метка, необходимо определить частоты которые передавались, с точностью до 3го знака после запятой. (DTMF - передается сумма двух гармоник, НЧ и ВЧ)

Что проделано:
выделяю из сигнала DTMF, fd=44100 Гц, для того, чтобы N была кратна двум, и равнялась 2048, немного обрезаю сигнал до 0.046 секунд (было 0.05 сек), таким образом T=1/fd=2.268*10^5.
Беру БПФ, строю график зависимости уровня(?)(20 log(sqwr(ReБПФ^2+ImБПФ^2)/fd-N), дБ от частоты, Гц (прикрепляю скрин, freqn=n*fd/N, n=0,...,N-1).

Пики приходятся на 947,46Гц и 1205,9Гц, хотя должны быть 941Гц и 1209Гц. такая разница по видимому из-за огромного ШАГА по частоте, в 21,55 Гц, то есть передаваемые частоты не кратны нашему "шагу".

Что делать дальше - не понимаю.
Пробовала инетрполяцию Лагранжа - чушь полная, не получается.
Есть предложение использовать коэффициенты корреляции, относительно синусоидального и косинусоидального сигналов, но результатов пока не дал, да и не до конца понятно, что с этими коэф делать.

вопрос: Как математически можно определить частоты?

[Kseniia]

Размер выборки кстати ограничен? Или можно копить до посинения?

размер выборки ограниченг самим сигналом, DTMF по стандарту идут от 40 мс сам сигнал и от 40 мс пауза между другим сигналом. В реальности же одна метка - 40-50 мс+перерыв 40-50мс.

вместо быстрого преобразования, посчитайте обычные интегралы Фурье s=\int{f(x)*sin(x)}; c=\int{f(x)*cos(x)}, на частоте 941 +- пара Гц с мелким шагом в 0.1Гц например чтобы не 3 в десяток-другой точек получился, на результат (s^2+c^2) натяните параболу ax^2+bx+c максимум будет в -b/2a.
ну или то же самое с автокорреляционной функцией можно попробовать.
исходные данные сюда выложить можете?

исходные данные:  DTMF.txt ( 37.08 килобайт ) Кол-во скачиваний: 106
, fd=44100 Гц

а на счет обычных интегралов Фурье, непонятно что туда подставлять, если у меня есть только отсчеты, что за f(x) и sin (x), что за x?
я параболическую интерполяцию применила, результат близок к истине, но не истина.

на счет автокоореляционной функции - можно поподробнее?
мне предлагают корреляционную функцию использовать, посчитать Rs=СУМ(sin(2pifn)*s(n), n=0..N и Rc=СУМ(cos(2pifn)*s(n). не непонятно становится, первое, у меня 2 частоты, значит не sin надо брать а сумму двух sin от ВЧ и НЧ. Я посчитала Rs и Rc, далее R=sqr(Rs^2+Rc^2), и ничего у меня не получается, да и не понятно, даже если я найду этот коэф корреляции, как мне это поможет найти искомую частоту? по идее я должна КАК-ТО построить график R от f и f взять близлежащие к искомой частоте, но! в формулах Rs и Rc присутсвует N, а оно определено и опять я кручусь на том, что выборка ограничена.

Подойдёт если всегда а приори известно в каком узком диапазоне будут лежать интересующие частоты. Если нет, то можно также дополнить выборку нулями и сделать БПФ, получится аналог интерполяции между частотными отсчётами. Но все эти операции размазывают пики. Есть гарантия что при таком подходе не будет ложного выделения пика частоты например в условиях шума? Ведь придётся выделять пик из группы отсчётов с отличием по амплитуде << 3 дБ. Как ни крути а разрешающая способность Фурье (любого, в т.ч. обычного ДПФ) ограничена только размерностью входной выборки и никакие операции и ухищрения её не поднимут. Сверхразрешение же даёт всегда выраженные однозначные пики, но они могут сместиться из-за шума/малого размера выборок/криворукости. Кароче смотреть надо Но до третьего знака после запятой это сильно. Обычно речь идёт о Гц, ну максимум десятых долях Герца. Для таких точностей потребуется большой сигнал-шум и выборки порядка не менее 1е4 отсчётов. Это даже для сверхразрешения. Для Фурье будут нужны миллионы отсчётов.

впринципе известно в каком диапазоне лежат частоты, можно вначале определить грубо частоты (по стандарту есть 4 НЧ и 4 ВЧ) и уже после грубого определения начинать ухищряться точность определить.
на счет шума - это следующая задача, помехоустойчивость метода определения частот проверить.

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

я делала: рассчет окна, перемножение окна с сигналом, БПФ от получившегося сигнала. и уже параболическая интерполяция. Вы предлагаете взять логарифм от |БПФ| и после этого аппрокимировать? правильно?

В самом начале рассуждений вы назначили что будет бпф. Почему?
Интереса ради, попробуйте взять интервал для анализа кратный не 2^н, а кратный периоду искомой частоты: результат вас порадует )

у меня 2 частоты, я не назначала, я просто предположила, что для достижения поставленной цели можно применить БПФ, если есть другие методы - с удовольствием послушаю, повторюсь: я только начинаю так плотно работать с сигналами и много еще не знаю...

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

прологарифмировала, результат порадовал, спасибо.
на сечт коэф, у нас же амплитуда меняется, как я могу пропорционально ей или ее степени выбрать? я поперебирала, и нашла что 0.1 - ближе всего результат дает.

[_pv]

размер выборки ограниченг самим сигналом, DTMF по стандарту идут от 40 мс сам сигнал и от 40 мс пауза между другим сигналом. В реальности же одна метка - 40-50 мс+перерыв 40-50мс.
исходные данные:  DTMF.txt ( 37.08 килобайт ) Кол-во скачиваний: 106
, fd=44100 Гц

а на счет обычных интегралов Фурье, непонятно что туда подставлять, если у меня есть только отсчеты, что за f(x) и sin (x), что за x?
я параболическую интерполяцию применила, результат близок к истине, но не истина.

на счет автокоореляционной функции - можно поподробнее?

впринципе известно в каком диапазоне лежат частоты, можно вначале определить грубо частоты (по стандарту есть 4 НЧ и 4 ВЧ) и уже после грубого определения начинать ухищряться точность определить.
на счет шума - это следующая задача, помехоустойчивость метода определения частот проверить.

про автокорреляцию это я погорячился, не поможет она тут.
а для фурье выборка пожалуй коротковата чтобы 0.001Гц вытянуть за 40мс.

похоже надо грубо определять частоту через БПФ или как обычно в декодерах герцелем чтобы просто понять какая цифра передаётся и знать частоту +-10Гц, а потом просто подгонять коэффициенты у суммы двух синусов наименьшими квадратами например.

Код

d = Import["d:\\downloads\\dtmf.dat"];
fit = FindFit[d, a0*Sin[2*Pi*f0*t + p0] + a1*Sin[2*Pi*f1*t + p1], {{a0, 1}, {f0, 940}, {p0, 0}, {a1, 1}, {f1, 1200}, {p1, 0}}, t]
{a0 -> 0.305163, f0 -> 941., p0 -> 1.76526*10^-6, a1 -> 0.305164, f1 -> 1209., p1 -> 7.98036*10^-8}