Всем доброго времени суток!
Для реализации распознавания отдельных слов потребовалось вычисление дискретного преобразования Фурье. Быстрое преобразование применить нельзя, поскольку число отсчётов не является степенью двойки (в противном случае ухудшаются результаты работы основного алгоритма). Восстановление сигнала после преобразования не требуется, всё, что нужно от ДПФ - амплитуда и частота.

Сейчас вычисляю ДПФ формулами по определению (формулы ниже, плюс ещё вычисление амплитуды). При обработке одного слова требуется выполнить 30 преобразований Фурье. Число отсчётов лежит в пределах 200 - 500 (конкретное значение заранее неизвестно), но одинаковое для всех тридцати преобразований.



В результате при реализации такого алгоритма на stm32f417 при частоте 168 МГц и арифметике с 32-битными float эти 30 преобразований (анализ одного произнесённого слова) вычисляются четыре-пять секунд, что составляет большую часть времени всех вычислений (80%).

Для вычисления коэффициентов ДПФ и квадратных корней использую функции DSP-ядра этого микроконтроллера (синус, косинус, квадратный корень). И большую часть времени (68%) занимает как раз вычисление этих самых коэффициентов.

Как ещё можно ускорить вычисление ДПФ?

Например, увеличить в 4 раза частоту выборок, и делать БПФ по 2048 точкам. И ничего не потеряете, а наоборот, я думаю, приобретете.

Для вычисления sin и cos можно использовать тот же подход, что используется в DDFS - небольшая таблица для фаз (pi/2/N, N - размер таблицы) + линейная аппроксимация между отсчетами из таблицы.

FFT требует ~100 тактов на отсчёт, соответственно не жалко потратить ещё десяток тактов на отсчёт и сделать resampling, чтобы точек было всегда 512, ну или хотя бы просто нулями буфер дозабить по длине до 512 точек.
соответственно для 30 FFT * 512 точек * 100 тактов надо 1.5-2МГц, или 10мс при 168МГц.

если же надо именно ДПФ только для определённых частот, а для не всего спектра, можно Герцелем посчитать, почти без вычисления тригонометрии (там надо только один раз косинус посчитать для каждой частоты, а не для каждого отсчёта).

будет несколько тысяч тактов на отсчёт, но если надо только несколько частот, а не весь спектр, может быть быстрее чем через FFT



только корень еще навеное взять надо и почему оно на длину выборки делилось, причём дважды, уже не вспомню.

Дело в том, что для распознавания необходим весь спектр, поэтому алгоритм Гёрцеля не подходит.
По поводу изменения частоты дискретизации. Для сравнения результатов в MATLAB построено два спектра:
– первый (чёрный) – исходный вариант, ДПФ на 348 точки;
– второй (красный) – после передискретизации, БПФ на 512 точек.
Графики построены в координатах частота – амплитуда.




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

Вот пример кода из MATLAB:

Глядя на первый график (на котором представлены спектры для 348 и 512 Гц) создается ощущение, что частоты соотносятся друг с другом с отношением 348/512.
Если это так, то делая БПФ Вам просто необходимо умножать каждый бин на эту константу.

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

Герцель это и есть ДПФ только без вычисления тригонометрии, соответственно вычисляется заметно быстрее простого ДПФ с синусами/косинусами, даже если посчитать его для всех частот.
но FFT всё равно быстрее.
если в формулу ДПФ из первого поста добавить еще K точек S[i] = 0, так, чтобы N+K стало 512, результат от этого поменяется?

То, что частоты сдвинуты - не страшно. Гораздо хуже, что амплитуды для двух вариантов для частот 0..4000 Гц не совпадают.

Так вы не интерполируйте свой файл, тем более, линейно. Считайте, что он взят с другой частотой, и ищите, соответственно другие.

этот файл "1.wav" сюда прикрепите

Взял примерные значения пиков в районе 1 кГц, 11,5 и 16,8. Их отношение 512(Гц)*11,5/16,8~=350(Гц).
Думаю у Вас более точные значения амплитуд еще больше совпадут с отношением частот дискретизации.
Учитывая имеющийся доступ к оригиналу - можете для всех точек пересчитать.
Еще лучше - возьмите максимальные значения в каждой реализации и относительно него пронормируйте полученные данные.

Когда то для ускорения использовали только uint16 - примерно 1-2 mS на lsi-11

на M4 не особо поможет, у него FPU есть, для быстрого преобразования Фурье разница в тактах между 16ти разрядным целым и 32х разрядным float меньше 20%.
а целочисленное 32х разрядое так еще и в полтора раза медленне чем float.

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


И вот подопытный файл для всех желающих.

Ну так, заранее вычислите эти коэффициенты и разместите их в таблице.

У меня на 100 МГц ДСП двойной МАК для преобразования Фурье вычисляется за 10 нс. Т.о., для 200-500 точек время счёта Re/Im одной палки составит порядка 2-5 мкс. Для 30 частот время счета составит 60-150 мкс. Но это прикид для 100 МГц тактовой, а у вас 168 МГц тактовая, так что о каких секундах может идти речь?

если буфер добить нулями и сделать FFT то при правильной нормировке (Sqrt[512/348]) разницы нет и быть собственно не должно.
небольшая разница есть только из-за того что точек всё-таки не сильно много, и сетки частот разные.
Нажмите для просмотра прикрепленного файла
зелёный график - ДФТ посчитанное Герцелем с шагом 1Гц. синий(348 точек) и желтый (512 точек = 348 + нули) график просто берут из него точки по разным сеткам частоты 8000/348 = 23Гц и 8000/512 = 15.6Гц соответственно

но если ваши алгоритмы распознавания чувствительны к такой разнице, то думаю можно сразу закапывать.

Я, наверное, неоптимально представляю, но для N отсчётов ДПФ требуется матрица размером N^2. Взяв по максимуму N=500 получим: 4*N^2=4*500^2/1024=976 КБ. Многовато.

1) считайте ДПФ Герцелем - пример в предыдущем посте на картинке. там только один косинус для каждой частоты считать надо. и можно и без таблицы.
2) делайте FFT.

Ну прямо в лоб-то считать не надо. Ваша таблица по рядам (и по столбцам) будет многократно повторяться.

Лучше составьте таблицу одного периода синуса на N=500 точек (для двоичного ЦПУ практичнее брать число N степени двойки, например N=512).
Для первой палки считаете МАКи, беря синусы: 0,1,2,..,499 (всего 500), для второй палки берете каждый второй: 0,2,4,..,498,0,2,4,..,498 (всего 500).
Для третьей палки берете каждый третий:0,3,6,.. (всего 500) и так далее.

Ну и не надо считать компоненты спектра (пресловутые "палки") выше частоты Найквиста, поскольку амплитудный спектр симметричен.

Для обработки речи что ДПФ, что БПФ - всё едино, т.к. любой микропроцессор успевает.
Например, в телефонном вокодере по 200 выборкам делают автокорреляцию,
по автокорреляции строят АР-фильтр, через фильтр пропускают синусоиды разной частоты,
выбирают синусоиды, на которые фильтр хорошо окликается
и делают интерполяцию частоты речевого сигнала.
На все это и многое другое типа вычисление периода основного тона, подбор сигнала возбуждения
тратится всего 20 млн.оп./с.