На главную Наши проекты:
Журнал   ·   Discuz!ML   ·   Wiki   ·   DRKB   ·   Помощь проекту
ПРАВИЛА FAQ Помощь Участники Календарь Избранное RSS
msm.ru
! ПРАВИЛА РАЗДЕЛА · FAQ раздела Delphi · Книги по Delphi
Пожалуйста, выделяйте текст программы тегом [сode=pas] ... [/сode]. Для этого используйте кнопку [code=pas] в форме ответа или комбобокс, если нужно вставить код на языке, отличном от Дельфи/Паскаля.
Следующие вопросы задаются очень часто, подробно разобраны в FAQ и, поэтому, будут безжалостно удаляться:
1. Преобразовать переменную типа String в тип PChar (PAnsiChar)
2. Как "свернуть" программу в трей.
3. Как "скрыться" от Ctrl + Alt + Del (заблокировать их и т.п.)
4. Как прочитать список файлов, поддиректорий в директории?
5. Как запустить программу/файл?
... (продолжение следует) ...

Вопросы, подробно описанные во встроенной справочной системе Delphi, не несут полезной тематической нагрузки, поэтому будут удаляться.
Запрещается создавать темы с просьбой выполнить какую-то работу за автора темы. Форум является средством общения и общего поиска решения. Вашу работу за Вас никто выполнять не будет.


Внимание
Попытки открытия обсуждений реализации вредоносного ПО, включая различные интерпретации спам-ботов, наказывается предупреждением на 30 дней.
Повторная попытка - 60 дней. Последующие попытки бан.
Мат в разделе - бан на три месяца...
Модераторы: jack128, D[u]fa, Shaggy, Rouse_
  
> Обратное преобразование Фурье (IFFT) , Не корректно работает
    Доброго времени суток. Разбираюсь с FFT и IFFT. Уже какой день парюсь над обратным преобразованием Фурье. Добрые люди, может объясните в чем может быть проблема?

    Звуковой сигнал 2 кГц, который имитирует программка справа, поступает на микрофон. Данные заносятся в массив и рисуется входной сигнал - верхняя диаграмма.

    Этот массив с данными раскладываю в FFT, получаю массив fft и рисую диаграмму - второй график сверху. Проверял на разных частотах, все здесь правильно.

    И наконец, я делаю обратное FFT, взяв массив fft. В итоге, 3-й график, вроде и частота сигнала та же, что и на входе, но почему-то блин амплитуда сигнала падает. И видно,что на стыке двух спектров (основного и зеркального) амплитуда минимальна.

    Делал эксперимент, если в спектре четко только одна гармоника, то при IFFT все нормально, амплитуда сигнала равномерна по всей длине, как на входном графике.
    Как только в спектре есть две гармоники и более, то сразу же график как на картинке. Как будто перемножение что ли где-то, или что может быть?

    Буду благодарен знатокам!

    ExpandedWrap disabled
      procedure TFFTBase.InitMem(WindowSize: Integer);
      begin
       fWindowSize:=WindowSize;
       tTlbSize:=(2*fWindowSize+1)*SizeOf(Single);
       fSinTbl:=AllocMem(tTlbSize);
       fCosTbl:=AllocMem(tTlbSize);
      end;
       
      procedure TFFTBase.InitSinCosTbl;
      var i :integer;
      begin
       for i :=1 to 2*fWindowSize do begin
        fSinTbl[i] := (-1)*Sin(PI/i);
        fCosTbl[i] := Cos(PI/i);
       end;
      end;
       
       NV2 :=N shr 1;
       NM1 :=N-1;
       J   :=1;
       if Inverse then
        for i :=0 to N-1 do
         TComplexArray(F^)[i].Im :=-TComplexArray(F^)[i].Im;
       
       for I :=1 to NM1 do begin
        if I<J then begin
         T      :=TComplexArray(F^)[J-1];
         TComplexArray(F^)[J-1] :=TComplexArray(F^)[I-1];
         TComplexArray(F^)[I-1] :=T;
        end;
        K :=NV2;
        while K < J do begin
         J :=J - K;
         K :=K shr 1;
        end;
        J :=J + K;
       end;
       for L :=1 to M do begin
        LE   :=2 shl (L-1);
        LE1  :=LE shr 1;
        U.Re :=1.0; U.Im :=0.0;
        W.Re :=fCosTbl[LE1];
        W.Im :=fSinTbl[LE1];
        for J :=1 to LE1 do begin
         I :=J;
         while I <= N do begin
          IP :=I + LE1;
          T.Re    :=TComplexArray(F^)[IP-1].Re * U.Re - TComplexArray(F^)[IP-1].Im * U.Im;
          T.Im    :=TComplexArray(F^)[IP-1].Re * U.Im + TComplexArray(F^)[IP-1].Im * U.Re;
          TComplexArray(F^)[IP-1].Re :=TComplexArray(F^)[I-1].Re - T.Re;
          TComplexArray(F^)[IP-1].Im :=TComplexArray(F^)[I-1].Im - T.Im;
          TComplexArray(F^)[I-1].Re:=TComplexArray(F^)[I-1].Re+T.Re;
          TComplexArray(F^)[I-1].Im:=TComplexArray(F^)[I-1].Im+T.Im;
          Inc(I,LE);
         end;
         Uo :=U;
         U.Re :=(Uo.Re * W.Re) - (Uo.Im * W.Im);
         U.Im :=(Uo.Re * W.Im) + (Uo.Im * W.Re);
        end;
       end;
       
       ImDummy :=1/Sqrt(N);
       if Inverse then ImDummy :=-ImDummy;
       ReDummy :=Abs(ImDummy);
       for I :=1 to N do
        begin
         TComplexArray(F^)[i-1].Re :=TComplexArray(F^)[i-1].Re*ReDummy;
         TComplexArray(F^)[i-1].Im :=TComplexArray(F^)[i-1].Im*ImDummy;
        end;
       
       DelMem;
      end;

    Прикреплённый файлПрикреплённый файл111.png (53,48 Кбайт, скачиваний: 555)
      Здесь не приведен заголовок процедуры FFT. Как она применяется? Каков размер массива?

      Какие частоты вызывают такой эффект?

      Получается ли то же самое с достоверно рабочим кодом FFT?
        Цитата MBo @

        Применен типовой код:

        ExpandedWrap disabled
          unit FFTBase;
          interface
          uses
            Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
            FFTFilterConst;
           
          {$R-}
          type
            TFFTBase = class(TComponent)
            private
              { Private declarations }
             fWindowSize, tTlbSize: Integer;
             fSinTbl, fCosTbl     : PTbl;
             procedure InitSinCosTbl;
            protected
              { Protected declarations }
            public
              { Public declarations }
             procedure FFT( var F        :Pointer;
                    N, M          :integer;
                    const Inverse :boolean;
                    const Window  :Integer ); //N -размер F, M - N =2^M
             procedure simpleFFT( var F        :Pointer;
                    N, M          :integer;
                    const Inverse :boolean;
                    const Window  :Integer ); //N -размер F, M - N =2^M
             procedure InitMem(WindowSize: Integer);
             procedure DelMem;
            published
              { Published declarations }
            end;
           
          procedure Register;
           
          implementation
           
          {+private++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++}
          procedure TFFTBase.InitSinCosTbl;
          var i :integer;
          begin
           for i :=1 to 2*fWindowSize do begin
            fSinTbl[i] := (-1)*Sin(PI/i);
            fCosTbl[i] := Cos(PI/i);
           end;
          end;
          {-private----------------------------------------------------------------------}
           
          {+PUBLIC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++}
          procedure TFFTBase.FFT( var F       :Pointer;
                              N, M          :integer;
                              const Inverse :boolean;
                      const Window  :Integer ); //N -размер F, M - N =2^M
          var U, Uo, W, T :TComplex;
              I,IP,J,K,L,LE,LE1,NV2,NM1 :integer;
              ReDummy, ImDummy :extended;
           
              procedure SetWindow(const Window: Integer);
              var i: Integer;
              begin
               if Not Inverse then //Если не обратный
                case Window of  //0-прямоугольное//1-треугольное//2-Хэминга//3-Ханна//4-Блэкмана
                 0: begin  //0-прямоугольное
                     {for i :=0 to N-1 do
                       F[i].Re :=F[i].Re*1.2}
                    end;
                 1: begin //1-треугольное
                     for i :=0 to N-1 do begin
                      if (i <= N div 2) then
                       TComplexArray(F^)[i].Re :=TComplexArray(F^)[i].Re*i/(N/2)
                      else TComplexArray(F^)[i].Re :=TComplexArray(F^)[i].Re*(N-i)/(N/2);
                     end;
                    end;
                 2: begin //2-Хэминга
                     for i :=0 to N-1 do
                   TComplexArray(F^)[i].Re :=TComplexArray(F^)[i].Re*(0.54-0.46*cos(2*Pi*i/N))
                    end;
                 3: begin //3-Ханна
                     for i :=0 to N-1 do
                  TComplexArray(F^)[i].Re :=TComplexArray(F^)[i].Re*(1+cos(2*Pi*i/N))/2;
                    end;
                 4: begin //4-Блэкмана
                 for i :=0 to N-1 do
                      TComplexArray(F^)[i].Re :=TComplexArray(F^)[i].Re*(0.42+0.5*cos(2*Pi*i/N)+0.08*cos(4*Pi*i/N))
                    end;
                end;
              end;
          begin
           
           InitMem(N);
           InitSinCosTbl;
           SetWindow(Window);
           
           NV2 :=N shr 1;
           NM1 :=N-1;
           J   :=1;
           if Inverse then
            for i :=0 to N-1 do
             TComplexArray(F^)[i].Im :=-TComplexArray(F^)[i].Im;
           
           for I :=1 to NM1 do begin
            if I<J then begin
             T      :=TComplexArray(F^)[J-1];
             TComplexArray(F^)[J-1] :=TComplexArray(F^)[I-1];
             TComplexArray(F^)[I-1] :=T;
            end;
            K :=NV2;
            while K < J do begin
             J :=J - K;
             K :=K shr 1;
            end;
            J :=J + K;
           end;
           for L :=1 to M do begin
            LE   :=2 shl (L-1);
            LE1  :=LE shr 1;
            U.Re :=1.0; U.Im :=0.0;
            W.Re :=fCosTbl[LE1];
            W.Im :=fSinTbl[LE1];
            for J :=1 to LE1 do begin
             I :=J;
             while I <= N do begin
              IP :=I + LE1;
              T.Re    :=TComplexArray(F^)[IP-1].Re * U.Re - TComplexArray(F^)[IP-1].Im * U.Im;
              T.Im    :=TComplexArray(F^)[IP-1].Re * U.Im + TComplexArray(F^)[IP-1].Im * U.Re;
              TComplexArray(F^)[IP-1].Re :=TComplexArray(F^)[I-1].Re - T.Re;
              TComplexArray(F^)[IP-1].Im :=TComplexArray(F^)[I-1].Im - T.Im;
              TComplexArray(F^)[I-1].Re:=TComplexArray(F^)[I-1].Re+T.Re;
              TComplexArray(F^)[I-1].Im:=TComplexArray(F^)[I-1].Im+T.Im;
              Inc(I,LE);
             end;
             Uo :=U;
             U.Re :=(Uo.Re * W.Re) - (Uo.Im * W.Im);
             U.Im :=(Uo.Re * W.Im) + (Uo.Im * W.Re);
            end;
           end;
           
           ImDummy :=1/Sqrt(N);
           if Inverse then ImDummy :=-ImDummy;
           ReDummy :=Abs(ImDummy);
           for I :=1 to N do
            begin
             TComplexArray(F^)[i-1].Re :=TComplexArray(F^)[i-1].Re*ReDummy;
             TComplexArray(F^)[i-1].Im :=TComplexArray(F^)[i-1].Im*ImDummy;
            end;
           
           if not Inverse then //Вычисление Мощности С=SQRT(А^2+B^2) при -прямом- БПФ
            for i :=0 to N-1 do
             TComplexArray(F^)[i].Re :=sqrt(sqr(TComplexArray(F^)[i].Re)
                                           +sqr(TComplexArray(F^)[i].Im));
           
           DelMem; //Освобождение ресурсов
          end;
           
          procedure TFFTBase.InitMem(WindowSize: Integer);
          begin
           fWindowSize:=WindowSize;
           tTlbSize:=(2*fWindowSize+1)*SizeOf(Single);
           fSinTbl:=AllocMem(tTlbSize);
           fCosTbl:=AllocMem(tTlbSize);
          end;
           
          procedure TFFTBase.DelMem;
          begin
           FreeMem(fCosTbl, tTlbSize);
           FreeMem(fSinTbl, tTlbSize);
          end;


        Так происходит при любых частотах. BufSize = N = 1024. Другие FFT еще не пробовал...
        Из моей картинки, как я понял, на втором и третьем графиках только половины сигналов спектра и IFFT, правые части - зеркальны. Восстановить целиком сигнал не выйдет, получается?...

        Запускается этот FFTBase отсюда, после комментария // 4-Блэкмана:

        ExpandedWrap disabled
          procedure MakeIFFT;               //Делаем IFFT
          var
           fftb: TFFTBase; //класс, который реализует БПФ
           fFFTComplBuf: ^TComplexArray;  //Буфер для хранения комплексных величин
           i: integer;
          begin
           GetMem(fFFTComplBuf, BufSize*SizeOf(TComplex)); //Выделение памяти под массив
           for i:=0 to BufSize-1 do //Заполняем данными массив
           begin
            fFFTComplBuf[i].Re := fDataBuf[i];
            fFFTComplBuf[i].Im := 0;
           end;
           fftb:=TFFTBase.Create(nil);
          //FFT - выполнение БПФ
          //Параметры:
          //указатель на массив данных (с комплексными числами)
          //N - количество данных (размерность массива)
          //2^X=N (степень числа два)
          //False – прямое преобразование; True – обратное
          //Тип окна:
          // 0-прямоугольное
          // 1-треугольное
          // 2-Хэминга
          // 3-Ханна
          // 4-Блэкмана
           fftb.FFT(Pointer(fFFTComplBuf), BufSize, 10, True, 2);
           for i:=0 to BufSize-1 do //Переносим результат БПФ в исходный массив
           begin
            //заполняем массив выходными значениями (предварительно масштабируем их)
            fDataBuf[i] := Round(fFFTComplBuf[i].Re / 1);
           end;
           fftb.Free;
           FreeMem(fFFTComplBuf, BufSize*SizeOf(TComplex)); //Освобождение памяти выделенной под массив
          end;


        И выполняются процедуры в первом коде по очереди:
        InitMem(N);
        InitSinCosTbl;
        SetWindow(Window);
        затем само Фурье -
        NV2 :=N shr 1;
        NM1 :=N-1;
        J :=1; и т.д.

        Результат переносится в fDataBuf[i] во втором коде.
        Сообщение отредактировано: ershovser -
          Так же на 10 кГц, видно на картинке...

          На всех частотах. Однако, к примеру на 5 кГц, я в массиве обнулил небольшие гармоники, оставил только одну - самую большую - все нормально, сигнал на выходе такой же как на входе.

          Может что из-за окна? При FFT применяю окно 2-Хэмминга. Очень похоже. А при IFFT окна не применяются, исходя из кода.
          Прикреплённый файлПрикреплённый файл222.png (40,23 Кбайт, скачиваний: 638)
          Прикреплённый файлПрикреплённый файл333.png (49,93 Кбайт, скачиваний: 620)
            Сделал сигнал длиной 1024 отсчета, ввел одну, две, три частоты.
            Не наблюдаю такого эффекта как с окном, так и без окна. Может быть, выводятся графики не того, что есть на самом деле?


            ExpandedWrap disabled
              var
                Re, Im: array of Single;
                i, N: Integer;
                Amp: Double;
              begin
                N := 1024;
                SetLength(Re, N);
                SetLength(Im, N);
                for i := 0 to N - 1 do begin
                  Re[i] := Sin(50 * i * 2 * Pi / N) +
                           0.3 * Sin(67 * i * 2 * Pi / N) +
                           0.3 * Sin(91 * i * 2 * Pi / N);
                 //окно, его можно отключить
                  Re[i] := Re[i] * (0.54-0.46*cos(2*Pi*i/N));
                  Im[i] := 0;
                  Series1.AddXY(i, Re[i]);
                end;
               
                FFT1D(Re, Im, N, 1);
                for I := 0 to N - 1 do begin
                  Amp := Hypot(Re[i], Im[i]);
                  Series2.AddXY(i, Amp);
                end;
               
                FFT1D(Re, Im, N, 0);  //обратное
                for I := 0 to N - 1 do begin
                  Series4.AddXY(i, Re[i] + 3.1); //сдвинем для удобства
                end;

            Прикреплённый файлПрикреплённый файлfft.png (13,99 Кбайт, скачиваний: 597)
            Сообщение отредактировано: MBo -
              Скорее всего у меня выводится сигнал зеркально и выходит, что левая и правая половинки вместе образуют по середине провал по амплитуде... У Вас амплитуда по краям минимальная, а у меня наоборот - посередине.
              И все же, если поставить окно прямоугольного типа, то эти провалы куда меньше. Найду другой FFT проверю как там.

              Спасибо. :good:
              Сообщение отредактировано: ershovser -
                Амплитуда минимальная по краям из-за применения окна
                Сообщение отредактировано: MBo -
                0 пользователей читают эту тему (0 гостей и 0 скрытых пользователей)
                0 пользователей:


                Рейтинг@Mail.ru
                [ Script execution time: 0,0407 ]   [ 21 queries used ]   [ Generated: 28.03.24, 08:02 GMT ]