[olegas 0]

Всем доброго времени суток!

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

 

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

 

 

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

 

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

 

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

[ViKo 1]

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

[andyp 9]

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

 

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

 

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

 

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

 

[_pv 50]

FFT требует ~100 тактов на отсчёт, соответственно не жалко потратить ещё десяток тактов на отсчёт и сделать resampling, чтобы точек было всегда 512, ну или хотя бы просто нулями буфер дозабить по длине до 512 точек.

соответственно для 30 FFT * 512 точек * 100 тактов надо 1.5-2МГц, или 10мс при 168МГц.

 

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

 

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

 

float SpecPowerGoertzel(float * y, int num, float dt, float w){
  const float K = 2.0 * cos(w * dt / num);
  float s1 = 0, s2 = 0;
  for (int i = 0; i < num; i++){
    float s0 = y[i] + K * s1 - s2;
    s2 = s1;
    s1 = s0;
  }
  return (s2*s2 + s1*s1 - K * s1 * s2) / ((float)num * num);
}

 

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

 

[olegas 0]

Дело в том, что для распознавания необходим весь спектр, поэтому алгоритм Гёрцеля не подходит.

По поводу изменения частоты дискретизации. Для сравнения результатов в MATLAB построено два спектра:

– первый (чёрный) – исходный вариант, ДПФ на 348 точки;

– второй (красный) – после передискретизации, БПФ на 512 точек.

Графики построены в координатах частота – амплитуда.

 

 

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

 

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

y = audioread('1.wav');
N = length(y);
T = 30;                %Число кадров
K = ceil(2*N/(T+1));    %Число отсчётов в кадре
Fs = 8000;            %Частота дискретизации

%Желаемые параметры сигнала
Kwant = 512;
Nwant = Kwant/2 * (T+1);    %Длина всего сигнала, чтобы число отсчётов в кадре было 512
Fswant= Nwant/N * Fs;

%Интерполируем
Y = interp1(1:N, y, 1:N/Nwant:N, 'linear');

frame1 = y(1:K);
fft_frame1 = abs(fft(frame1));

frame1Y = Y(1:Kwant);
fft_frame1Y = abs(fft(frame1Y));

fky = 1:Fs/K:Fs;
fkY = 1:Fswant/Kwant:Fswant;

fft_frame1y = zeros(size(frame1Y));
fft_frame1y(1:K) = fft_frame1;

plot(fkY, fft_frame1Y,'r', fkY,fft_frame1y,'k');
xlabel 'частота, Гц'

 

[Александр77 1]

Глядя на первый график (на котором представлены спектры для 348 и 512 Гц) создается ощущение, что частоты соотносятся друг с другом с отношением 348/512.

Если это так, то делая БПФ Вам просто необходимо умножать каждый бин на эту константу.

[_pv 50]

интерполировать надо нормально. от линейной интерполяции спектр портится особенно на высоких частотах.

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

 

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

но FFT всё равно быстрее.

если в формулу ДПФ из первого поста добавить еще K точек S = 0, так, чтобы N+K стало 512, результат от этого поменяется?

 

[olegas 0]

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

[ViKo 1]

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

[_pv 50]

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

 

[Александр77 1]

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

Взял примерные значения пиков в районе 1 кГц, 11,5 и 16,8. Их отношение 512(Гц)*11,5/16,8~=350(Гц).

Думаю у Вас более точные значения амплитуд еще больше совпадут с отношением частот дискретизации.

Учитывая имеющийся доступ к оригиналу - можете для всех точек пересчитать.

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

[x893 33]

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

[_pv 50]

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

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

а целочисленное 32х разрядое так еще и в полтора раза медленне чем float.

 

[olegas 0]

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

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

 

 

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

1.wav

Изменено 27 сентября, 2016 пользователем Olegas

[=GM= 0]

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

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

 

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