2014 dxdy logo

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

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




Начать новую тему Ответить на тему
 
 fft с параметрами FORTRAN
Сообщение12.07.2019, 13:09 


08/03/11
186
Не могу понять как изменить обычный 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 
Заслуженный участник


27/04/09
28128
Ну было бы наверно неплохо выписать и то, что вычисляет
и сопоставить эти два станет куда легче даже тем, кто в фортране не разбирается и не хочет искать тот текст и в том тексте.

 Профиль  
                  
 
 Re: fft с параметрами FORTRAN
Сообщение12.07.2019, 13:46 


08/03/11
186
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 
Заслуженный участник


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

 Профиль  
                  
 
 Re: fft с параметрами FORTRAN
Сообщение15.07.2019, 13:57 


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

 Профиль  
                  
 
 Re: fft с параметрами FORTRAN
Сообщение15.07.2019, 16:57 
Заслуженный участник


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

 Профиль  
                  
 
 Re: fft с параметрами FORTRAN
Сообщение16.07.2019, 08:22 


08/03/11
186
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 
Аватара пользователя


30/04/19
235
sithif в сообщении #1404726 писал(а):
С fp1 проблем нет. Если просто домножить показатель экспоненты, то результаты не совпадают (но совпадают для нечетных fp2), нужно чтобы совпадали для дюбых значений.
Как правильно добавить fp2?


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

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 8 ] 

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group