Improve interface of discrete cosine transforms

[cval26]

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.

[zoziha]

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

[cval26]

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.

[certik]

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?

[cval26]

@zoziha I suggest that we keep this issue open until the rest of the checklist items have been completed too.

[ivan-pi]

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: fortran-lang/benchmarks#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.