2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




 
 fft с параметрами FORTRAN
Сообщение12.07.2019, 13:09 
Не могу понять как изменить обычный 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
 

 
 
 
 Re: fft с параметрами FORTRAN
Сообщение12.07.2019, 13:23 
Ну было бы наверно неплохо выписать и то, что вычисляет
и сопоставить эти два станет куда легче даже тем, кто в фортране не разбирается и не хочет искать тот текст и в том тексте.

 
 
 
 Re: fft с параметрами FORTRAN
Сообщение12.07.2019, 13:46 
arseniiv в сообщении #1404730 писал(а):
Ну было бы наверно неплохо выписать и то, что вычисляет
и сопоставить эти два станет куда легче даже тем, кто в фортране не разбирается и не хочет искать тот текст и в том тексте.


Вычисляет 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]

 
 
 
 Re: fft с параметрами FORTRAN
Сообщение12.07.2019, 14:28 
Ну это прекрасно, но если б вы написали формулу, которая этим кодом реализуется, я бы вам мигом написал, что остаётся сделать, а так леееень. :wink:

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

 
 
 
 Re: fft с параметрами FORTRAN
Сообщение15.07.2019, 16:57 
Пока $b=\pm1$, это (плюс-минус) обычное DFT — а вам нужны были именно другие $b$?

 
 
 
 Re: fft с параметрами FORTRAN
Сообщение16.07.2019, 08:22 
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]]

 
 
 
 Re: fft с параметрами FORTRAN
Сообщение26.07.2019, 15:16 
Аватара пользователя
sithif в сообщении #1404726 писал(а):
С fp1 проблем нет. Если просто домножить показатель экспоненты, то результаты не совпадают (но совпадают для нечетных fp2), нужно чтобы совпадали для дюбых значений.
Как правильно добавить fp2?


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

 
 
 [ Сообщений: 8 ] 


Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group