Войти
ФлеймФорумОбщее

На замену FFT!

#0
14:06, 17 июня 2010

Нашёл интересное обсуждение: http://www.dsprelated.com/showmessage/7716/1.php К автору (Dmitry Terez) отнеслись с недоверием, хотя он предоставил исходник для Mathcad.

Я попробовал реализовать его метод на Delphi и оно заработало, весьма даже неплохо. Значения лишь немного отличаются от аналогичных FFT, хотя есть вопросы по "настройке метода". С шумом метод тоже неплохо справляется. Размер окна можно применять любой. Есть и недостатки - по моим догадкам можно выловить одну-две самые мощные частоты, полный спектр как с FFT не получится. Но метод несомненно имеет и свои преимущества, в том числе и применительно к DSP. Особенно когда нужно выловить один пик частоты.

Как известно, на DSP любят считать FFT "бабочками", комбинируя их разными способами, чтобы оптимизировать число умножений. Но умножени всё равно остаётся дофига. А если в DSP без плавающей точки, то там продолжают шаманить с расширением битности и эмулировать эту плавающую точку. Ну и как вы понимаете, пытаясь увеличить точность расчёта, берут FFT с окном по-больше... 1024, 2048, 4096 и т.д., при этом внушительно растёт число умножений, и растёт ошибка округления. Можно представить как хреново будет выглядеть посчитанный на 8 битах integer FFT с окном 4096. По этому и по другим причинам на 8-битах считают FFT возможно только с малюсеньким окном.

Что же за метод предлагает Dmitry Terez?
А метод прост как пробка и интуитивно понятен.

1. Считываем кусочек из WAV файла в массив wave[j]

2. Нормализуем (приводим значения в массиве wave[j] к диапазону -1..+1)

3. Перебрасываем все данные из массива wave[j] в массив vector[k] слудеющим образом:
vector[k].x=wave[j]  -  каждое значение входного массива рассматривается как координата вектора X в одномерном пространстве. Также м
vector[k].time=j - запоминаем "время" этого вектора относительно начала файла, измеренное в сэмплах

4. Задаём некое расстояние равное к примеру RADIUS=0.15

5. Для всех пар  vector[k] vector[m], для которых abs(vector[k].x-vector[m].x)<RADIUS
увеличиваем на единицу значение в финальном масиве hist[n]:=hist[n]+1,
где n = abs(vector[k].time-vector[m].time) - период между сэмплами
Проще говоря, мы смотрим сэмплы с примерно одинаковыми значениями по оси Y и увеличиваем столбик в гистограмме, соответствующий периоду между сэмплами.

6. Находим максимумы n_max в массиве hist[n]. Эти максимумы соответствуют периоду сигнала
Например для WAV файла с частотой 44100kHz частоту получаем по формуле
freq=44100/n_max.

Удивительно простой метод определения частоты. Причём возможно реализовать вовсе без умножений, если создать lookup таблицу значений "44100/n".

Пример-сравнение анализа реального сигнала в AUDACITY с помощью FFT и методом Terez-а.
Dmitry Terez Frequency Estimation Method | На замену FFT!
FFT с окном 8192 даёт пик  806 Герц
Метод Terez-а с окном 512(!!!) даёт 801 Герц


Набросок алгоритма:

 AssignFile(F, '0.wav');
  Reset(F);
  BlockRead(F, Header, SizeOf(Header));
  BlockRead(F, Buffer, SizeOf(Buffer), Readed);
   for i:=1 to N do
     signal[i-1]:= (Buffer[i]);
//    signal[i-1]:=sin(i*(2*Pi*2000)/44100) {+sin(i*(2*Pi*150)/44100)};

  //  +2.1*sin(i*(2*Pi*8944)/44100);

  CloseFile(F);


  min:=0;
  max:=0;
  for i := 0 to N - 1 do
   begin
    if signal[i]>max then max:=signal[i];
    if signal[i]<min then min:=signal[i];
   end;

  for i := 0 to N - 1 do
   begin
     signal[i]:=(signal[i]-min)/(max-min)-0.5;
   end;



   tau:=scrollbar3.Position;

  setlength(vectors,N);
  for i := 0 to N - 1 do
   begin
     Vectors[i].x:=signal[i];
{     Vectors[i].y:=signal[i+tau];
     Vectors[i].z:=signal[i+tau*2];
     Vectors[i].w:=signal[i+tau*3];}


{
     Vectors[i].y:=0;
     Vectors[i].z:=0;
     Vectors[i].w:=0;
 }

     Vectors[i].time:=i;
     timez1[i]:=0;
     timez2[i]:=0;     
   end;


   for I := 0 to N - 1 do
    for J := I to N - 1 do
    if i<>j then

     if sqrt((Vectors[i].x-Vectors[j].x)*(Vectors[i].x-Vectors[j].x)+
     (Vectors[i].y-Vectors[j].y)*(Vectors[i].y-Vectors[j].y)+
     (Vectors[i].z-Vectors[j].z)*(Vectors[i].z-Vectors[j].z)+
     (Vectors[i].w-Vectors[j].w)*(Vectors[i].w-Vectors[j].w))
     <radius then
//     if abs(vectors[i].x-vectors[j].x)<radius then

      begin
       inc(timez1[abs(Vectors[i].time-Vectors[j].time)]);

       if timez1[abs(Vectors[i].time-Vectors[j].time)]>maxbig then
       maxbig:=timez1[abs(Vectors[i].time-Vectors[j].time)];



      end;

      for i := 1 to N - 2 do
       timez2[i]:=round(
                  timez1[i-1]*0.25+
                  timez1[i]*0.5+
                  timez1[i+1]*0.25);



  for i := 0 to N - 1 do
   begin
    k1:=i-20;
    k2:=i+20;
    if k1<0 then k1:=0;
    if k2>N-1 then k2:=N-1;

    k3:=0;
    for j := k1 to k2 do
     if timez2[i]<timez2[j] then k3:=1;

//     if k3=1 then timez1[i]:=0 else timez1[i]:=timez2[i];
                                      timez1[i]:=timez2[i]




   end;

  for i := 1 to N - 2 do
       timez2[i]:=timez1[i];

   max:=0;
      for i := 1 to N - 2 do
     if timez2[i]>max then
        if i>15 then
        begin
         max:=timez2[i];
         max_i:=i;
        end;






             {
        k1:=max_i-5;
        k2:=max_i+5;
        if k1<0 then k1:=0;
        if k1>N-1 then k1:=N;


        sum1:=0;
        sum2:=0;
        for I := k1 to k2 do
         begin
          inc(sum1,i*timez2[i]);
          inc(sum2,timez2[i]);
         end;

         max_i:=round(sum1/sum2);
              }




image3.Picture.Bitmap.Width:=N;
image3.Picture.Bitmap.Height:=round(maxbig)+10;
image3.Picture.Bitmap.Canvas.Brush.Color:=clBlack;
        image3.Picture.Bitmap.Canvas.FillRect(image3.Picture.Bitmap.Canvas.ClipRect);

image3.Refresh;
qp3.Attach(image3.Picture.Bitmap);


    for I := 0 to N - 1 do
     begin
     for j := 0 to timez2[i] do
     qp3.SetPixel(i, j ,rgb(255,0,0));


     for j := timez2[i] to round(max) do
     qp3.SetPixel(i, j ,rgb(0,0,0));

     if i = max_i then
     for j := 0 to timez2[i]-100 do
     qp3.SetPixel(i, j ,rgb(255,255,0));

     if i = ifreq then
     for j := 0 to round(max) do
     qp3.SetPixel(i, j ,rgb(255,0,255));

     if i = ifreq div 2  then
     for j := 0 to round(max) do
     qp3.SetPixel(i, j ,rgb(255,0,255));


     end;

    image3.Refresh;



  exit;

Как видно в моей программе есть небольшие отличия: а именно вектор имеет аж четыре координаты x,y,z,w. Terez советует их использовать, "запихивая" таким образом сигнал в пространство N-мерных векторов (в данном случае четырёхмерных). Выгоды особой я от этого не заметил, кроме как уменьшения амплитуды несуществующих частот.

также в одномерном случае конструкция

   if sqrt((Vectors[i].x-Vectors[j].x)*(Vectors[i].x-Vectors[j].x)+
     (Vectors[i].y-Vectors[j].y)*(Vectors[i].y-Vectors[j].y)+
     (Vectors[i].z-Vectors[j].z)*(Vectors[i].z-Vectors[j].z)+
     (Vectors[i].w-Vectors[j].w)*(Vectors[i].w-Vectors[j].w))
     <radius then
превращается в abs(vector[k].x-vector[m].x)<RADIUS.


Но по методу остаются некоторые непонятности:
при выборе радиуса RADIUS частота определяется стабильно при определённом значении скачок на 50Гц. Хотелось бы в этом разобраться и выжать из алгоритма максимум.

Кого заинтересовал метод, прошу подключиться к исследованию нового направления.
p.s. Есть идея добавить локальный вевлет в окрестности найденной частоты для получения высокого разрешения по времени. CWT,DWT,FFT,TerezMethod всё проработано, без раздутой теории =)


#1
14:16, 17 июня 2010

tmtlib
Правильно я понимаю, что этот метод может только находить гармонику с максимальной амплитудой? Это действительно настолько важная и распространенная задача?

И да, кто-нибудь смог обнаружить его пейпер по приведенной ссылке? Я потыкался по сайту, что-то не понял, куда смотреть, а теоретическое основание хотелось бы заценить.

#2
14:24, 17 июня 2010

Manticore

> Это действительно настолько важная и распространенная задача?
Определение частоты основного тона (если о ней идет речь) это вообще центральная задача в анализе/кодировании речи.

#3
14:29, 17 июня 2010

Ghost2
> Определение частоты основного тона (если о ней идет речь) это вообще центральная задача в анализе/кодировании речи.

А можно чуть поподробнее, если не трудно?

#4
15:33, 17 июня 2010

Manticore

> А можно чуть поподробнее, если не трудно?
Не, трудно ) У нас на кафедре защищается один по этой теме, так вот там только про это и говорят.

#5
16:03, 17 июня 2010

Manticore
> Правильно я понимаю, что этот метод может только находить гармонику с
> максимальной амплитудой?

По идее и двух, возможно более.
Возможно имеют место следующие проблемы:
От частоты скажем 800Гц ещё возникают пики с меньшей амплитудой 400, 200, 100, 50 и т.д., так как периоду T соответствуют синусоиды 2*T, 4*T,... Можно предположить что они будут конфликтовать с реальными частотами 400,200 и т.п..
Вторая вероятная проблема заключается в том, что чем больше в сигнале различных частот, тем меньше будут по высоте их пики. Так как длина сигнала в семплах одинаковая, то в гистограмме сумма столбцов = константе (числу самплов). Встаёт вопрос какая частота съест больше.

> Это действительно настолько важная и распространенная
> задача?
Задача в принципе достаточно распространённая. Помимо упомянутого распознания речи:
- уровнеметры, и измерители расстояний, основанные на ЛЧМ (линейно-частотной модуляции).
- измерение доплеровского смещения (скорость объектов)
- во многих сигналах нам достаточно знать главный пик
Не совсем ясна применимость к примеру в фурье-спектрометрии. Там нужен целый спектр.
А вот как основная частота, в окрестности которой сделать непрерывное вейвлет преобразование - может сгодиться.

> И да, кто-нибудь смог обнаружить его пейпер по приведенной ссылке? Я потыкался
> по сайту, что-то не понял, куда смотреть, а теоретическое основание хотелось бы
> заценить.
Как видно из обсуждения на www.dsprelated.com пейпер ещё тот (фиг разберёшь). Сайт автора к сожалению уже не существует (домен занят по сути спамерской фигнёй по ключевым словам FFT и т.п.). Оригинальный пейпер не нашёл, разве что по имени Terez что-то на эту тему.
Можешь поискать его патент в гугле, там есть номер патента, но он сильно углублён в сторону распознания речи, чем про сам метод.

В принципе метод достаточно прозрачен:
В некотором выдуманном пространстве, например трёхмерном, рисуем точки с координатами (x,y,z).
Координаты x,y,z берём из WAV файла в цикле:

x:=файл.wav[i];
y:=файл.wav[i+tau];
z:=файл.wav[i+tau*2];
где tau берём от балды.
Понятно, что в этом трёхмерном вымышленном кубе у нас образуются некоторые яркие скопления точек (векторов).
Также не забываем при добавлении яркости точки запомнить её положение в wav-файле (индекс i).
Затем выбираем "допуск" - некий радиус, и перебираем все пары точек, и находим только те точки, расстояние между которыми меньше радиуса RADIUS. При нахождении пары таких точек - берём по модулю разность их положения в wav-файле, что соответствует периоду. Найденные периоды заносятся в гистограмму (растут вертикальные столбики, соотствующие найденным периодам). По простейшей формуле, данной выше постом, находим из периода частоту.

Если взять отдельно синусоиду - у неё есть места, которые часто повторяются. Особенно часто будут повторятьтся "верхушки", так как Y в окрестности верхушки меняется медленно, так как мы помним их координаты в wav файле - найдутся и периоды. В сложных запутанных случаях этот метод может выявить то, на чём FFT пролетит из-за большого окна и фазы, а wavelet будет долго и упорно делать свой continious. Эти случаи частные, но интересные.

#6
5:20, 18 июня 2010

update:
подробное описание алгоритма от Terez-а (pdf с картинками можно скачать, сначала жмите download pdf, а потом на дискетку):
http://www.freepatentsonline.com/y2003/0088401.html
http://www.freepatentsonline.com/7124075.html
Описание не такое уж и сложное, оппоненты на www.dsprelated.com оказались слабоваты, видимо не читывали геймдевелоперских пейперов по графике.
Из беглого просмотра видно, что многомерность его идеи я несколько недооценил, так что исходник можно существенно улучшить. Результаты выложу на igrodel.ru.

#7
15:32, 18 июня 2010

update:
terez method 2 | На замену FFT!
Тесты:
FFT1024 сэмплов - главный пик 2466Гц
FFT128 сэмплов - главный пик 2488Гц

Terez 256 сэмплов - главный пик 2450Гц
Terez 128 сэмплов - главный пик 2450Гц

Улучшено распознавание пиков - берётся первый пик слева, независимо от амплитуды:

min:=hist[0];
max_index:=0;
thr:=strtoint(Edit10.text);
for i := 0 to len - 1 do
 begin
  if hist[i]<min then min:=hist[i];

  if hist[i]-min>thr then
   begin
    max:=hist[i];
    max_index:=i;
    break
   end;

 end;

    for j := max_index to max_index+5 do
     if hist[j]>max then
      begin
        max_index:=j;
        max:=hist[j];
      end;

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

#8
15:02, 20 июня 2010

update: повышена точность определения частоты. Для этого находится центр масс пика.
Для зашумленных сигналов и сигналов с глючной DC введен простой и быстрый ВЧ фильтр.

Провёл тесты на различных сигналах с умеренным шумом, и пришёл к следующим выводам:

- метод даёт точность сопоставимую с FFT и большим окном
- в отличие от FFT даёт высокое временное разрешение по амплитуде, практически аналог CWT
- не имеет присущее FFT реакцию с запозданием на измение параметров сигнала
- во много раз превосходит по скорости как FFT и CWT
- метод чувствителен к глюкам DC, нужен ВЧ фильтр
- ошибочное определение частоты можно компенсировать огромным числом измерений и построением гистограммы наиболее вероятной частоты

Из графиков сравнения двух методов: FFT-4096 и TEREZ-256 видно, что в основном кривые совпадают, а детальный анализ мест неправильного определения частоты в основном соответствует спорным кускам сигнала. Ошибки в FFT и TEREZ часто совпадают, но есть случаи когда ошибается только FFT или только TEREZ.
Dmitry Terez method | На замену FFT!
красным цветом показано FFT-4096, белым - TEREZ-256.
Новый метод выглядит даже несколько стабильнее, чем FFT. Но здесь нужно учитывать отсутствие сглаживающего окна у FFT и его простейшей реализации.

Ну и самое интересное -
Практически вся суть метода умещается в этих строчках:

for i := 0 to len - 1 do
 for j := i to len - 1 do
  if abs(signal2[i]-signal2[j])<radius then
   begin
    dt:=abs(i-j);
    hist[dt]:=hist[dt]+1;
   end;

В добавок к этому все вычисления ведутся в INTEGER без использования плавающей точки!
Единственное место, где есть плавающая точка - это вычисление центра масс пика, которое выполняется один раз в программе (7 строчек кода).
Это означает практически отсутствие ошибки округления.

Так же в программе теперь практически нет умножений. Единственное место - перевод периода в частоту!

ФлеймФорумОбщее

Тема в архиве.