Вопрос по алгоритму БПФ [Закрыто]

dimon72

Требуется получить спектр некоего сигнала.
Я имею алгоритм, который делает БПФ.
На вход алгоритма подаю динамический массив SQRe.
В результате работы алгоритма получаю выходные векторы - динамические массивы
SpRe и SpIm.
Вопрос: SpIm меня не интересует, мне нужна информация об амплитуде спектра, т.е.,
меня интересует массив SpRe. Почему-то его значения получаются в каких-то непонятных
единицах.
Я что-то не могу понять. Вроде я должен получить набор коэффициентов 0...1 для заданной
сетки частот. Например, беру входной сигнал 10Гц, 16 бит со значениями, изменяющимися на полную шкалу от -32768 до +32768. После БПФ я, по идее, должен получить пик на 10 Гц, который должен быть амплитудой +32768. (Я имею в виду действительную часть, её модуль; считаю также, что я попал в сетку частот преобразования).
В результате же использования приведённого ниже алгоритма я получаю с виду вроде бы нормальный спектр, но его значения по Y выражены в каких-то огромных числах. Непонятно откуда они взялись. Кто-нибудь знает что это за хрень?
P.S. Просьба не пинать, т.к. я гуманитарий и основы ПФ знаю очень смутно.
Вот сам алгоритм:
Переменные:


type

  TDigits = array of Double;

var

   // Массивы для хранения входных данных, выходных данных, входного и выходного спектров

    SQRe, SQIm, SpRe,SpIm,SpMod: TDigits;

    SQReO, SQImO, SpReO,SpImO,SpModO: TDigits;

  

//////////////////////////////////////////////////////////////

// Реализует алгоритм быстрого прямого преобразования Фурье //

//////////////////////////////////////////////////////////////

procedure TLowPassFilter.Bpf;

var

  xC1re,xC1im,xC2re,xC2im,xVre,xVim : double;

  i,j,k : integer;

  xMm,xLl,xJj,xKk,xNn,xNv2,xNm1: integer;

  xCounter, xNecessary: Integer;

begin

// copy

  for i := 0 to FSpectrumCount - 1 do

  begin

    SpIm[i] := 0;

    SpMod[i] := 0;

    if i < FSampleCount then

      SpRe[i]:=SQRe[i]

    else

      SpRe[i]:=0;

  end;



  //Begin

  xCounter := 0;

  xNecessary := Round(FSpectrumCount * log2(FSpectrumCount;

  xMm:=1;

  xLl:=FSpectrumCount;



  // Внешний цикл для слоев Nexp

  for k := 1 to Nexp do

  begin

    xNn:=xLl div 2;

    xJj:=xMm+1;

   // Переупорядочивание и предварительные вычисления

    i:=1;

    while i <= FSpectrumCount do

    begin

      xKk := i + xNn;

      xC1re := SpRe[i-1] + SpRe[xKk-1];

      xC1im := SpIm[i-1] + SpIm[xKk-1];

      SpRe[xKk-1] := SpRe[i-1] - SpRe[xKk-1];

      SpIm[xKk-1] := SpIm[i-1] - SpIm[xKk-1];

      SpRe[i-1] := xC1re;

      SpIm[i-1] := xC1im;

      i := i + xLl;

    end;



    if xNn <> 1 then

    begin

     // Двухточечное преобразование Фурье

      for j := 2 to xNn do

      begin

        xC2re := Cos(TwoPi * (xJj-1) / FSpectrumCount);

        xC2im := -Sin(TwoPi * (xJj-1)/ FSpectrumCount);

        xCounter := xCounter + 2;

        i := j;

        while i <= FSpectrumCount do

        begin

          xKk := i + xNn;

          xC1re := SpRe[i-1] + SpRe[xKk-1];

          xC1im := SpIm[i-1] + SpIm[xKk-1];

          xVre := (SpRe[i-1] - SpRe[xKk-1]) * xC2re - (SpIm[i-1] - SpIm[xKk-1]) * xC2im;

          xVim := (SpRe[i-1] - SpRe[xKk-1]) * xC2im + (SpIm[i-1] - SpIm[xKk-1]) * xC2re;

          SpRe[xKk-1] := xVre;

          SpIm[xKk-1] := xVim;

          SpRe[i-1] := xC1re;

          SpIm[i-1] := xC1im;

          i := i+xLl;



          xCounter := xCounter + 2;

          if (xCounter mod (xNecessary div 20) = 0) then

            begin

              Application.ProcessMessages;

              if FBreakExecute then exit;

              if Assigned(FOnStepExecute) then

                FOnStepExecute(Self, xCounter / xNecessary * 33);

            end;

        end;



        xJj := xJj + xMm;

      end;



      xLl := xNn;

      xMm := xMm * 2;

    end;



  end;



  // Восстановление нужной последовательности

  xNv2 := FSpectrumCount div 2;

  xNm1 := FSpectrumCount - 1;

  j := 1;



  for i := 1 to xNm1 do

  begin

    if i < j then

    begin

      xC1re := SpRe[j-1];

      xC1im := SpIm[j-1];

      SpRe[j-1] := SpRe[i-1];

      SpIm[j-1] := SpIm[i-1];

      SpRe[i-1] := xC1re;

      SpIm[i-1] := xC1im;

    end;

    k := xNv2;

    while k < j do

    begin

      j := j - k;

      k := k div 2;

    end;

    j := j + k;



  end;



  // Вычисление модуля

  SumSpectrum := 0;

  SquareSumSpectrum :=0;

  for i := 0 to FSpectrumCount - 1 do

  begin

    SpMod[i] := Sqrt(SpRe[i] * SpRe[i] + SpIm[i] * SpIm[i]);

    if i<>0 then

    begin

      SumSpectrum := SumSpectrum + SpMod[i];

      SquareSumSpectrum := SquareSumSpectrum + SpMod[i] * SpMod[i];

    end;

  end;

end;

SPARTAK3959

Если сигнал идет с амплитудой A, то на выходе будет максимум величины A*n, где n - размер массива. По крайней мере большинство реализаций ведут себя так.

dimon72

Вот пример.
Имеется тестовый сигнал, являющийся суммой синусоид 10Гц и 10.5Гц.
Максимум амплитуды сигнала отмечен маркером и показан слева.
Полученный спектр виден справа. Видно два пика, соответствующие 10 и 10.5 Гц. Здесь всё ok.
Однако непонятно, в каких единицах и что означают числа по оси ординат?

dimon72

Если сигнал идет с амплитудой A, то на выходе будет максимум величины A*n, где n - размер массива. По крайней мере большинство реализаций ведут себя так.

Пробовал подавать на вход сигналы с разной частотой. n был одинаков, а выходные значения
амплитуды спектра различались.

dimon72

Если сигнал идет с амплитудой A, то на выходе будет максимум величины A*n, где n - размер массива. По крайней мере большинство реализаций ведут себя так.

Что-то не сходится.
В приведённом выше примере имеется сигнал со следующими парамерами:
кол-во отсчётов - 15000 (60 сек
при частоте дискретизации, Fs=250Гц;
максимум амплитуды сигнала: +32734;
минимум амплитуды сигнала: -32734;
Максимальное отмеченное значение амплитуды спектра получилось равным величине:
119 670 305, 89.

dimon72

Вообще, похоже, что на выходе получается A*n/2, где A-амплитуда входного сигнала, а n-размер массива, дополненный до степени двойки, т.е. если берётся 15000 отсчётов, получается A*16384/2=A*8192
Надо проверить, чтобы тестовый сигнал не "расплывался" по выбранной сетке частот, чтобы вся спектральная мощность ушла в один пик.
Есть ещё мнения?
Если сигнал идет с амплитудой A, то на выходе будет максимум величины A*n, где n - размер массива. По крайней мере большинство реализаций ведут себя так.

dimon72

Всё, разобрался!
Итак, величина по Y выходного спектра равна произведению амплитуды входного сигнала на
половину эпохи БПФ преобразования.
Y=A*FFT/2
Всем спасибо!
Отдельное спасибо !
ЗЫ: А разработчики могли бы пояснить всё это, например в описании компонента. Тем более, что контора-то известная (www.basegroup.ru).
Оставить комментарий
Имя или ник:
Комментарий: