Фурье спектр изображения

[triest]

Необходимо сделать фурье преобразование изображения. Вначале прохожу по строкам, затем по столбцам промежуточного результата.

Для одномерного массива:

For j:=0 to M-1 do begin
For i:=0 to M-1 do begin
f_x_Re[j]:=f_x_Re[j]+Re[i]*cos((2*Pi*i*j)/M);
f_x_Im[j]:=f_x_Im[j]+Im[i]*sin((2*Pi*i*j)/M);
end;
end;

Для двухмерного

забираем пиксели в массив

 for i:=0 to Hi-1 do  //по строкам
    begin
    for j:=0 to Wi-1 do   //столбцам
      begin
        fxy[i,j] :=image.Canvas.Pixels[j,i];     //
         fxy[i,j] :=fxy[i,j]*power(-1,i+j); //    умножаем на -1 в степени i+j
    end;

Сначала делаю преобразования по строкам
//первое число массива строка
// второе столбец

  for i:=0 to Hi-1 do
        begin
            for j:=0 to Wi-1 do
              begin
                  for k:=0 to Wi-1 do
                    begin
                       fxv_re[i,j]:=fxv_re[i,j]+round(fxy[i,k]*cos((2*pi*j*k)/Wi));
                       fxv_im[i,j]:=fxv_im[i,j]+round(fxy[i,k]*sin((2*pi*j*k)/Wi));
                    end;
              end;
        end;

Затем по слолбцам промежуточного результата, но с учетом комплексной части

 for j:=0 to Wi-1 do
    begin
      for i:=0 to Hi-1 do
        begin
          for k:=0 to Hi-1 do
            begin
                fuv_re[i,j]:=fuv_re[i,j]+fxv_re[k,j]*cos((2*pi*i*k)/Hi)+fxv_im[k,j]*sin((2*pi*j*k)/Hi);
                fuv_im[i,j]:=fuv_im[i,j]+fxv_re[k,j]*sin((2*pi*i*k)/Hi)-fxv_im[k,j]*cos((2*pi*j*k)/Hi);
            end;
        end
    end;

Затем выделяю действительную часть по формуле Эйлера.

        for i:=0 to Hi-1 do
        begin
            for j:=0 to Wi-1 do
              begin
                fuv[i,j]:=sqrt(sqr(fuv_re[i,j])+sqr(fuv_im[i,j]));
              end;
        end;

Выводим результат

form8.visible:=true;
  for i:=0 to Hi-1 do  //
    begin
    for j:=0 to Wi-1 do   //
      begin
      form8.image1.Canvas.Pixels[j,i]:=round(fuv[i,j]);   //выводим результат
      end;
    end;

Но вот результат не совпадает с тем. что предполагается (онлайн преобразование)
http://www.ejectamenta.com/Imaging-E...filtering.html
Код функции, Spectr, запускающейся при нажатии кнопки "прямое гауса" (Unut1).

Полный проект.
http://file.sampo.ru/dw2rvv/

Где ошибка?

procedure FilterSetupClass.Spectr (image: TImage);
var
  i,j,u,v,x,y,M,N,pixeli,pixelj ,t,p,k,z,Height,Width,sizex,sizey,Wi,Hi:integer;
  pwr,D,H,H1,H2: real;
  summa:Longint;
  tempr, tempg, tempb, ftemp: extended;
  filteredPixel: TColor;
  fxy, fuv, fuv_im, fuv_re, fxv, fxv_re, fxv_im: array of array of real;
 const nn=3;
begin

 SetLength(fxy, image.Height);    //массив предварительной обработки ?
  SetLength(fuv, image.Height);    //массив преобразования фурье
  SetLength(fuv_re, image.Height);  //массив действительной части
  SetLength(fuv_im, image.Height);  //массив мнимой части
  SetLength(fxv, image.Height);
  SetLength(fxv_re, image.Height);
  SetLength(fxv_im, image.Height);
  for i:=0 to image.Height-1 do
  begin
    SetLength(fxy[i], image.Width);
    SetLength(fuv[i], image.Width);
    SetLength(fuv_re[i], image.Width);
    SetLength(fuv_im[i], image.Width);
    SetLength(fxv[i], image.Width);
    SetLength(fxv_re[i], image.Width);
    SetLength(fxv_im[i], image.Width);
  end;


  for i:=0 to image.Height-1 do
    begin
         SetLength(fxy[i], image.Width);
    end;

    sizex:=image.height;
    Hi:=image.height;
    sizey:=image.Width;
    Wi:=image.Width;

    // умножаем изображение на -1 в степени x+y
   for i:=0 to Hi-1 do  //высота
    begin
    for j:=0 to Wi-1 do   //ширина
      begin
        fxy[i,j] :=image.Canvas.Pixels[j,i];     //забираем пиксели в массив
         fxy[i,j] :=fxy[i,j]*power(-1,i+j); //    умножаем функцию на -1 в степени 3}
      end;
    end;



    //первое число массива-строка, второе-столбец
     // по стокам

      for i:=0 to Hi-1 do
        begin
            for j:=0 to Wi-1 do
              begin
                  for k:=0 to Wi-1 do
                    begin
                       fxv_re[i,j]:=fxv_re[i,j]+round(fxy[i,k]*cos((2*pi*j*k)/Wi));
                       fxv_im[i,j]:=fxv_im[i,j]+round(fxy[i,k]*sin((2*pi*j*k)/Wi));
                    end;
              end;
        end;
    //

  // по столбцам

  for j:=0 to Wi-1 do
    begin
      for i:=0 to Hi-1 do
        begin
          for k:=0 to Hi-1 do
            begin
                fuv_re[i,j]:=fuv_re[i,j]+fxv_re[k,j]*cos((2*pi*i*k)/Hi)+fxv_im[k,j]*sin((2*pi*j*k)/Hi);
                fuv_im[i,j]:=fuv_im[i,j]+fxv_re[k,j]*sin((2*pi*i*k)/Hi)-fxv_im[k,j]*cos((2*pi*j*k)/Hi);
            end;
        end
    end;

      //  выделяем действ часть
        for i:=0 to Hi-1 do
        begin
            for j:=0 to Wi-1 do
              begin
                fuv[i,j]:=sqrt(sqr(fuv_re[i,j])+sqr(fuv_im[i,j]));
              end;
        end;

       form8.visible:=true;
  for i:=0 to Hi-1 do  //высота
    begin
    for j:=0 to Wi-1 do   //ширина
      begin
      form8.image1.Canvas.Pixels[j,i]:=round(fuv[i,j]);   //выводим результат в форму 8
      end;
    end;
   //   Form7.Visible:=true;
  
end;