fftpack
fftpack copied to clipboard
Improve interface of discrete cosine transforms
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 moderndct
correspond to the DCT type-1 (DCT-I); -
dcosqb
and moderniqct
correspond to DCT type-2 (DCT-II); -
dcosqf
and modernqct
correspond 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.
- [x] Combine the modern interfaces
dct(x,n)
andqct(x,n)
under one interfacedct(x,n,type)
wheretype=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 inscipy.fft
. This will have no performance overhead, it's simply combining the two existing interfaces so that they're in line with current terminology. - [x] 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 existingdcost
, 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.
Thank you very much for your investigation @cval26 , I understand your content!
Combining to create a unified dct
and idct
interface is a great decision.
But I'm worried that changing the old routine names (such as dcost
, dcosqb
) might require some discussion, @certik ? However, I know that Fortran can easily overload function names to preserve the original interface, such as:
module fftpack
interface dct
module procedure :: dcost ! or just `procedure :: dcost` for external/global routine.
end interface dct
contains
subroutine dcost()
print *, "hello from subroutine dcost"
end subroutine dcost
end module fftpack
program main
use fftpack
call dct()
call dcost()
end program main
Thank you for the illustration of overloading procedure names, @zoziha. I think this is actually a better way of providing a modern interface for the in-place cosine transforms, while also preserving the classic interface. It also allows us to reuse the existing code, without duplicating code to achieve the result.
For example, if we were to include an interface for DCT types 1, 2, 3, we could do so in the following way, if I understand you correctly.
module fftpack
interface dct_t1 ! DCT type-1
module procedure :: dcost
end interface dct_t1
interface dct_t2 ! DCT type-2
module procedure :: dcosqb
end interface dct_t2
interface dct_t3 ! DCT type-3
module procedure :: dcosqf
end interface dct_t3
contains
! where dcost, dcosqb, dcosqf are the existing, "classic" subroutines
end module fftpack
And similarly, when we extend, say, dcost
to 2D dcost2
, we will be able to create an interface for a 2D DCT-I in the same way.
I believe this is an excellent suggestion. I'm going to edit my initial comment to propose this as the suggested change, instead of renaming the original subroutines. I'll wait for feedback from others too, to see if they agree.
Also, just to explain my motivation for providing an interface to the in-place cosine transforms with names based on their type, I often use the in-place subroutines (say, dcost
) to avoid the dynamic memory allocations of the modern wrapper dct
. But because of the lack of clear documentation and the fact that their names are not based on their transform type, it wasn't initially clear to me which transform does what. By now I know what type they correspond to, but we should make it easier for new users to find the transform they want, and be able to take advantage of the performance of the in-place subroutines. I believe the above interface suggestion achieves that in the most straightforward way, in addition to improving the documentation for the DCTs.
It seems like a good idea, especially if we preserve (at least for now) the old names as indicated by @zoziha.
@jacobwilliams, what do you think?
@zoziha I suggest that we keep this issue open until the rest of the checklist items have been completed too.
Sorry, it was the PR merger that caused the closure of this issue, which has now been reopened.
In case anyone is interested - Apple has a DCT method for the default float type in their Accelerate framework: https://developer.apple.com/documentation/accelerate/discrete_cosine_transforms?language=objc
The interface in Fortran would look something like this:
enum, bind(c)
enumerator :: DCT_II = 2
enumerator :: DCT_III = 3
enumerator :: DCT_IV = 4
end enum
interface
function DCT_CreateSetup(previous,length,type) &
bind(c,name="vDSP_DCT_CreateSetup")
type(c_ptr), value, optional :: previous
integer(c_long), value :: length
integer(c_int), value :: type
end function
subroutine DCT_Execute(setup,input,ouput) &
bind(c,name="vDSP_DCT_Execute")
type(c_ptr), value :: setup
real(c_float), intent(in) :: input(*)
real(c_float), intent(out) :: output(*)
end subroutine
subroutine DCT_DestroySetup(setup) &
bind(c,name="vDSP_DFT_DestroySetup")
type(c_ptr), value :: setup
end subroutine
end interface
The DCT_Execute
routine can be used in-place (a Fortran wrapper with a single dummy argument would be needed to map the intent
attribute correctly). There are some restrictions on the sizes.
It could be nice to compare performance of the FFTPACK DCT, versus the one in Accelerate. This could go to the Fortran-Lang benchmarks repository, see: https://github.com/fortran-lang/benchmarks/issues/18
A usage example would look something like:
integer :: n
real(c_float), allocatable :: x(:), y(:)
type(c_ptr) :: dct_setup
n = ...
allocate(x(n),y(n))
x = ...
dct_setup = DCT_CreateSetup(length=n,type=DCT_IV)
call DCT_Execute(dct_setup,input=x,output=y)
call DCT_DestroySetup(dct_setup)
To link it, add -framework Accelerate
to the command.