18+ Некоторые материалы сайта могут содержать информацию, запрещенную для детей.

Запись от AZM на субдомене electronics-and-mechanics
Все записи на субдомене: Электроника и механика (записки от AZM)

Для микроконтроллеров и не только: Алгоритм Герцеля, БПФ, и другие цифровые способы определить наличие искомой частоты в оцифрованном сигнале (цифровые фильтры)
Сразу скажу - я не математик, более того, я только вчера узнал что: sin(x)^2+cos(x)^2=1, потому что благополучно прогуливал уроки в школе.
В общем, я не буду сыпать формулами, а расскажу всё как ни будь по проще, постараюсь привести наглядные примеры.

Алгоритмы детекторов наличия частоты в оцифрованном сигнале
Если в общих чертах, то что бы определить наличие некой частоты в оцифрованном сигнале нужно:
Массив_синусов в одну волну
такой как если бы мы создавали синусоиду искомой частоты на нашей частоте дискретизации.
Массив_косинусов в одну волну
такой, как если бы мы создавали синусоиду искомой частоты на нашей частоте дискретизации.
Естественно, что Массив_синусов и Массив_косинусов будут иметь одну длину, а массив с данными нашего оцифрованного сигнала должен быть не меньшей длины.
Значения в массивах должны быть знаковыми (не важно "signed char" или "signed int", важно что "signed").
Затем нужно проделать следующее:
a=0;
b=0;
for (n=0; n < length(Массив_синусов); n++){
a = a + Оцифровка[n] * Массив_синусов[n];
b = b + Оцифровка[n] * Массив_косинусов[n];
}
Теперь получаем амплитуду искомой частоты в сигнале:
Амплитуда = sqrt(a*a + b*b)/length(Массив_синусов)/(Максимальная_амплитуда_оцифровки/2);

Где "Максимальная_амплитуда_оцифровки": для 8 бит будет 128, для 16 бит 32768.

Естественно, что можно эту процедуру в одном цикле проделать для множества частот сразу, для этого потребуются наборы массивов синусов и косинусов для всех искомых частот и по массиву a[] и b[], как то так:
for (m=0; m <число_частот; m++){
for (n=0; n < длина_Массива_синусов; n++){
a[m] = a[m] + Оцифровка[n] * Массив_синусов[m][n];
b[m] = b[m] + Оцифровка[n] * Массив_косинусов[m][n];
}
}
for (m=0; m <число_частот; m++){
Амплитуда[m] = sqrt(a[m]*a[m] + b[m]*b[m]);
}

Длины массивов косинусов и синусов должны быть равны и их длина будет такой, какой длины будет массив синусов и косинусов для самой нижней искомой частоты, то есть для анализа нескольких частот таблицы для высоких частот будут фактически длиннее чем необходимо.

Для достижения большей точности анализа, нужно повышать число проходов, то есть выбирать участок сигнала для анализа не 1 период колебания, а 10 или более (таблицы синусов и косинусов естественно при этом не нужно делать длиннее, просто проходить по ним снова и снова).

Для таблиц синусов и косинусов можно вообще пойти на ухищрения, а именно - если посмотреть на график:
График sin(x) и cos(x)

очевидно, что косинус, это синус сдвинутый на 1/4 периода частоты, а значит достаточно иметь только таблицу синуса длиннее на 1/4 периода и начинать проход по ней для получения косинуса сдвинув указатель на 1/4 периода от начала.

Как сгенерировать массив значений синуса для заданной частоты дискретизации:
// Например, для частоты 88.5 Гц и частоты дискретизации 11025
#define GeneratorK 0.008027210884353741 // = Нужная_частота/Частота_дискретизации
#define GeneratorAmplitude 128 // если оцифровка 16 бит, то здесь 32768
#define GeneratorLenOnePeriod 125 // = Частота_дискретизации / Нужная_частота

int n;
signed char ArraySin[GeneratorLenOnePeriod];
for(n=0; n<GeneratorLenOnePeriod; n++){
ArraySin[n]=(signed char)(GeneratorAmplitude * sin(2*PI * n * GeneratorK));
}
Как сгенерировать массив значений синуса заданой длины и максимальных значений
Исходный код на СИ (C):
// Например, для 127 в качестве максимума и длиной 64 элемента
#define GeneratorMaxAmp 127
#define GeneratorLenArray 64

int n;
signed char ArraySin[GeneratorLenArray];
for(n=0; n<GeneratorLenArray; n++){
ArraySin[n]=(signed char)(GeneratorAmplitude * sin(2*PI * n * (1/GeneratorLenArray)));
}
Исходный код на perl:
#!/usr/bin/perl
$PI=3.141592653589793;

$GeneratorMaxAmp=127;
$GeneratorLenArray=64;

for($n=0; $n<$GeneratorLenArray; $n++){
print int($GeneratorMaxAmp * sin(2 * $PI * $n * (1/$GeneratorLenArray ))).", ";
}
print "\n";
exit;
Теперь самое интересное:
этот метод обнаружения частоты работает даже если у вас дискретность всего 1 бит на выборку!
Проще говоря, если вам не нужно точно измерять соотношение амплитуд сигналов в оцифровке, а достаточно лишь ограничиться результатом вида "есть/нет" то таблицу синусов/косинусов можно сделать однобитную и входной сигнал принимать не с АЦП (ADC), а с обычного компаратора или даже с логического входа микроконтроллера.
При однобитном сигнале и таблице синус/косинус можно вместо умножения применять XOR, что значительно ускорит алгоритм.

Алгоритм Герцеля
Говорят, что алгоритм Герцеля это частный случай от того, что изложено выше, просто оптимизированный для одной частоты, но для меня это магия, по этому просто приведу реально работающий код (наверно это не лучшая реализация, но рабочая):
float constantaA = 2 * cos(PI * 2 * (нужная частота / частота дискретизации));
float TMPVAR_v0=0;
float TMPVAR_v1=0;
float TMPVAR_v2=0;
float resultat;
float prev_resultat=0;
bool flag_frq_yes=1;
for (int n=0; n<100; n++){
TMPVAR_v0 = (InSample[n] + constantaA * TMPVAR_v1 - TMPVAR_v2);
TMPVAR_v2=TMPVAR_v1;
TMPVAR_v1=TMPVAR_v0;
resultat = (TMPVAR_v1*TMPVAR_v1 + TMPVAR_v2*TMPVAR_v2 - constantaA * TMPVAR_v1 * TMPVAR_v2);
if (prev_resultat < resultat){prev_resultat=resultat*0.85;}else{flag_frq_yes=0}
}
соответственно на выходе из цикла "flag_frq_yes" либо 1 если в массиве оцифровки "InSample" есть сигнал с нужной частотой или 0, если такового сигнала там нет.
Алгоритм работает так, что значение "resultat" при наличии в оцифровке нужной частоты, при каждом новом проходе будет больше предыдущего значения "resultat", иначе оно будет изменяться от шага к шагу (например в случае если в оцифровке гармоника нужной частоты).
Допустимые колебания заданы как 0.85, чем меньше допуск, тем более точна должна быть искомая частота и тем меньше шума должно быть в исходном сигнале.
Длина "окна" здесь выбрана 100 (алгоритм ориентирован на поиск частот порядка 1500Гц при оцифровке 15кГц), при большем окне точность будет выше, ширину более 200 наверно вообще нет смысла выбирать для таких приложений как декодирование DTMF или CTCSS при оцифровке 15кГц или ниже.

Алгоритм ФНЧ (фильтра низких частот)
ФНЧ и реализуется крайне просто: надеюсь, все знают, что максимальная частота, которая может быть оцифрована, равна частоте дискретизации / 2.
Соответственно если у нас есть оцифровка на 12000 Гц, а нам надо отрезать всё, что выше 300Гц, то всё что надо сделать - преобразовать частоту дискретизации из 12000 в 600.
То есть нужно вычислить среднее арифметическое 20 выборок и это значение использовать для этих 20 выборок, так как:
12000 / 600 = 20
Пример в коде:
define DownsampleK 20 // = текущая_большая_частота/нужная_меньшая; (пусть условно нам надо опустить частоту семплирования в 20 раз)
define InSampleLeng 1700 // пусть условно у нас 1700 выборок

signed char ArrayInSample[InSampleLeng]=(наши,результаты,оцифровки);
signed int tmpval;
for (int n=0; n<(InSampleLeng-DownsampleK); n++){
tmpval=0;
for (int m=0; m<DownsampleK; m++){tmpval=tmpval+(signed int)ArrayInSample[n+m];}
ArrayInSample[n]=(signed char)(tmpval/DownsampleK)
}
Результат работы такого алгоритма будет аналогичен результату действия интегрирующей RC цепи, если требуется алгоритм с более крутой АЧХ (Блэкмена, Кайзера, Хемминга, Хеннинга), то придётся применить более сложный алгоритм требующий больше ресурсов, построение такого алгоритма описано например здесь:
habrahabr.ru - Построение цифрового фильтра с конечной импульсной характеристикой

Алгоритм ФВЧ (фильтра высоких частот)
По сути, обратный алгоритм к алгоритму ФНЧ фильтра низких частот.
То есть нам нужно проделать тот же процесс что и для реализации ФНЧ, только полученные значения вычесть из наших данных оцифровки.

Полосовой фильтр
Комбинация ФНЧ и ФВЧ, просто ФНЧ - определяет верхнюю частоту которую пропустит фильтр, а ФВЧ - нижнюю частоту.
По иному говоря:
ФНЧ с частотой среза 4000Гц, затем ФВЧ с частотой среза 500Гц = фильтр с полосой 500 ... 4000 Гц

Изложенное здесь написано для тех, кто никак не мог понять, как же обрабатывать оцифровку, а что-то сделать с оцифровкой ему нужно было "ещё вчера".
Алгоритмы, приведённые здесь вероятно кривые, медленные, но они работают и приведены в качестве примера. Улучшать, оптимизировать, вникнуть в суть ещё глубже, каждый сможет сам, по мере надобности.

Что ещё почитать по теме:
Программное дeкодирование DTMF по принципу АОН, на базе микроконтроллера PIC16F628
DEMO.DESIGN - Subj : быстpая частотная обpаботка сигнала

Задать вопрос или оставить свой комментарий можно на форуме:
Форум на 27kb.ru - Гражданская-радиосвязь.РФ

Добавлено: 3963 дн 17 час 20 мин 48 сек назад | Внесений правок: 6 | Последняя правка: 3365 дн 0 час 45 мин 1 сек назад



Электроника и механика (записки от AZM)