Автоматическая передача данных через радиоканал

[ru0aog]

У нас расстояние между разделяемыми сигналами - 170 Гц, значит, хорошо бы установить ширину фильтра 170/2 = 85 Гц.
Примем частоту дискретизации вдвое выше полосы пропускания тракта - 6 кГц. Кол-во отсчётов - 132. Длительность кванта - 22 мс.
Обрабатываем 22 периода сигнала 1 кГц.
.

.
ширина характеристики - 70 Гц
.

.
спектр поближе,
селективность 65/15 = 12 дБ
.

.
Будто бы неплохо. Но длительность одного кванта занимает всё время посылки.
А если поднять частоту сэмплирования вдвое, до 12 кГц, сохранив количество отсчётов 132?
.

.
Нет, ширина характеристики сразу расползается. Ведь вместо 22 периодов сигнала теперь осталось только 11.

.

Гёрцель_05b.rar

(137.51 КБ) 275 скачиваний

[ru0aog]

С другой стороны, если не менять параметры оцифровки, а двигать только частоту настройки, то ширина характеристики и селективность почти не меняются.
Например, 10 кГц сэмплирование, 73 отсчёта (длительность чтения - 7,3 мс) - настройка 1 кГц
.

.
то же самое, но настройка 2 кГц
.

.
Хм, если же поднять частоту семплирования до 50 кГц и количество отсчётов до 367 (длительность чтения - 7,3 мс), то ширина характеристики не изменится, селективность та же...
.

.
Собственно говоря, спектральное разрешение равно отношению частоты дискретизации к количеству отсчётов.
В обоих случаях оно одинаково:
10кГц/73 = 136,9 Гц
50кГц/367 = 136,2 Гц,
как и длительность чтения:
1/10кГц = 100 мкс, 100*73 = 7,3 мс,
1/50кГц = 20 мкс, 20*367 = 7,34 мс.

Выходит, за 7,3 мс не получить ширину фильтра 85 Гц...

[ru0aog]

Хорошо, пойдём в другую сторону.
Теорема Котельникова ставит нам ограничение: частота сэмплирования не должна быть меньше, чем удвоенное значение наивысшей частоты обрабатываемого сигнала.
Пусть частота сэмплирования будет 6 кГц. При длительности чтения 7,3 мс это 37 отсчётов.
Смотрим.
.

.
Не обманул Котельников - выше половины частоты дискретизации получается зеркальное отражение характеристики.
Посмотрим ширину - нет изменений.
.

[ru0aog]

Продолжаем.
Без использования оконных функций мы получали от силы 13 дБ соотношения сигнал/шум после обработки цифровым фильтром.
Очень мешали выбросы по краям основного пика характеристики
.

.
Для их устранения придумали различные оконные функции.

Как это выглядит?
Без оконной функции осциллограмма сигнала была вот такая
.

.
Оконная функция выполняется над таблицей отсчётов ДО вычисления алгоритмом Гёрцеля,
и вот осциллограмма после её работы (обработано функцией HAMMING)
.

.
Такая осциллограмма чем-то напоминает классическую колоколообразную телеграфную посылку.
Собственно, так оно и есть. Спектр обработки таких данных практически теряет боковые всплески
.

.
Подавление отстоящих частот составляет более 40 дБ.
За это приходится заплатить расширением характеристики и падением уровня основного пика - примерно вдвое.
И уровень упал вдвое, и характеристика расширилась тоже вдвое.
.

[ru0aog]

Пора переходить к ручным расчётам, чтобы прочувствовать всю прелесть работы алгоритма

Исходные данные:
Частота дискретизации Fs = 4 кГц,
Частота настройки фильтра Fq = 1 кГц,
Количество отсчётов N = 20.

Максимальная рабочая частота - Fmax = Fs/2 = 2 кГц
Спектральное разрешение Fs/N = 200 Гц

Пусть чтение будет на пределе разрешения 10-разрядного АЦП: 0...1023
.

.
Вычисляем коэффициенты:
k = N * Fq / Fs = 20*1000/4000 = 5,
omega = (2 * Pi * k) / N = 2*3,14159*5/20 = 1,571
sine = Sin(omega) = Sin(1,571) = 1
cosine = Cos(omega) = Cos(1,571) = 0
coeff = 2 * cosine = 2*0 = 0
q0 = 0, q1 = 0, q2 = 0

Как хорошо, что k получился целочисленный и синус омеги = 1 (в радианах, естественно)!
Даже на осциллограмме видно, что точки отсчётов попали на пики графика.
.

.
Считаем в цикле от 0 до N:
0:
q0 = coeff * q1 - q2 + A(i) = 0*0 - 0 + 512 = 512,
q2 = q1 = 0
q1 = q0 = 512

1:
q0 = coeff * q1 - q2 + A(i) = 0*512 - 0 + 1023 = 1023,
q2 = q1 = 512
q1 = q0 = 1023

2:
q0 = coeff * q1 - q2 + A(i) = 0*1023 - 512 + 512 = 0,
q2 = q1 = 1023
q1 = q0 = 0

3:
q0 = coeff * q1 - q2 + A(i) = 0*0 - 1023 + 0 = -1023,
q2 = q1 = 0
q1 = q0 = -1023

4:
q0 = coeff * q1 - q2 + A(i) = 0*(-1023) - 0 + 512 = 512,
q2 = q1 = -1023
q1 = q0 = 512

5:
q0 = coeff * q1 - q2 + A(i) = 0*512 - (-1023) + 1023 = 2 046,
q2 = q1 = 512
q1 = q0 = 2046

6:
q0 = coeff * q1 - q2 + A(i) = 0*2046 - 512 + 512 = 0,
q2 = q1 = 2046
q1 = q0 = 0

7:
q0 = coeff * q1 - q2 + A(i) = 0*0 - 2046 + 0 = -2046,
q2 = q1 = 0
q1 = q0 = -2046

8:
q0 = coeff * q1 - q2 + A(i) = 0*(-2046) - (-2046) + 512 = -1534,
q2 = q1 = -2046
q1 = q0 = -1534

9:
q0 = coeff * q1 - q2 + A(i) = 0*(-1534) - (-2046) + 1023 = 3069,
q2 = q1 = -1534
q1 = q0 = 3069

10:
q0 = coeff * q1 - q2 + A(i) = 0*3069 - (-1534) + 512 = 2046,
q2 = q1 = 3069
q1 = q0 = 2046

11:
q0 = coeff * q1 - q2 + A(i) = 0*2046 - 3069 + 0 = -3069,
q2 = q1 = 2046
q1 = q0 = -3069

12:
q0 = coeff * q1 - q2 + A(i) = 0*(-3069) - 2046 + 512 = -1534,
q2 = q1 = -3069
q1 = q0 = -1534

13:
q0 = coeff * q1 - q2 + A(i) = 0*(-1534) - (-3069) + 1023 = 4092,
q2 = q1 = -1534
q1 = q0 = 4092

14:
q0 = coeff * q1 - q2 + A(i) = 0*4092 - (-1534) + 512 = 2046,
q2 = q1 = 4092
q1 = q0 = 2046

15:
q0 = coeff * q1 - q2 + A(i) = 0*2046 - 4092 + 0 = -4092,
q2 = q1 = 2046
q1 = q0 = -4092

16:
q0 = coeff * q1 - q2 + A(i) = 0*(-4092) - 2046 + 512 = -1534,
q2 = q1 = -4092
q1 = q0 = -1534

17:
q0 = coeff * q1 - q2 + A(i) = 0*(-1534) - (-4092) + 1023 = 5115,
q2 = q1 = -1534
q1 = q0 = 5115

18:
q0 = coeff * q1 - q2 + A(i) = 0*5115 - (-1534) + 512 = 2046,
q2 = q1 = 5115
q1 = q0 = 2046

19:
q0 = coeff * q1 - q2 + A(i) = 0*2046 - 5115 + 0 = -5115,
q2 = q1 = 2046
q1 = q0 = -5115

20:
q0 = coeff * q1 - q2 + A(i) = 0*(-5115) - 2046 + 512 = -1534,
q2 = q1 = -5115
q1 = q0 = -1534

Реальная часть:
real = q1 - q2 * cosine = -1534 - (-5115)*0 = -1534

Мнимая часть:
imag = q2 * sine = -5115 * 1 = -5115

Магнитуда:
magn = Sqr(real * real + imag * imag) = корень((-1534)^2+(-5115)^2)) = корень(2353156+26163225) = 5340

Сверяем. Однако - сошлось
.

.

Гёрцель_05d.rar

(162.52 КБ) 1055 скачиваний

[ru0aog]

И для контроля, проверим таким же способом таблицу отсчётов, но для частоты 800 Гц:
.

.
Таблица отсчётов сигнала 800 Гц:
.

.
Настройка фильтра не изменилась, поэтому коэффициенты те же:
k = N * Fq / Fs = 20*1000/4000 = 5,
omega = (2 * Pi * k) / N = 2*3,14159*5/20 = 1,571
sine = Sin(omega) = Sin(1,571) = 1
cosine = Cos(omega) = Cos(1,571) = 0
coeff = 2 * cosine = 2*0 = 0
q0 = 0, q1 = 0, q2 = 0

Считаем в цикле от 0 до N:
0:
q0 = coeff * q1 - q2 + A(i) = 0*0 - 0 + 512 = 512,
q2 = q1 = 0
q1 = q0 = 512

1:
q0 = coeff * q1 - q2 + A(i) = 0*512 - 0 + 998 = 998,
q2 = q1 = 512
q1 = q0 = 998

2:
q0 = coeff * q1 - q2 + A(i) = 0*998 - 512 + 813 = 301,
q2 = q1 = 998
q1 = q0 = 301

3:
q0 = coeff * q1 - q2 + A(i) = 0*301 - 998 + 211 = -787,
q2 = q1 = 301
q1 = q0 = -787

4:
q0 = coeff * q1 - q2 + A(i) = 0*301 - 301 + 24 = -277,
q2 = q1 = -787
q1 = q0 = -277

5:
q0 = coeff * q1 - q2 + A(i) = 0*(-277) - (-787) + 510 = 1297,
q2 = q1 = -277
q1 = q0 = 1297

6:
q0 = coeff * q1 - q2 + A(i) = 0*1297 - (-277) + 998 = 1275,
q2 = q1 = 1297
q1 = q0 = 1275

7:
q0 = coeff * q1 - q2 + A(i) = 0*1275 - 1297 + 814 = -483,
q2 = q1 = 1275
q1 = q0 = -483

8:
q0 = coeff * q1 - q2 + A(i) = 0*(-483) - 1275 + 213 = -1062,
q2 = q1 = -483
q1 = q0 = -1062

9:
q0 = coeff * q1 - q2 + A(i) = 0*(-1062) - (-483) + 24 = 507,
q2 = q1 = -1062
q1 = q0 = 507

10:
q0 = coeff * q1 - q2 + A(i) = 0*507 - (-1062) + 508 = 1570,
q2 = q1 = 507
q1 = q0 = 1570

11:
q0 = coeff * q1 - q2 + A(i) = 0*1570 - 507 + 997 = 490,
q2 = q1 = 1570
q1 = q0 = 490

12:
q0 = coeff * q1 - q2 + A(i) = 0*490 - 1570 + 816 = -754,
q2 = q1 = 490
q1 = q0 = -754

13:
q0 = coeff * q1 - q2 + A(i) = 0*(-754) - 490 + 214 = -276,
q2 = q1 = -754
q1 = q0 = -276

14:
q0 = coeff * q1 - q2 + A(i) = 0*(-276) - (-754) + 23 = 777,
q2 = q1 = -276
q1 = q0 = 777

15:
q0 = coeff * q1 - q2 + A(i) = 0*777 - (-276) + 507 = 783,
q2 = q1 = 777
q1 = q0 = 783

16:
q0 = coeff * q1 - q2 + A(i) = 0*783 - 777 + 997 = 220,
q2 = q1 = 783
q1 = q0 = 220

17:
q0 = coeff * q1 - q2 + A(i) = 0*220 - 783 + 817 = 34,
q2 = q1 = 220
q1 = q0 = 34

18:
q0 = coeff * q1 - q2 + A(i) = 0*34 - 220 + 215 = -5,
q2 = q1 = 34
q1 = q0 = -5

19:
q0 = coeff * q1 - q2 + A(i) = 0*(-5) - 34 + 23 = -11,
q2 = q1 = -5
q1 = q0 = -11

20:
q0 = coeff * q1 - q2 + A(i) = 0*(-11) - (-5) + 505 = 510,
q2 = q1 = -11
q1 = q0 = 510

Реальная часть:
real = q1 - q2 * cosine = 510 - (-11)*0 = 510

Мнимая часть:
imag = q2 * sine = -11 * 1 = -11

Магнитуда:
magn = Sqr(real * real + imag * imag) = корень(510^2+(-11)^2)) = 510

Вот и результат:
Сигнал полного размаха с частотой 800 Гц показал магнитуду 510,
сигнал того же размаха, но на частоте настройки 1000 Гц дал магнитуду 5340.
Магия цифр!

[ru0aog]

Так, пора переходить к практическим опытам.
Беру ардуиновскую плату UNO на контроллере MEGA328P,
вход А0 притянул резистором к земле, чтобы не ловить случайные наводки.
.

.
АЦП включен на опорное напряжение 1,1 В
analogReference(INTERNAL);

частота дискретизации 6 кГц,
снимаем 20 отсчётов и выводим их на СОМ-порт.
.
Подаю сигнал 1 кГц, 1Вр-р
.

.
получаем

Отсчёты
138
477
388
0
0
0
138
470
401
0
0
0
116
468
400
0
0
0
116
460
---

Хм, почему остались только полуволны?
Смотрим осциллографом, что на пине А0
.

.
а там - переход через ноль.
Всё правильно - отрицательные значения АЦП не читает.
.

0001.rar

(2.93 МБ) 809 скачиваний

[ru0aog]

Делаю смещение в половину Vref
.

.
Выставляю подстроечником половину полной шкалы АЦП
.

.
Отсчёты читаются
.

.

0002.rar

(107.94 КБ) 707 скачиваний

[ru0aog]

Продолжаем.
Будто бы всё работает.
Пора уже задействовать алгоритм Гёрцеля.
Выше уже посчитаны вручную примеры с частотой дискретизации 4 кГц, количеством отсчётов 20 и частотой настройки фильтра 1 кГц.
Так что это и проверим:
.

0003.rar

(131.76 КБ) 858 скачиваний

.
Подаём 1 кГц, размах 1,3 В р-р - судя по отсчётам, как раз входим в ограничение (значения 0 и 1023)
Магнитуда 5129.62
.

.
Смещаемся вниз, подаём на вход частоту 800 Гц, магнитуда - 5,83
.

.
Смещаемся вверх, на входе 1200 Гц - магнитуда 14,21
.

.
Неплохо работает фильтр!
Но это если не вспомнить спектрограмму, приведённую выше. 800 и 1200 Гц как раз попадают в провалы графика.
Боковые пики на 700 и 1300 Гц дают магнитуду 1500...

На случай, если вложения станут недоступными:

Код: Выделить всё

/*
Реализация алгоритма Гёрцеля на микроконтроллере ATMega328p

Основано на коде из блога Мастера Ласто
http://lasto.com/blog/index_post_1699686000.htm
Математику Ардуино хорошо объяснил Алекс Гайвер
https://alexgyver.ru/lessons/compute/
https://alexgyver.ru/arduino-math/
*/

// определим глобальные константы и переменные
const uint8_t  ON  = 1; // старт таймера
const uint8_t  OFF = 0;	// стоп таймера

const uint16_t samplingRate = 4000; // частота сэмплирования, Гц
const uint16_t samples = 20;        // количество отсчётов
      uint16_t frame[samples];      // таблица отсчётов
volatile uint16_t frameIndex;       // номер в таблице отсчётов

      uint16_t Fq = 1000;           // частота настройки фильтра
      float    k, omega, fsin, fcos, coef, mgnt; // коэффициенты Гёрцеля для частот настройки фильтра и результат обработки

#include <GyverTimers.h>            // библиотека таймера
volatile bool     timer1;           // состояние таймера


void analogReadTimer1(uint8_t onoff) {
// процедура запуска/останова таймера для чтения последовательности отсчётов с частотой сэмплирования
  switch(onoff) {
    case ON:
      frameIndex = 0;
      timer1 = true;
      Timer1.setFrequency(samplingRate);
      Timer1.enableISR(CHANNEL_B);
    break;
    case OFF:
      Timer1.stop();
      Timer1.setDefault();
    break;
  }
}

ISR(TIMER1_B) {
  // чтение отсчёта в таблицу с аналогового входа А0
    if (frameIndex < samples) frame[frameIndex] = analogRead(A0); else timer1 = false;
    frameIndex++;
}

void goertzelInit(void) {
  // определение коэффициентов для алгоритма Гёрцеля
   /* используемые внешние глобальные переменные:
    uint16_t Fq;              // частота настройки фильтра
    uint16_t samples;         // количество отсчётов
    uint16_t samplingRate;    // частота сэмплирования, Гц
   */
  // объявление локальных переменных
    float Fs = samplingRate;
    float N = samples;
  // вычисления коэффициентов для каждой из частот настройки
      k = floor(0.5 + (N * Fq) / Fs);
      omega = (2.0 * M_PI * k) / N;      // M_PI = 3.141592654
      fsin = sinf(omega);                // float sinf(float x) возвращает результат в формате float
      fcos = cosf(omega);
      coef = fcos * 2.0;
}

void goertzelLogic(void) {
  /* используемые внешние глобальные переменные:
    uint16_t samples;                   // количество отсчётов
    uint16_t frame[samples];            // таблица отсчётов
    float    fsin, fcos, coef;          // коэффициенты Гёрцеля для частот настройки фильтра
    float    mgnt;                      // результат обработки
  */
  // объявление локальных переменных
    float    real, imag;                // вещественная и мнимая составляющая результата
    float    q0, q1, q2;                // переменные алгоритма Гёрцеля

    q0 = 0.0;
    q1 = 0.0;
    q2 = 0.0;

    for (uint16_t index = 0; index < samples; index++) {
    //цикл Гёрцеля по таблице отсчётов
      q0 = coef * q1 - q2 + frame[index];
      q2 = q1;
      q1 = q0;
    }
    // вычисляем вещественную и мнимую составляющие результата
    real = q1 - q2 * fcos;
    imag = q2 * fsin;
    // вычисляем магнитуду (результат обработки)
    mgnt = sqrtf(real * real + imag * imag);
}



void GetPeak() {
  // инициализация алгоритма Гёрцеля
    goertzelInit();
  // заполнить массив отсчётов с частотой сэмплирования
    analogReadTimer1(ON);
    while (timer1) {}
    analogReadTimer1(OFF);
  // алгоритм Гёрцеля
    goertzelLogic();
  // вывести данные в СОМ-порт
    // Отсчёты
    Serial.println("Отсчёты");
      for (uint8_t m = 0; m < samples; m++) {
        Serial.println(frame[m]);
      }
    Serial.println("---");
    // Коэффициенты
    Serial.println("Коэффициенты");
        Serial.print("Fq="); Serial.println(Fq);
        Serial.print("k    ="); Serial.println(k);
        Serial.print("Omega="); Serial.println(omega);
        Serial.print("fsin ="); Serial.println(fsin);
        Serial.print("fcos ="); Serial.println(fcos);
        Serial.print("coef ="); Serial.println(coef);
        Serial.println("");
    Serial.println("---");
    // Результаты
    Serial.println("Результаты");
        Serial.print("mgnt="); Serial.println(mgnt);
    Serial.println("---");
    Serial.println("*********");
}


void setup() {
  // глобальные установки
   analogReference(INTERNAL);  // назначить режим работы АЦП с опорным напряжением 1.1 В
   pinMode(A0,INPUT);          // пин А0 - назначить входом
   Serial.begin(115200);       // инициировать СОМ-порт
}


void loop() {
  //основной цикл программы
   GetPeak();
   delay(5000);
}

[ru0aog]

Далее.
Введём полученные отсчёты в эксель и пересчитаем.

на 1000 Гц
хорошо видно постепенное накопление результата
.

.
теперь на 800 Гц
а тут длина последовательности имеет значение - подавление получилось только на последнем шаге
.

.

0003b.rar

(377.07 КБ) 825 скачиваний

.

Теперь с коэффициентами.
Их можно вычислить заранее:
k = N * Fq / Fs
omega = (2 * Pi * k) / N
sine = Sin(omega)
cosine = Cos(omega)
coeff = 2 * cosine

Определяющим является угол omega. Он зависит от соотношения частот настройки и сэмплирования
omega = (2 * Pi) * Fq / Fs
целочисленные значения косинуса (и коэффициента coeff) получаются при шаге соотношения частот через 1/4,
т.е. 0, 1/4, 1/2, 3/4, 1, 1+1/4, 1+1/2, 1+3/4, 2, 2 и 1/4 и т.д.

Выполним перебор в диапазоне звуковых частот радиотракта (от 300 до 3000 Гц) в цикле по частотам сэмплирования от 2 до 10 кГц.
Будем искать такие частоты, у которых при разности в 170 Гц получаются коэффициенты, близкие к целым числам.
При этом удвоенная наивысшая частота настройки не должна превышать частоту сэмплирования.
(макрос стартует по комбинации Ctrl+q)
.

.

Подбор.rar

(140.42 КБ) 910 скачиваний

.
Например, при частоте сэмплирования 2040 Гц, частоты 340 и 510 Гц дают коэффициенты, близкие к 1 и 0