Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add type-based DCT procedure names #38

Merged
merged 5 commits into from
Mar 23, 2023
Merged
Show file tree
Hide file tree
Changes from 3 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 22 additions & 7 deletions doc/specs/fftpack.md
Original file line number Diff line number Diff line change
Expand Up @@ -809,13 +809,16 @@ end program demo_dzfftb

## DCT type-1 (DCT-1)

### Initialize DCT-1: `dcosti`
### Initialize DCT-1: `dcosti` or `dct_t1i`

#### Description

Initializes the array `wsave` which is used in subroutine `dcost`.
The prime factorization of `n` together with a tabulation of the trigonometric functions are computed and stored in `wsave`.

The two procedures are completely equivalent and expect the same arguments.
It is a matter of personal preference which one you choose to use.

#### Status

Experimental
Expand Down Expand Up @@ -851,18 +854,21 @@ program demo_dcosti
end program demo_dcosti
```

### Compute DCT-1: `dcost`
### Compute DCT-1: `dcost` or `dct_t1`

#### Description

Computes the DCT-1 of the input real data.
The transform is defined below at output parameter `x`.

The two procedures are completely equivalent and expect the same arguments.
It is a matter of personal preference which one you choose to use.

For real input data `x` of length `n`, the DCT-1 of `x` is equivalent, up to a
scaling factor, to the DFT of the even extension of `x` with length `2*(n-1)`,
where the first and last entries of the original data are not repeated in the
extension. For example, the DCT-1 of input data *abcde* (size \[5\]) is
equivalent to the DFT of data *abcdedcb* (size \[2*4=8\]).
extension. For example, the DCT-1 of input data *abcde* (size \(5\)) is
equivalent to the DFT of data *abcdedcb* (size \(2*4=8\)).

Also, `dcost` is the unnormalized inverse of itself. This means that a call of
`dcost` followed by another call of `dcost` will multiply the input sequence `x`
Expand Down Expand Up @@ -932,7 +938,7 @@ end program demo_dcost

## DCT of types 2, 3 (DCT-2, 3), a.k.a "Quarter" cosine transforms

### Initialize DCT-2, 3: `dcosqi`
### Initialize DCT-2, 3: `dcosqi` or `dct_t23i`

#### Description

Expand All @@ -941,6 +947,9 @@ The prime factorization of `n` together with
a tabulation of the trigonometric functions are computed and
stored in `wsave`.

The two procedures are completely equivalent and expect the same arguments.
It is a matter of personal preference which one you choose to use.

#### Status

Experimental
Expand Down Expand Up @@ -978,13 +987,16 @@ program demo_dcosqi
end program demo_dcosqi
```

### Compute DCT-3: `dcosqf`
### Compute DCT-3: `dcosqf` or `dct_t3`

#### Description

Computes the DCT-3 of the input real data.
The transform is defined below at output parameter `x`.

The two procedures are completely equivalent and expect the same arguments.
It is a matter of personal preference which one you choose to use.

Also, `dcosqf` (DCT-3) is the unnormalized inverse of `dcosqb` (DCT-2), since a
call of `dcosqf` followed by a call of `dcosqb` will multiply the input sequence
`x` by `4*n`.
Expand Down Expand Up @@ -1049,13 +1061,16 @@ program demo_dcosqf
end program demo_dcosqf
```

### Compute DCT-2: `dcosqb`
### Compute DCT-2: `dcosqb` or `dct_t2`

#### Description

Computes the DCT-2 of the input real data.
The transform is defined below at output parameter `x`.

The two procedures are completely equivalent and expect the same arguments.
It is a matter of personal preference which one you choose to use.

For real input data `x` of length `n`, the DCT-2 of `x` is equivalent, up to a
scaling factor, to the DFT of the even extension of `x` with length `4*n`,
where all the even-frequency entries are zero.
Expand Down
42 changes: 42 additions & 0 deletions src/fftpack.f90
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,8 @@ module fftpack
public :: dcosqi, dcosqf, dcosqb
public :: dcosti, dcost
public :: dct, idct
public :: dct_t1i, dct_t1
public :: dct_t23i, dct_t2, dct_t3

public :: rk

Expand Down Expand Up @@ -268,6 +270,46 @@ pure module function idct_rk(x, n, type) result(result)
end function idct_rk
end interface idct

!> Version: experimental
!>
!> Initialize DCT type-1
!> ([Specification](../page/specs/fftpack.html#dct_t1i))
cval26 marked this conversation as resolved.
Show resolved Hide resolved
interface dct_t1i
procedure :: dcosti
end interface dct_t1i

!> Version: experimental
!>
!> Perform DCT type-1
!> ([Specification](../page/specs/fftpack.html#dct_t1))
interface dct_t1
procedure :: dcost
end interface dct_t1

!> Version: experimental
!>
!> Initialize DCT types 2, 3
!> ([Specification](../page/specs/fftpack.html#dct_t23i))
interface dct_t23i
procedure :: dcosqi
end interface dct_t23i

!> Version: experimental
!>
!> Perform DCT type-2
!> ([Specification](../page/specs/fftpack.html#dct_t2))
interface dct_t2
procedure :: dcosqb
end interface dct_t2

!> Version: experimental
!>
!> Perform DCT type-3
!> ([Specification](../page/specs/fftpack.html#dct_t3))
interface dct_t3
procedure :: dcosqf
end interface dct_t3

!> Version: experimental
!>
!> Shifts zero-frequency component to center of spectrum.
Expand Down
36 changes: 31 additions & 5 deletions test/test_fftpack_dct.f90
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
module test_fftpack_dct

use fftpack, only: rk, dcosti, dcost, dct, idct, dcosqi, dcosqf, dcosqb
use fftpack
use testdrive, only: new_unittest, unittest_type, error_type, check
implicit none
private
Expand All @@ -24,27 +24,53 @@ end subroutine collect_dct

subroutine test_classic_dct(error)
type(error_type), allocatable, intent(out) :: error
real(kind=rk) :: w(3*4 + 15)
real(kind=rk) :: w(3*4 + 15), w2(3*4 + 15)
real(kind=rk) :: x(4) = [1, 2, 3, 4]
real(kind=rk) :: x2(4)
real(kind=rk) :: eps = 1.0e-10_rk


x2 = x
call dcosti(4, w)
call dcost(4, x, w)
call check(error, all(x == [real(kind=rk) :: 15, -4, 0, -1.0000000000000009_rk]), "`dcosti` failed.")
call check(error, sum(abs(x - [real(kind=rk) :: 15, -4, 0, -1.0000000000000009_rk])) < eps, &
"`dcost` failed.")
if (allocated(error)) return

call dct_t1i(4, w2)
call dct_t1(4, x2, w2)
call check(error, maxval(abs(x2-x)) < eps, "dct_t1 failed")
if (allocated(error)) return

call dcost(4, x, w)
call check(error, all(x/(2*3) == [real(kind=rk) :: 1, 2, 3, 4]), "`dcost` failed.")
call check(error, sum(abs(x/(2*3) - [real(kind=rk) :: 1, 2, 3, 4])) < eps, &
"2nd `dcost` failed.")
if (allocated(error)) return

call dct_t1(4, x2, w2)
call check(error, maxval(abs(x2-x)) < eps, "2nd dct_t1 failed")
if (allocated(error)) return

x = [1, 2, 3, 4]
x2 = x
call dcosqi(4, w)
call dcosqf(4, x, w)
call check(error, sum(abs(x - [11.999626276085150_rk, -9.1029432177492193_rk, &
2.6176618435106480_rk, -1.5143449018465791_rk])) < eps, &
"`dcosqf` failed.")
if (allocated(error)) return

call dct_t23i(4, w2)
call dct_t3(4, x2, w2)
call check(error, maxval(abs(x2-x)) < eps, "dct_t3 failed")
if (allocated(error)) return

call dcosqb(4, x, w)
call check(error, sum(abs(x/(4*4) - [real(kind=rk) :: 1, 2, 3, 4])) < eps, &
"`dcosqb` failed.")
if (allocated(error)) return

call dct_t2(4, x2, w2)
call check(error, maxval(abs(x2-x)) < eps, "dct_t2 failed")

end subroutine test_classic_dct

Expand Down