fortran-lang / fftpack Goto Github PK

Related to #9, #22.
EDIT: See PR #35 for the implementation of item 1.
EDIT: See PR #38 for the implementation of item 2.

I suggest the following changes, which I argue will bring the interface of the FFTPACK cosine transforms in line with modern terminology, and thus make the interface easier to understand and use. I am certainly willing to implement those changes myself, I just post them here first to gather feedback.

The most common discrete cosine transforms (DCT) are types 1, 2, 3 and 4. Up to a scaling factor, the FFTPACK procedures

• dcost and modern dctcorrespond to the DCT type-1 (DCT-I);
• dcosqb and modern iqct correspond to DCT type-2 (DCT-II);
• dcosqf and modern qctcorrespond to DCT type-3 (DCT-III).

Note also that DCT-I is the unnormalized inverse of itself; while DCT-III is the unnormalized inverse of DCT-II, and conversely. Wikipedia has a nice summary of the definitions of the most common types of DCT.

Because FFTPACK calls the DCT types 2 & 3 transforms quarter-wave transforms, it's difficult for a new user of the library (like me) to find out what they correspond to; the quarter terminology is not in use anymore. This is also why @zoziha in #9 couldn't find their corresponding transforms in scipy.

Based on the above, I propose the following changes.

• Combine the modern interfaces dct(x,n) and qct(x,n) under one interface dct(x,n,type) where type=1,2,3 (and similarly for the inverse ones). In this way, dct(x,n,type=1) corresponds to the DCT-I of the input data, and similarly for the others. This is the way the cosine transforms are also handled in scipy.fft. This will have no performance overhead, it's simply combining the two existing interfaces so that they're in line with current terminology.
• Create an interface to overload the existing in-place DCT subroutines, so that they can be called with names based on their type. For example, create a subroutine dct_t1 for DCT-I that simply overloads the existing dcost, and similarly for DCT types 2, 3. See comments below for code example; thanks to @zoziha for the suggestion.
• Write the definition formulas of the DCTs in LaTeX in the documentation.

These changes will allow us to potentially include more sine/cosine transforms in the future consistently, without having to break the interface later on. They will also allow us to improve the documentation to make it clear what each procedure computes.

I am willing to implement the above and change the documentation and tests accordingly.

We should also change the DST interface in the analogous way in a standalone PR for consistency.

This FFTPACK version seems to be released in the public domain. I think it is important to note in the README that contributions to this repository require to disclaim copyright or dedicate the contributed code to the public domain (if possible).

It looks like FFTPACK has "double precision" hard coded in. Would the maintainers be open to accepting a pull request that made this configurable? I would:

• change all instances of "double precision" to "real(rk)" and do the same for complex
• add a new file, that contains "integer, parameter :: rk = kind(1.0d0)"
• include that new file in all the old fortran77 files, so that the kind being used could be changed at compile time in exactly one place.

We could use a module instead of "include", but I'm not familiar enough with old-style loose functions to know if that would work.

The end result would be that the functions are still only provided for one kind, but rather than that kind being hard-coded to double precision, the user could choose the kind themselves. And just to be clear, I'm offering to do this, not requesting that someone else do it (though that's fine too, of course).

So what do we think? Does that involve too many changes to the old code? Or is it ok?

Do you think it would be useful to tag the initial version of this project? Maybe the use the first commit that works with fpm for this purpose. Even if it is only used a fixed reference point to compare against future version, for reasons of API and ABI compatibility.

You recommended me (see fortran-lang/fpm-registry#40 (comment)) to check whether the code in this repository(certik/fftpack) is consistent with netlib/dfftpack1.0 (fftpack4.0). I confirm that they are consistent.
And I want to learn from brocolis/fftpack to make dfftpack1.0(fftpack4.0) also support fpm. I will do this immediately.

Later, We update dfftpack1.0(fftpack4.0) to fftpack5.1:
netlib/dfftpack1.0(fftpack4.0) with a easy-to-use interface (-> ncar/fftpack5.1 with a easy-to-use interface -> john/fftpack5.1 with a easy-to-use interface) -> fftpack5.1 with a easy-to-use interface.

I got inspiration from stdlib, I think we can use fypp as a preprocessor to deal with the multi-precision problem of fftpack package.
Do we need to use fypp to refactor netlib/dffpack1.0(fftpack4.0) and nacr/fftpack5.1?

Tasks

• Let's collect the original sources from netlib. Commit. This might be similar to https://github.com/certik/fftpack, but we should check we got everything.
• we should add tests, the only test I see there is in https://github.com/certik/fftpack/blob/030cbf4fd35c35ef19387b4835bde3ac73019b64/test.f
• We should add Fortran 90 interface to the package (add fpm support, keep make support)
• Add 2D and 3D transforms
• Check any other functionality from the repositories you mentioned and see if there is anything else that we should include
  Such as fftshift, czt, hilbert, filter.
• Add benchmarks
• Some documentation how to use
• Multiple precision
• Use Fortran SAVE feature to save wsave(:)

Links

1. netlib/dfftpack1.0(fftpack4.0) (License : public domain)
   Fortran77
2. nacr/fftpack5.1 (License: look like MIT and BSD-3, see comment)
   Fortran77
3. John Burkardt/fftpack5.1 (License: LGPL/GPL)
   Fortran90 from nacr/fftpack5.1
4. QcmPlab/SciFortran (License: shall be GPL)
   From John Burkardt fortran90, fftpack5.1 with interface.
5. keurfonluu/FFTPack (License: shall be GPL)
   From John Burkardt fortran90, fftpack5.1 with interface.
6. jlokimlin/modern_fftpack (License: FFTPACK license look like MIT and BSD-3 ?)
   Fortran90, look like from nacr/fftpack5.1 ?

In the process of fftpack for example and testing, routines such as check, is_close, and disp are needed (for convenience). I personally think that it is possible to introduce stdlib-fpm and forlab dependencies only as dev-dependencies.

stdlib-fpm

I have submitted a PR (see fortran-lang/stdlib#508) in stdlib to support the stdlib-fpm branch through github-ci, and now I use LKedward/stdlib-fpm, I think I can switch to fortran-lang/stdlib-fpm in the future.
Since fortran-lang/stdlib-fpm is not yet supported, I introduced forlab. When stdlib-fpm is implemented, I will remove forlab.

hilbert transform

I implemented the hilbert transform function in my branch (see add_hilbert), which introduced stdlib-fpm and forlab as dev-dependencies.
Using the dev-dependencies feature of fpm means that it is not convenient for us to support the more make test of the make tool in the future.

See fortran-lang/fpm#421 for details.

Please collect references in @fortran-lang projects to files inside this repository such that we can fix those links after renaming.

Problem 1: The API name (i)qct?

I improved the interfaces and its documentation of dcosq routines in this branch (I may consider submitting a PR tomorrow?). They compute the fast fourier transform of quarter wave data, but I don't know enough about them ("quarter wave data"), this involves more professional terms and I can't even find their corresponding routines in scipy.fftpack.
What should their simplified modern interface be called? qct(Quarter-Cosine Transform)?
(see https://www2.cisl.ucar.edu/resources/legacy/fft5/documentation)

Problem 2: more complex document (fftpack.md) that is difficult to read for users?

fftpack.md is getting more complicated and may be difficult to read?
It contains original fftpack routines (dfft/dzfft/...), modern interfaces (fft/rfft/...), and utility functions (fftshift).

Maybe we can divide the existing fftpack.md into three markdown documents: fftpack_classic_API.md, fftpack_modern_API.md, fftpack_utilities.md.

This issue is not urgent, I am just worried.

Propose creating and adding a FFTPACK logo.

Could someone help me understand the storage sequence of the rfft function result and how to get the actual complex result from the forward FFT of a real data? I wrote the test program below to transform a 3 + cos(x). I'm unclear why the rfft result is real so I tried using transfer to get the corresponding complex values. The result is close to what I'd expect but not quite.

program main
  use fftpack, only: rfft, irfft
  implicit none
  integer m
  integer, parameter :: N = 8
  double precision, parameter :: two_pi = 2.D0*acos(-1.D0), tolerance = 1.0D-06, x_avg = 3.D0
  double precision, parameter :: x(*) = x_avg + [(cos(two_pi*dble(m)/dble(N)), m=0,N-1)]
  double precision, allocatable, dimension(:) :: rfft_x, x_round_trip
  double complex, allocatable :: x_hat(:)

  rfft_x = rfft(x)/dble(N)
  if (size(rfft_x) /= size(x)) error stop "rfft_x/x size mismatch"

  x_round_trip = real(irfft(rfft_x))

  if (size(x_round_trip) /= size(x)) error stop "x_round_trip/x: size mismatch"
  if (any(abs(x_round_trip - x) > tolerance)) error stop "round trip out of tolerance"

  print '(3(10x,a,10x))',"x", "x_round_trip", "rfft_x"
  do m = 1, size(x)
    print *, x(m), x_round_trip(m), rfft_x(m)
  end do
  print *

  x_hat = transfer(rfft_x, mold=x_hat, size=size(x)/2)
  print *, "x_hat = ", x_hat
end program

which gives the output

          x                    x_round_trip                    rfft_x
   4.0000000000000000        4.0000000000000000        3.0000000000000000     
   3.7071067811865475        3.7071068286895752       0.50000000000000000     
   3.0000000000000000        3.0000000000000000       -2.7755575615628914E-017
   2.2928932188134525        2.2928931713104248        0.0000000000000000     
   2.0000000000000000        2.0000000000000000       -0.0000000000000000     
   2.2928932188134521        2.2928931713104248        0.0000000000000000     
   3.0000000000000000        3.0000000000000000       -2.7755575615628914E-017
   3.7071067811865475        3.7071068286895752        0.0000000000000000     

 x_hat =               (3.0000000000000000,0.50000000000000000)        (-2.77555756156289135E-017,0.0000000000000000)              (-0.0000000000000000,0.0000000000000000)        (-2.77555756156289135E-017,0.0000000000000000)

where I would have expected like

 x_hat = (3., 0.) (0.5, 0.) (0.,0.) (0.,0.)

and a corresponding reordering of rfft_x.