Алгоритм поиска основной частоты, Метод максимального правдоподобия Опции

[QuickNick]

Здравствуйте, товарищи.

Мне нужно улучшить алгоритм поиска основной частоты по методу максимального правдоподобия. Был пример у меня, как делать, я его модифицировал в нескольких местах (добавил настраиваемость оконной функции и фильтра). Возможно, некоторые места вам покажутся совершенно абсурдными - но я пока новичок в DSP, некоторые места в старом алгоритме вызывают у меня сильное удивление, но пока не хватает опыта для однозначного решения: или я не понимаю момент, или предыдущий программист допустил косяк.

- Вход: массив байтов byteData длиной sampleLength, количество разбиений sampleRate.
- Выход: фундаментальная частота fundamentalFrequency и громкость fundamentalVolume.

Алгоритм:
1. Массив байтов преобразуется в вешественный массив амплитуд.
2. Проверяется максимальная амплитуда - если она ниже определённого значения, то это квалифицируется как почти полная тишина (выход из алгоритма)
3. Получение спектра из массива амплитуд:
3.1. Умножение массива амплитуд на оконную функцию;
3.2. Преобразование Фурье.
4. Фильтрация массива спектральных значений.
5. Расчёт массива громкости (модуль от спектральных значений).
6. Нормализация громкости:
6.1. Поиск максимальной громкости (maxVolume).
6.2. Нормализация (volume[i] = volume[i]/maxVolume);
7. Квалификация спектра как информативного или зашумлённого.
8. Инициализация 3 пиков громкости и соответствующих им частот.
9. Поиск фундаментальной частоты.
9.1. Если один из 2 главных пиков громкости ниже некоторого порога громкости, то выбрать частоту другого главного пика, выход из алгоритма.
9.2. Поиск фундаментальной частоты по 2 главным пикам согласно известным пропорциям 1/2, 1/3, 2/3, 2/5, 3/4, 3/5.
9.3. Коррекция фундаментальной частоты по 3-му пику.

Мне нужно пробежаться по некоторым пунктам, чтобы удостовериться в валидности алгоритма.
Пункты 1, 2, 6, 8 и 9 вопросов не вызывают.

Самые горячие вопросы сейчас такие:
1. Правильно ли считается громкость? Массив спектральных значений имеет следующую структуру: в ячейках spectre[2*i] лежат коэффициэнты при косинусах, в spectre[2*i+1] лежат коэффициэнты при синусах. Соответственно, громкость считаю так: volume[i] = sqrt(sqr(spectre[2*i])+sqr(spectre[2*i+1])); Может, таким методом не громкость считается, а что-то другое?
2. Как рассчитать шум? И как рассчитать порог шума (если его превысить, то выход из алгоритма)? Сейчас шум рассчитывается как среднее арифметическое по массиву громкости.

Сообщение отредактировал QuickNick - Jul 1 2011, 07:00

[thermit]

QuickNick:
Вот и нарвались на ступеньку sad.gif
Так что теперь нужны более точные методы.

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

Явление гиббса еще никто не отменял. Методы борьбы с ним широко известны. Возьмите окно кайзера с параметром ~30. Будет Вам равенство значений соседних точек, если частота попадает точно в середину между бинами.

[QuickNick]

Явление гиббса еще никто не отменял. Методы борьбы с ним широко известны. Возьмите окно кайзера с параметром ~30. Будет Вам равенство значений соседних точек, если частота попадает точно в середину между бинами.

Уже расправился - пока использую интерполяцию Ньютона по 5 точкам, затем с шагом по 0.1 Гц прохожу по интервалу (он получается длиной примерно 3.9 Гц), нахожу максимум.
Возможно, есть какой-то хитрый способ интерполяции наоборот, когда отображение не инъективно, но за полдня найти его не успел.