Skip to main content
Mathematica

Return to Answer

edited body
Source Link
I.M.
I.M.
  • 3.2k
  • 1
  • 14
  • 19

Here is an implementation based on Bailey's paper.

ClearAll[ffrft$set] ;
ffrft$set[length_, span_] := Block[
{trig, work},
trig =Exp[I*Pi*span/length*Range[0.0, length-1]^2] ;
work = Fourier[Flatten[{Conjugate[trig] , Exp[-I*Pi*span/length*Range[-length, -1.0]^2]}] , FourierParameters -> {1, 1} ];
{trig, work}
] ;

ClearAll[ffrft$get] ;
ffrft$get[length_, signal_, trig_, work_] := Block[
{ffrft},
ffrft = ConstantArray[0.0, 2*length]  ;
ffrft[[1;;length]] = trig*signal ;
ffrft = work*Fourier[ffrft, FourierParameters -> {1, 1} ] ;
ffrft =Fourier[ffrft, FourierParameters -> {1, -1} ][[1;;length]] ;
ffrft*trig/(2*length)
] ;

Example:

length = 1024 ;
signal = Sin[2*Pi*0.12345*Range[length]] ;
b = 1/Sqrt[2.0]  ;
{trig, work} = ffrft$set[length, b]  ;
ref = Fourier[signal, FourierParameters -> {1, b}] ;
res = ffrft$get[length, signal, trig, work] ;
ref-res // Chop // DeleteDuplicates
(* {0} *)

Timing:

length = 89128192 ;
signal = Sin[2*Pi*0.12345*Range[length]] ;
b = 1/Sqrt[2.0]  ;
t$ref = First@RepeatedTiming[Fourier[signal, FourierParameters -> {1, b}] ];
t$res = First@RepeatedTiming[{trig, work} = ffrft$set[length, b]  ; ffrft$get[length, signal, trig, work] ];
t$ref/t$res

For FourierParameters -> {1, 1} WM is most likely uses FFT from MKL, while for FourierParameters -> {a, b} it is not clear. While it is possible to commutecompute with n*log(n), I'm not sure it is the case in WM.

Here is an implementation based on Bailey's paper.

ClearAll[ffrft$set] ;
ffrft$set[length_, span_] := Block[
{trig, work},
trig =Exp[I*Pi*span/length*Range[0.0, length-1]^2] ;
work = Fourier[Flatten[{Conjugate[trig] , Exp[-I*Pi*span/length*Range[-length, -1.0]^2]}] , FourierParameters -> {1, 1} ];
{trig, work}
] ;

ClearAll[ffrft$get] ;
ffrft$get[length_, signal_, trig_, work_] := Block[
{ffrft},
ffrft = ConstantArray[0.0, 2*length]  ;
ffrft[[1;;length]] = trig*signal ;
ffrft = work*Fourier[ffrft, FourierParameters -> {1, 1} ] ;
ffrft =Fourier[ffrft, FourierParameters -> {1, -1} ][[1;;length]] ;
ffrft*trig/(2*length)
] ;

Example:

length = 1024 ;
signal = Sin[2*Pi*0.12345*Range[length]] ;
b = 1/Sqrt[2.0]  ;
{trig, work} = ffrft$set[length, b]  ;
ref = Fourier[signal, FourierParameters -> {1, b}] ;
res = ffrft$get[length, signal, trig, work] ;
ref-res // Chop // DeleteDuplicates
(* {0} *)

Timing:

length = 8912 ;
signal = Sin[2*Pi*0.12345*Range[length]] ;
b = 1/Sqrt[2.0]  ;
t$ref = First@RepeatedTiming[Fourier[signal, FourierParameters -> {1, b}] ];
t$res = First@RepeatedTiming[{trig, work} = ffrft$set[length, b]  ; ffrft$get[length, signal, trig, work] ];
t$ref/t$res

For FourierParameters -> {1, 1} WM is most likely uses FFT from MKL, while for FourierParameters -> {a, b} it is not clear. While it is possible to commute with n*log(n), I'm not sure it is the case in WM.

Here is an implementation based on Bailey's paper.

ClearAll[ffrft$set] ;
ffrft$set[length_, span_] := Block[
{trig, work},
trig =Exp[I*Pi*span/length*Range[0.0, length-1]^2] ;
work = Fourier[Flatten[{Conjugate[trig] , Exp[-I*Pi*span/length*Range[-length, -1.0]^2]}] , FourierParameters -> {1, 1} ];
{trig, work}
] ;

ClearAll[ffrft$get] ;
ffrft$get[length_, signal_, trig_, work_] := Block[
{ffrft},
ffrft = ConstantArray[0.0, 2*length]  ;
ffrft[[1;;length]] = trig*signal ;
ffrft = work*Fourier[ffrft, FourierParameters -> {1, 1} ] ;
ffrft =Fourier[ffrft, FourierParameters -> {1, -1} ][[1;;length]] ;
ffrft*trig/(2*length)
] ;

Example:

length = 1024 ;
signal = Sin[2*Pi*0.12345*Range[length]] ;
b = 1/Sqrt[2.0]  ;
{trig, work} = ffrft$set[length, b]  ;
ref = Fourier[signal, FourierParameters -> {1, b}] ;
res = ffrft$get[length, signal, trig, work] ;
ref-res // Chop // DeleteDuplicates
(* {0} *)

Timing:

length = 8192 ;
signal = Sin[2*Pi*0.12345*Range[length]] ;
b = 1/Sqrt[2.0]  ;
t$ref = First@RepeatedTiming[Fourier[signal, FourierParameters -> {1, b}] ];
t$res = First@RepeatedTiming[{trig, work} = ffrft$set[length, b]  ; ffrft$get[length, signal, trig, work] ];
t$ref/t$res

For FourierParameters -> {1, 1} WM is most likely uses FFT from MKL, while for FourierParameters -> {a, b} it is not clear. While it is possible to compute with n*log(n), I'm not sure it is the case in WM.

Source Link
I.M.
I.M.
  • 3.2k
  • 1
  • 14
  • 19

Here is an implementation based on Bailey's paper.

ClearAll[ffrft$set] ;
ffrft$set[length_, span_] := Block[
{trig, work},
trig =Exp[I*Pi*span/length*Range[0.0, length-1]^2] ;
work = Fourier[Flatten[{Conjugate[trig] , Exp[-I*Pi*span/length*Range[-length, -1.0]^2]}] , FourierParameters -> {1, 1} ];
{trig, work}
] ;

ClearAll[ffrft$get] ;
ffrft$get[length_, signal_, trig_, work_] := Block[
{ffrft},
ffrft = ConstantArray[0.0, 2*length]  ;
ffrft[[1;;length]] = trig*signal ;
ffrft = work*Fourier[ffrft, FourierParameters -> {1, 1} ] ;
ffrft =Fourier[ffrft, FourierParameters -> {1, -1} ][[1;;length]] ;
ffrft*trig/(2*length)
] ;

Example:

length = 1024 ;
signal = Sin[2*Pi*0.12345*Range[length]] ;
b = 1/Sqrt[2.0]  ;
{trig, work} = ffrft$set[length, b]  ;
ref = Fourier[signal, FourierParameters -> {1, b}] ;
res = ffrft$get[length, signal, trig, work] ;
ref-res // Chop // DeleteDuplicates
(* {0} *)

Timing:

length = 8912 ;
signal = Sin[2*Pi*0.12345*Range[length]] ;
b = 1/Sqrt[2.0]  ;
t$ref = First@RepeatedTiming[Fourier[signal, FourierParameters -> {1, b}] ];
t$res = First@RepeatedTiming[{trig, work} = ffrft$set[length, b]  ; ffrft$get[length, signal, trig, work] ];
t$ref/t$res

For FourierParameters -> {1, 1} WM is most likely uses FFT from MKL, while for FourierParameters -> {a, b} it is not clear. While it is possible to commute with n*log(n), I'm not sure it is the case in WM.