fft с параметрами FORTRAN

[sithif]

Не могу понять как изменить обычный fft (из Numerical Recipes in FORTRAN 77), что бы он работал с параметрами как в Mathematica.
В Mathematica используется определение:

$v_s = \frac{1}{n^{(1-a)/2}}\sum _{r=1}^n u_r e^{(2 \pi b i (r-1) (s-1))/n}$

$u_r$ -- элемент входных данных, $v_s$ -- преобразование $s$-того элемента, $a$ и $b$ -- параметры (в коде fp1 и fp2).

Простой пример в Mathematica:

sig = Sin[2 Pi 0.12345 Range[128]] ;
fou = Fourier[sig,FourierParameters -> {0,2}] ;
fou[[1;;10]]//TableForm

Ответ:

$ \left( \begin{array}{l} 0.031827\, +0. i \\ 0.0314988\, -0.0177738 i \\ 0.0304466\, -0.0372855 i \\ 0.02843\, -0.0609235 i \\ 0.0248761\, -0.0929647 i \\ 0.0183513\, -0.143123 i \\ 0.00433582\, -0.2411 i \\ -0.0418564-0.547854 i \\ 0.851101\, +5.26196 i \\ 0.119906\, +0.495437 i \\ \end{array} \right) $

С fp1 проблем нет. Если просто домножить показатель экспоненты, то результаты не совпадают (но совпадают для нечетных fp2), нужно чтобы совпадали для дюбых значений.
Как правильно добавить fp2?

Используется синтаксис Fortran

program test
    implicit none
    integer,  parameter       :: rk=selected_real_kind(15,307)
    real(rk), dimension(256)  :: dat
    real(rk), dimension(128)  :: sig
    integer                   :: i
        
    do i=1,128,1
        sig(i) = sin(6.28318530717959_rk*0.12345_rk*real(i,rk))
    end do

    dat = 0._rk
    dat(1:256:2) = sig
    ! fft call with a = 0 and b = 2
    call fft(128,dat,0._rk,2.0_rk)
    do i=1,2*10,2
        write(*,*) dat(i), dat(i+1)
    end do
    
    
contains
subroutine fft(num,arr,fp1,fp2)
    implicit none
    integer,  parameter     :: rk=selected_real_kind(15,307)
    integer,  intent(in)    :: num
    real(rk), intent(in)    :: fp1
    real(rk), intent(in)    :: fp2
    real(rk), intent(inout) :: arr(2*num)
    integer                 :: n,i,j,m,lim,ste
    real(rk)                :: pim, pre, ang, wr, wi, wpr, wpi, wx
    n=2*num
    j=1
    do i=1,n,2
        if(j.gt.i)then
            pre=arr(j)
            pim=arr(j+1)
            arr(j)=arr(i)
            arr(j+1)=arr(i+1)
            arr(i)=pre
            arr(i+1)=pim
        endif
        m=n/2
1       if ((m.ge.2).and.(j.gt.m)) then
            j=j-m
            m=m/2
            goto 1
        endif
        j=j+m
    enddo
    lim=2
2   if (n.gt.lim) then
        ste=2*lim
        ang=6.28318530717959_rk/real(lim,rk)*fp2   ! -- just mult 2 pi
        wpr=-2.0_rk*sin(0.5_rk*ang)**2             ! -- ? *sign(1._rk,ang)
        wpi=sin(ang)
        wr=1.0_rk
        wi=0.0_rk
        do m=1,lim,2
            do i=m,n,ste
                j=i+lim
                pre=wr*arr(j)-wi*arr(j+1)
                pim=wr*arr(j+1)+wi*arr(j)
                arr(j)=arr(i)-pre
                arr(j+1)=arr(i+1)-pim
                arr(i)=arr(i)+pre
                arr(i+1)=arr(i+1)+pim
            enddo
            wx=wr
            wr=wr*wpr-wi*wpi+wr
            wi=wi*wpr+wx*wpi+wi
        enddo
        lim=ste
        goto 2
    endif
    arr = arr/real(num,rk)**((1.0_rk-fp1)/2.0_rk)
    return
end subroutine fft

end program test
 

[arseniiv]

Ну было бы наверно неплохо выписать и то, что вычисляет

sithif в сообщении #1404726 писал(а):

fft (из Numerical Recipes in FORTRAN 77)

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

[sithif]

arseniiv в сообщении #1404730 писал(а):

Ну было бы наверно неплохо выписать и то, что вычисляет

sithif в сообщении #1404726 писал(а):

fft (из Numerical Recipes in FORTRAN 77)

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

Вычисляет subroutine fft(num,arr,fp1,fp2), в исходном сообщении, почти без изменений из Numerical Recipes in FORTRAN 77

Вот код в Mathematica для $b=1$, для $b=2$ не совпадают уже.

fp1 = 0. ; (* -- set a *)
fp2 = 1. ; (* -- set b *)

Block[
   {},
   
   sig = ConstantArray[0.,128] ;
   Do[
      sig[[i]] = Sin[6.28318530717959*0.12345*i]
      ,{i,1,128,1}
   ] ;

   dat = ConstantArray[0.,2 128] ;
   dat[[1;;256;;2]] = sig ;
   num = Length[sig] ;
   n=2*num ;
   j = 1 ;

   Do[
      If[
         j > i,
         pre = dat[[j]] ;
         pim = dat[[j+1]] ;
         dat[[j]] = dat[[i]] ;
         dat[[j+1]] = dat[[i+1]] ;
         dat[[i]] = pre ;
         dat[[i+1]] = pim ;
      ] ;
      m = n/2 ;
      Label[1] ;
      If[
         m >= 2 && j > m,
         j = j - m ;
         m = m/2 ;
         Goto[1] ;
      ] ;
      j = j + m
      ,{i,1,n,2}
   ] ;
   lim = 2 ;
   Label[2] ;
   If[  
      n > lim,
      ste=2*lim ;
      ang=6.28318530717959/lim*fp2 ; (* -- mut pi by b *)
      wpr=-2.0*Sin[0.5*ang]^2 ;
      wpi=Sin[ang] ;
      wr=1.0 ;
      wi=0.0 ;
      Do[
         Do[
            j = i + lim ;
            pre = wr*dat[[j]]-wi*dat[[j+1]] ;
            pim = wr*dat[[j+1]]+wi*dat[[j]] ;
            dat[[j]]=dat[[i]]-pre;
            dat[[j+1]]=dat[[i+1]]-pim ;
            dat[[i]]=dat[[i]]+pre ;
            dat[[i+1]]=dat[[i+1]]+pim
            ,{i,m,n,ste}
         ] ;
         wx=wr ;
         wr=wr*wpr-wi*wpi+wr ;
         wi=wi*wpr+wx*wpi+wi ;
      ,{m,1,lim,2}
      ] ;
      lim=ste ;
      Goto[2] ;
   ] ; 
   dat = dat/num^((1.0-fp1)/2.0) ;
] ;
res = Map[Apply[Complex],Partition[dat,2]] ; res[[1;;5]]
ref = Fourier[sig,FourierParameters->{0,1}] ; ref[[1;;5]]
res-ref//Chop//DeleteCases[0]

[arseniiv]

Ну это прекрасно, но если б вы написали формулу, которая этим кодом реализуется, я бы вам мигом написал, что остаётся сделать, а так леееень.

[sithif]

Оказывается, функция Fourier в WM с FourierParameters вычисляет linear fractional Fourier.
И нельзя "обобщить" обычный fft, но fft может быть использован для вычисления (см., например, Bailey).

[arseniiv]

Пока $b=\pm1$, это (плюс-минус) обычное DFT — а вам нужны были именно другие $b$?

[sithif]

arseniiv в сообщении #1405191 писал(а):

Пока $b=\pm1$, это (плюс-минус) обычное DFT — а вам нужны были именно другие $b$?

Мне нужно конкретное значение для зума между модами, как в примерах WM по Fourier (Applications/Frequency Identification).
Но, в принципе, должно работать для любых значений $b$.
WM код для вычисления ffrft через fft на основе статьи Bailey, его легко переписать в FORTRAN:

(* ffrft with normal fft *)
(* https://www.davidhbailey.com/dhbpapers/fracfft.pdf *)
(* Set data length *)
m = 1024 ;
(* Generate test data *)
dat = Sin[2 Pi EulerGamma/Pi Range[m]] ;
(* Set fractional Fourier parameters *)
a = 1 ;   (* -- scaling *)
b = 2/m ; (* -- any complex *)
(* Compute with Fourier[] *)
ref = Fourier[dat,FourierParameters->{1,b}] ;
(* Bailey's method p = m & flip sign for i to match Fourier[] notation with direct fft defined with {a,b} = {1,1} *)
alpha = b/m ;
(* Generate extended arrays *)
arr1 = arr2 = ConstantArray[0.,2 m] ;
arr1[[1;;m]] = dat N[Exp[Pi I Range[0,m-1]^2 alpha]] ;
arr1[[m+1;;-1]] = 0. ;
arr2[[1;;m]] = N[Exp[- Pi I Range[0,m-1]^2 alpha]] ;
arr2[[m+1;;-1]] = N[Exp[-Pi I (Range[m,2 m-1]- 2m )^2 alpha]] ;
(* Use normal fft (emulated with Fourier with {a,b} = {1,1}) *)
fou1 = Fourier[arr1,FourierParameters->{1,1}] ;
fou2 = Fourier[arr2,FourierParameters->{1,1}] ;
arr3 = fou1 fou2 ;
(* Inverse fft (only 1st m points make sense) *)
fou3 = Fourier[arr3,FourierParameters->{1,-1}] ;
fou3 = Take[fou3,m] ;
(* Result  *)
res = N[Exp[Pi I Range[0,m-1]^2 alpha]] fou3 /(2 m);
(* Compare *)
ref-res // Chop // DeleteCases[0]
ListPlot[Abs[ref]]
ListPlot[Abs[res]]

[Snegovik]

sithif в сообщении #1404726 писал(а):

С fp1 проблем нет. Если просто домножить показатель экспоненты, то результаты не совпадают (но совпадают для нечетных fp2), нужно чтобы совпадали для дюбых значений.
Как правильно добавить fp2?

Как то не совсем понятно описали проблему. У Вас и fp1 и fp2 это параметры которые передаются в подпрограмму и оба используются. В формуле это параметры $a$ и $b$, насколько я понял. Что значит - "как правильно добавить fp2"?