Гёрцель - как убрать разброс (болтанку) в нч ?

Мечта идиота сделать диммер для мотора с обратной связью по звуку вращения. Принцип задуман простой - вручную устанавливаешь диммером нужную скорость, нажимаешь кнопку - он считывает "главный звук" и старается его держать неизменным. Видно что для 1500 Гц входных - болтанка +- 3 Герца. Это в допуске. Для 100 Гц входных - болтанка уже +-15 Гц. Чем ниже частота - болтанка еще возрастёт. Число оцифровок увеличить нельзя - мало памяти.
Частоты заданы как float NNN.xxxx - c 4-мя знаками после запятой.
Под болтанкой - речь идет о специфике алгоритма Гёрцеля, а не механики.

#define SAMPLE_RATE = 8800 ;
#define SampleNum = 100 ;

F_CPU = 16500000 ;

const float Freq_Min = 10.0000;
const float Freq_Max = 2000.0000;
const float Freq_Step = 5.0000;

Вопрос - как победить. С уважением.

 

Не представляю что ЭТО, но может ЭТО умножить на 1000 или 10000 и работать с целыми числами, если возможно. Как направление мысли.

 

А зачем здесь Герцель, если нужен ПИД?
И таки да: если частоты отличаются на порядок, то фильтры, построенные на идеях дискретного преобразования Фурье, требуют децимации.

 

Понятия не имею, как автор считает гармоники - какие, сколько, но для сглаживания биений боковых лепестков применяют окна. На мой взгляд самое лучшее - Blackman. Но формула там адовая). Как обратный эффект от применения окна мы получим расширение главного лепестка.

 

Ойё. Не надеялся что вообще кто-то ответит, судя по "пустоте в зале". Здравствуйте коллеги. Разрешите поздравить с матушкой - зимой !

 

Шановне)

Я так понимаю, что алгоритмом Герцеля делается попытка вычислить гармонику с максимальной амплитудой, а через неё скорость вращения?
Смотрим:

Исходя из этих вводных, фундаментальная гармоника расположена на частоте 88 Гц. Первая ближайшая целая - вторая - на частоте 176 Гц, третья - 264. В окрестности частоты 1,5 кГц это будут 17 и 18 гармоники. При этом дискретность спектра составит всё те же 88 Гц. Т.е. если решать задачу в лоб, то нам надо одновременно считать сразу 20 гармоник. Очевидно, что дискретность в 88 Гц в области частот выше килогерца дадут гораздо меньший разброс, чем в области 100 Гц. Решение тут может быть только одно: алгоритм Герцеля позволяет вычислять не только целые гармоники, но и дробные. Если не децимироваться, то необходимо как-то грубо прикинуть тот диапазон частот, в котором мы можем находиться и внутри него вести расчёт. Если нет, то брать ПЛИС и синтезировать внутри неё схему, которая будет считать сразу хоть сто гармоник (ну, насколько хватит ячеек). Если сто гармоник, крайняя частота 1,5 кГц, то первая гармоника должна быть равной 15 Гц. И дискретность всей системы будет также 15 Гц.

 

Сановнэ. Да именно так. Но - Нам нет необходимости знать частоту вращения в цифрах, нам нужно ее просто держать постоянной. Не зависимо от того орет ли сам двигатель на частоте вращения или трещат зубья фрезы на валу. Обычно сам движок особо не шумит.

Увеличил число отсчетов до 150-ти. Болтанка съузилась до приемлемой. ОЗУ еще позволяет, но будем искать.

 

Всё твои частоты так или иначе должны совпадать с частотами вычисляемых гармоник. Если количество отсчётов увеличен до 150, то первая фиксированная частота будет 58,5, вторая 117, третья 176 и так далее. Поддерживать частоты между ними ты не сможешь - просто по ним нет данных.

 

Честно говоря, не понимаю о чем вы. В цикле провожу анализ с шагом 5 Герц. Конечно это много циклов. Но оно быстро. Все равно бугор проявится и прога приклеит его к ближайшей 5ти-герцно-кратной частоте. И всё. Запомнит частоту и будет регулировать так чтоб пик был на этой запомненной частоте. Двигателю точность "на слух" лучше 5-герц не нужна. Ничто не запрещает сужать границы поиска и уменьшать шаг обсчета.

Спойлер

//==================================
void loop() {

TONE_MAX = 0.0000;
FREQ_MAX = 0.0000;

char ValStr[17];

// Производим оцифровку звука в буфер
for (int n = 0; n < numSampl; n++) { // Срисовываем волну звука
Samples[n] = analogRead(inputADC); // Get samples from ADC
}

// Производим преобразование Герцеля для очереди фиксированных частот
for (float Freq = Freq_Min; Freq <= Freq_Max; Freq += Freq_Step) {

Goertzel TONE(Freq, SAMPLE_RATE); //Запускаем преобразование

TONE_MAG = TONE.Mag(Samples, numSampl); // Получаем магнитуду частоты

if (TONE_MAG > TONE_MAX) {
TONE_MAX = TONE_MAG;
FREQ_MAX = (Freq * 1.0670); // Специфика цифрового Герцеля/Фурье
}
// delay(10);
}

prints_P(PSTR(" FREQ=")); // Печать строковой константы
itoa(FREQ_MAX, ValStr, 10); // Конвертация числа в строку
prints(ValStr); // Печать строковой переменной
prints_P(PSTR(" MAG=")); // Печать строковой константы
itoa(TONE_MAX, ValStr, 10); // Конвертация числа в строку
prints(ValStr); // Печать строковой переменно
prints_P(PSTR("\r\n")); // Перевод строки

// delay(500);
// prints_P(PSTR("Next\r\n\r\n")); // Печать строковой константы
}

Тут фишка именно в алгоритме.

 

5 Герц откуда взялись? По твоим входным данным никакими 5 ти герцами не пахнет.

У алгоритма Герцеля ровно две фишки:
1. Можно считать дробные гармоники, но всегда больше единицы. Субгармоники он не считает.
2. Коэффициенты косинуса и синуса для всех вычисляемых гармоник можно посчитать заранее. В большинстве применений это даёт прирост в вычислениях по сравнению с Бпф.
Есть и минус - фильтры на основе алгоритма Герцеля не устойчивы, поэтому общего решения их создания на последовательности отсчётов нет. Каждый новый расчёт всегда начинается с 2 мя обнуленными последними отсчётами.
Я пмсал алгоритм Герцеля на языке Julia - для набора чистых синусов, которые чётко встают на рассчитаные гармоники, никакой болтанки нет. Если частоты попадают меж гармоник, то их вычисление невозможно. ДПФ, ГЕРЦЕЛЬ - это не вычисление частоты, это вычисление гармоник. Фундаментальная гармоника (она же первая) считается заранее по частоте дискретизации и количеству отсчётов на период. По гармоникам можно оценить спектр, а можно восстановить исходный сигнал.

НАПРИМЕР
В сети 50 Гц. Мы берём стандартные 40 отсчётов на период и типовую частоту дискретизации 2 кГц. Первая гармоника встаёт ровно на 50 Гц. По алгоритму Герцеля мы посчитали первые типовые 7 гармоник и получили представление о качестве сети. Но вот первая гармоника просела по амплитуде. Что это? Частота уплыла или просело напряжение? Без дополнительных измерений мы это сказать не сможем. Ни по Герцелю, ни по ДПФ/БПФ.

 

В самом начале показал. Ну да, надо было код полностью выложить.
const float Freq_Min = 10.0000;
const float Freq_Max = 2000.0000;
const float Freq_Step = 5.0000;

 

Что подразумевать под дополнительными измерениями ? В данном случае измерение это - оцифровка. Имея "образцовый" массив как константу - и один "текущий" массив мы можем сделать все выводы какие бывают. Например видите что просела магнитуда главного пика. Два варианта - просело напряжение или ушла частота. Делаете два симметричных обсчета ниже и выще по частоте. Если они равны - ухода частоты нет. Делаете обсчет боковых - если магнитуды просели пропорционально главной - просело напряжение. Если изменение не пропорционально - искажение синусоиды. Конечно у идеальной синусоиды нет боковых совсем. Но это у идеальной.
Наверно так, согласно здравому размышлизму.
.

 

Нет, не сделаешь. Имея ацп, который я описал выше (40 отсчётов и 2кГц), ниже 50 Гц ничего не посчитаешь. Ну или посчитаешь, но с ошибками. Посчитать частоту 51 Гц можно, но для этого надо задать гармонику с номером 1.02 (50 * 1,02 = 51). Но если поплыла основная частота, то поплыли и гармоники. И их номера тоже надо смещать. Это все сложно и поэтому частоту электросети считают по-другому.

Ширина главного лепестка и наличие боковых лепестков в спектре определяется не искажениями синусоиды, а корреляционной функцией синуса и прямоугольного сигнала. Синус откуда берётся понятно, а прямоугольный сигнал - это набор наших отсчётов с неизменной амплитудой. Как убрать боковые лепестки я писал в самом начале. В любой книжке по цос все это описано.

Ещё раз - дпф и все его производные не оперируют частотами. Только гармоники. У гармоники есть свойства: её номер, амплитуда и фаза. Сопоставление номера гармоники к частоте происходит снаружи, вне всех этих алгоритмов (для этого надо знать период (или частоту) дискретизации и количество отсчётов на один период фундаментальной гармоники. Разумеется, мы это знаем, т. к. это параметры нашего ацп).

 

Оцифровка - это переход от аналогового сигнала к дискретному с неизбежной потерей информации. К измерениям никакого отношения не имеет.

 

Что такое образцовый массив? Можно перечислить все эти выводы и способы их получения.

 

Хорошо. Оставим на время тему. Ответьте или ткните носом на что-нибудь для чайников. Вопрос - оцифровка входной частоты 1КГц.
1. Почему максимум магнитуды получается на 911 Гц ?
2. Почему на плавно-спадающих "крыльях" наблюдаются провалы на (примерно) 911- 60 и на 911+60 Гц. Можно начинать цикл с любой частоты и на любой заканчивать - картинка не меняется. С уважением.

Спойлер

#define SAMPLE_RATE = 8800 ;
#define SampleNum = 150 ;

F_CPU = 16500000 ;

const float Freq_Min = 810.0000;
const float Freq_Max = 1010.0000;
const float Freq_Step = 1.0000;

ПС. Где-то видел объяснение, но увы, не могу найти заново.
ППС. Нашёл эмпирически - расстояние от пика до провала (и между соседними провалами) равно -
SAMPLE_RATE / SampleNum Гц.

 

У меня под рукой только Julia, но на си перекладывается очень просто. Но даже можно ничего не перекладывать - алгоритм реализован на всех языках мира)
Задаём частоту дискретизации:

Код (Text):

f = 51200   # оптимизация под выборку в 512 отсчётов

У меня целочисленная арифметика, поэтому с оптимизациями. Целочисленная, так как планирую загнать её в ПЛИС.
Находим период дискретизации

Код (Text):

dt = 1 / f

Считаем первую гармонику:

Код (Text):

N = length(buf) # длина буфера 512 отсчётов
freq1 = 1 / (dt * N) #частота первой гармоники

Теперь сам Герцель

Код (Text):

function goertzel(buf, k)
    N = length(buf)
    if N == 0 || k == 0.0
        return
    end
    u0 = 0
    u1 = 0
    u2 = 0
    ω = 2 * pi * k / N
    α = 2 * cos(ω)

    for i = 1:N
        u0 = α * u1 - u2 + buf[i]
        u2 = u1
        u1 = u0
    end

    return α * u1 / 2 - u2, sin(ω) * u1
end

Здесь buf - это вектор. На языке Си - это массив с последовательными отсчётами. Количество отсчётов, как мы договорились ранее - 512. k - номер (номер!) вычисляемой гармоники. Соответственно, если мы хотим найти 10 каких-то гармоник, нам эту функцию надо вызвать 10 раз. 100 гармоник - 100 раз с разными к. Здесь мы сразу видим, что поворачивающие коэффициенты 2 * cos(ω) и sin(ω) мы можем посчитать заранее для всего множества к и просто записать в память. Герцель нам выдаёт результат в комплексном виде. Запишем его в удобном нам виде:

Код (Text):

Re, Im = goertzel(buf, curK)

curK - номер рассчитываемой гармоники
Считаем амплитуду (фазу не считаем - она нам не нужна):

Код (Text):

a = ((Re * Re + Im * Im)^(1 / 2)) * 2 / N

Таким образом, для всех к мы нашли амплитуды a. Теперь перебираем все наши a и находим максимальное значение. Сопоставляем a номеру к и находим частоту:

Код (Text):

freq = freq1 * к

Вот как бы и всё.
Теперь картинки.
Вот я подаю сигнал. Это косинус + 3-я гармоника + этот же сигнал, но сдвинутый на pi/2 и с небольшой амплитудой + шум.

Видно, что здесь ровно 512 отсчётов. Сигнал очень кривой.
Посчитал 650 гармоник. Построил график (на графике гармоники уже приведены к частоте)

За кадром нашёл пиковую гармонику и её частоту

Код (Text):

Частота точно: 354.0

Заданная частота:

Код (Text):

f1 = 351

Разбег всего 3 Герца, что для такого сигнала и такой выборки считаю отличным результатом.

Постулат: гармоника с наибольшей амплитудой не находится сама. Нам необходимо внешним алгоритмом сделать предположение, где она может находится и самостоятельно перебирать все номера в окрестности.

 

что-то я это пропустил. Тогда нам вообще не надо считать частоту. И даже думать в её сторону. Оперировать только номерами гармоник. Это не сильно облегчит вычисления, но по крупице, по зернышку...