From 13a72e4d55708d6c89f106045dcf9943c81a273e Mon Sep 17 00:00:00 2001 From: Chris Val Date: Thu, 2 Mar 2023 12:25:34 -0500 Subject: [PATCH 1/5] combine dct/qct interfaces --- src/CMakeLists.txt | 2 -- src/Makefile | 4 --- src/fftpack.f90 | 40 ++++++--------------- src/fftpack_dct.f90 | 85 +++++++++++++++++++++++++++++++++++++++----- src/fftpack_iqct.f90 | 36 ------------------- src/fftpack_qct.f90 | 36 ------------------- 6 files changed, 86 insertions(+), 117 deletions(-) delete mode 100644 src/fftpack_iqct.f90 delete mode 100644 src/fftpack_qct.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f27225a..a58de7d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -31,9 +31,7 @@ set(FFTPACK_SOURCES ${dir}/fftpack_fftshift.f90 ${dir}/fftpack_ifft.f90 ${dir}/fftpack_ifftshift.f90 - ${dir}/fftpack_iqct.f90 ${dir}/fftpack_irfft.f90 - ${dir}/fftpack_qct.f90 ${dir}/fftpack_rfft.f90 ${dir}/fftpack_utils.f90 ${dir}/passb.f90 diff --git a/src/Makefile b/src/Makefile index 3ad8ac0..884a6f1 100644 --- a/src/Makefile +++ b/src/Makefile @@ -57,8 +57,6 @@ SRCF90 = \ fftpack_irfft.f90\ fftpack_fftshift.f90\ fftpack_ifftshift.f90\ - fftpack_qct.f90\ - fftpack_iqct.f90\ fftpack_dct.f90\ rk.f90\ fftpack_utils.f90 @@ -82,8 +80,6 @@ fftpack_fft.o: fftpack.o rk.o fftpack_ifft.o: fftpack.o rk.o fftpack_rfft.o: fftpack.o rk.o fftpack_irfft.o: fftpack.o rk.o -fftpack_qct.o: fftpack.o rk.o -fftpack_iqct.o: fftpack.o rk.o fftpack_dct.o: fftpack.o rk.o fftpack_fftshift.o: fftpack.o rk.o fftpack_ifftshift.o: fftpack.o rk.o diff --git a/src/fftpack.f90 b/src/fftpack.f90 index d69283f..3fff9d6 100644 --- a/src/fftpack.f90 +++ b/src/fftpack.f90 @@ -15,8 +15,6 @@ module fftpack public :: dzffti, dzfftf, dzfftb public :: dcosqi, dcosqf, dcosqb - public :: qct, iqct - public :: dcosti, dcost public :: dct, idct @@ -246,46 +244,28 @@ end function irfft_rk !> Version: experimental !> - !> Forward transform of quarter wave data. - !> ([Specifiction](../page/specs/fftpack.html#qct)) - interface qct - pure module function qct_rk(x, n) result(result) - real(kind=rk), intent(in) :: x(:) - integer, intent(in), optional :: n - real(kind=rk), allocatable :: result(:) - end function qct_rk - end interface qct - - !> Version: experimental - !> - !> Backward transform of quarter wave data. - !> ([Specifiction](../page/specs/fftpack.html#iqct)) - interface iqct - pure module function iqct_rk(x, n) result(result) - real(kind=rk), intent(in) :: x(:) - integer, intent(in), optional :: n - real(kind=rk), allocatable :: result(:) - end function iqct_rk - end interface iqct - - !> Version: experimental - !> - !> Discrete fourier cosine (forward) transform of an even sequence. + !> Dsicrete cosine transforms. !> ([Specification](../page/specs/fftpack.html#dct)) interface dct - pure module function dct_rk(x, n) result(result) + pure module function dct_rk(x, n, t) result(result) real(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n + integer, intent(in), optional :: t real(kind=rk), allocatable :: result(:) end function dct_rk end interface dct !> Version: experimental !> - !> Discrete fourier cosine (backward) transform of an even sequence. + !> Inverse discrete cosine transforms. !> ([Specification](../page/specs/fftpack.html#idct)) interface idct - module procedure :: dct_rk + pure module function idct_rk(x, n, t) result(result) + real(kind=rk), intent(in) :: x(:) + integer, intent(in), optional :: n + integer, intent(in), optional :: t + real(kind=rk), allocatable :: result(:) + end function idct_rk end interface idct !> Version: experimental diff --git a/src/fftpack_dct.f90 b/src/fftpack_dct.f90 index 88a46d9..e8569ca 100644 --- a/src/fftpack_dct.f90 +++ b/src/fftpack_dct.f90 @@ -2,10 +2,11 @@ contains - !> Discrete fourier cosine transform of an even sequence. - pure module function dct_rk(x, n) result(result) + !> Discrete cosine transforms of types 1, 2, 3. + pure module function dct_rk(x, n, t) result(result) real(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n + integer, intent(in), optional :: t real(kind=rk), allocatable :: result(:) integer :: lenseq, lensav, i @@ -23,14 +24,80 @@ pure module function dct_rk(x, n) result(result) result = x end if - !> Initialize FFT - lensav = 3*lenseq + 15 - allocate (wsave(lensav)) - call dcosti(lenseq, wsave) - - !> Discrete fourier cosine transformation - call dcost(lenseq, result, wsave) + ! Default to DCT-2 + if (.not.present(t)) then + lensav = 3*lenseq + 15 + allocate (wsave(lensav)) + call dcosqi(lenseq, wsave) + call dcosqb(lenseq, result, wsave) + return + end if + if (t == 1) then ! DCT-1 + lensav = 3*lenseq + 15 + allocate (wsave(lensav)) + call dcosti(lenseq, wsave) + call dcost(lenseq, result, wsave) + else if (t == 2) then ! DCT-2 + lensav = 3*lenseq + 15 + allocate (wsave(lensav)) + call dcosqi(lenseq, wsave) + call dcosqb(lenseq, result, wsave) + else if (t == 3) then ! DCT-3 + lensav = 3*lenseq + 15 + allocate (wsave(lensav)) + call dcosqi(lenseq, wsave) + call dcosqf(lenseq, result, wsave) + end if end function dct_rk + !> Inverse discrete cosine transforms of types 1, 2, 3. + pure module function idct_rk(x, n, t) result(result) + real(kind=rk), intent(in) :: x(:) + integer, intent(in), optional :: n + integer, intent(in), optional :: t + real(kind=rk), allocatable :: result(:) + + integer :: lenseq, lensav, i + real(kind=rk), allocatable :: wsave(:) + + if (present(n)) then + lenseq = n + if (lenseq <= size(x)) then + result = x(:lenseq) + else if (lenseq > size(x)) then + result = [x, (0.0_rk, i=1, lenseq - size(x))] + end if + else + lenseq = size(x) + result = x + end if + + ! Default to t=2; inverse DCT-2 is DCT-3 + if (.not.present(t)) then + lensav = 3*lenseq + 15 + allocate (wsave(lensav)) + call dcosqi(lenseq, wsave) + call dcosqf(lenseq, result, wsave) + return + end if + + if (t == 1) then ! inverse DCT-1 is DCT-1 + lensav = 3*lenseq + 15 + allocate (wsave(lensav)) + call dcosti(lenseq, wsave) + call dcost(lenseq, result, wsave) + else if (t == 2) then ! inverse DCT-2 is DCT-3 + lensav = 3*lenseq + 15 + allocate (wsave(lensav)) + call dcosqi(lenseq, wsave) + call dcosqf(lenseq, result, wsave) + else if (t == 3) then ! inverse DCT-3 is DCT-2 + lensav = 3*lenseq + 15 + allocate (wsave(lensav)) + call dcosqi(lenseq, wsave) + call dcosqb(lenseq, result, wsave) + end if + end function idct_rk + end submodule fftpack_dct diff --git a/src/fftpack_iqct.f90 b/src/fftpack_iqct.f90 deleted file mode 100644 index 1ae2a20..0000000 --- a/src/fftpack_iqct.f90 +++ /dev/null @@ -1,36 +0,0 @@ -submodule(fftpack) fftpack_iqct - -contains - - !> Backward transform of quarter wave data. - pure module function iqct_rk(x, n) result(result) - real(kind=rk), intent(in) :: x(:) - integer, intent(in), optional :: n - real(kind=rk), allocatable :: result(:) - - integer :: lenseq, lensav, i - real(kind=rk), allocatable :: wsave(:) - - if (present(n)) then - lenseq = n - if (lenseq <= size(x)) then - result = x(:lenseq) - else if (lenseq > size(x)) then - result = [x, (0.0_rk, i=1, lenseq - size(x))] - end if - else - lenseq = size(x) - result = x - end if - - !> Initialize FFT - lensav = 3*lenseq + 15 - allocate (wsave(lensav)) - call dcosqi(lenseq, wsave) - - !> Backward transformation - call dcosqb(lenseq, result, wsave) - - end function iqct_rk - -end submodule fftpack_iqct diff --git a/src/fftpack_qct.f90 b/src/fftpack_qct.f90 deleted file mode 100644 index ceb2f6b..0000000 --- a/src/fftpack_qct.f90 +++ /dev/null @@ -1,36 +0,0 @@ -submodule(fftpack) fftpack_qct - -contains - - !> Forward transform of quarter wave data. - pure module function qct_rk(x, n) result(result) - real(kind=rk), intent(in) :: x(:) - integer, intent(in), optional :: n - real(kind=rk), allocatable :: result(:) - - integer :: lenseq, lensav, i - real(kind=rk), allocatable :: wsave(:) - - if (present(n)) then - lenseq = n - if (lenseq <= size(x)) then - result = x(:lenseq) - else if (lenseq > size(x)) then - result = [x, (0.0_rk, i=1, lenseq - size(x))] - end if - else - lenseq = size(x) - result = x - end if - - !> Initialize FFT - lensav = 3*lenseq + 15 - allocate (wsave(lensav)) - call dcosqi(lenseq, wsave) - - !> Forward transformation - call dcosqf(lenseq, result, wsave) - - end function qct_rk - -end submodule fftpack_qct From a7a74d9c868a39097c39dd885edb72b75d0a4ad9 Mon Sep 17 00:00:00 2001 From: Chris Val Date: Thu, 2 Mar 2023 12:26:45 -0500 Subject: [PATCH 2/5] combine & add dct tests --- test/CMakeLists.txt | 1 - test/Makefile | 9 ++--- test/test_fftpack.f90 | 2 - test/test_fftpack_dct.f90 | 81 ++++++++++++++++++++++++++++++++------- test/test_fftpack_qct.f90 | 77 ------------------------------------- 5 files changed, 71 insertions(+), 99 deletions(-) delete mode 100644 test/test_fftpack_qct.f90 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 471697b..4f36b69 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -14,7 +14,6 @@ endmacro() set(FFTPACK_TEST_SOURCES test_fftpack_dct.f90 test_fftpack_fft.f90 - test_fftpack_qct.f90 test_fftpack_rfft.f90 test_fftpack_utils.f90 test_fftpack.f90 diff --git a/test/Makefile b/test/Makefile index cabb8a3..5429e23 100644 --- a/test/Makefile +++ b/test/Makefile @@ -3,12 +3,11 @@ FETCH = curl -L SRC = \ test_fftpack_fft.f90 \ test_fftpack_rfft.f90 \ - test_fftpack_qct.f90 \ test_fftpack_dct.f90 \ test_fftpack_utils.f90 \ test_fftpack.f90 \ testdrive.F90 - + OBJ = $(SRC:.f90=.o) OBJ := $(OBJ:.F90=.o) @@ -24,10 +23,10 @@ tstfft: tstfft.f test_fftpack: $(OBJ) $(FC) $(FFLAGS) $(OBJ) -L../src -l$(LIB) -I../src -o $@.x ./test_fftpack.x - + testdrive.F90: $(FETCH) https://github.com/fortran-lang/test-drive/raw/v0.4.0/src/testdrive.F90 > $@ - + %.o: %.F90 $(FC) $(FFLAGS) -c $< @@ -36,14 +35,12 @@ testdrive.F90: test_fftpack.o: test_fftpack_fft.o \ test_fftpack_rfft.o \ - test_fftpack_qct.o \ test_fftpack_dct.o \ test_fftpack_utils.o \ testdrive.o test_fftpack_fft.o: testdrive.o test_fftpack_rfft.o: testdrive.o -test_fftpack_qct.o: testdrive.o test_fftpack_dct.o: testdrive.o test_fftpack_utils.o: testdrive.o diff --git a/test/test_fftpack.f90 b/test/test_fftpack.f90 index 3396421..b3b4f47 100644 --- a/test/test_fftpack.f90 +++ b/test/test_fftpack.f90 @@ -3,7 +3,6 @@ program test_fftpack use testdrive, only: run_testsuite, new_testsuite, testsuite_type use test_fftpack_fft, only: collect_fft use test_fftpack_rfft, only: collect_rfft - use test_fftpack_qct, only: collect_qct use test_fftpack_dct, only: collect_dct use test_fftpack_utils, only: collect_utils implicit none @@ -16,7 +15,6 @@ program test_fftpack testsuites = [ & new_testsuite("fft", collect_fft), & new_testsuite("rfft", collect_rfft), & - new_testsuite("qct", collect_qct), & new_testsuite("dct", collect_dct), & new_testsuite("utils", collect_utils) & ] diff --git a/test/test_fftpack_dct.f90 b/test/test_fftpack_dct.f90 index e435923..0256cd5 100644 --- a/test/test_fftpack_dct.f90 +++ b/test/test_fftpack_dct.f90 @@ -1,6 +1,6 @@ module test_fftpack_dct - use fftpack, only: rk, dcosti, dcost, dct, idct + use fftpack, only: rk, dcosti, dcost, dct, idct, dcosqi, dcosqf, dcosqb use testdrive, only: new_unittest, unittest_type, error_type, check implicit none private @@ -32,21 +32,54 @@ subroutine test_classic_dct(error) call dcost(4, x, w) call check(error, all(x == [real(kind=rk) :: 15, -4, 0, -1.0000000000000009_rk]), "`dcosti` failed.") if (allocated(error)) return - call dcost(4, x, w) - call check(error, all(x/(2.0_rk*(4.0_rk - 1.0_rk)) == [real(kind=rk) :: 1, 2, 3, 4]), "`dcost` failed.") + call check(error, all(x/(2*3) == [real(kind=rk) :: 1, 2, 3, 4]), "`dcost` failed.") + + x = [1, 2, 3, 4] + 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 dcosqb(4, x, w) + call check(error, sum(abs(x/(4*4) - [real(kind=rk) :: 1, 2, 3, 4])) < eps, & + "`dcosqb` failed.") end subroutine test_classic_dct subroutine test_modernized_dct(error) type(error_type), allocatable, intent(out) :: error + real(kind=rk) :: eps = 1.0e-10_rk real(kind=rk) :: x(3) = [9, -9, 3] - call check(error, all(dct(x, 2) == [real(kind=rk) :: 0, 18]), "`dct(x, 2)` failed.") + ! DCT-1 + call check(error, sum(abs(dct(x,2,1) - [0.0_rk, 18.0_rk])) < eps, "`dct(x,2,1)` failed.") + if (allocated(error)) return + call check(error, sum(abs(dct(x,3,1) - dct(x,t=1))) < eps, "`dct(x,3,1)` failed.") + if (allocated(error)) return + call check(error, sum(abs(dct(x,4,1) - [real(kind=rk) :: -3, -3.0000000000000036_rk, 15, 33])) < eps,& + "`dct(x,4,1)` failed.") + !DCT-2 + x = [9, -9, 3] + call check(error, sum(abs(dct(x,3,2) - [12.0_rk, 20.784609690826525_rk, 60.0_rk])) < eps,& + "`dct(x,3,2)` failed.") + call check(error, sum(abs(dct(x,3) - [12.0_rk, 20.784609690826525_rk, 60.0_rk])) < eps,& + "`dct(x,3)` failed.") + call check(error, sum(abs(dct(x,4,2) - [12.0_rk, 14.890858416882008_rk, 42.426406871192853_rk,& + 58.122821125684993_rk])) < eps, "`dct(x,4,2)` failed.") + + ! DCT-3 + x = [9, -9, 3] + call check(error, sum(abs(dct(x,2,3) - [-3.7279220613578570_rk, 21.727922061357859_rk])) < eps, & + "`dct(x,2,3)` failed.") if (allocated(error)) return - call check(error, all(dct(x, 3) == dct(x)), "`dct(x, 3)` failed.") + call check(error, sum(abs(dct(x,3,3) - dct(x,t=3))) < eps, & + "`dct(x,3,3)` failed.") if (allocated(error)) return - call check(error, all(dct(x, 4) == [real(kind=rk) :: -3, -3.0000000000000036_rk, 15, 33]), "`dct(x, 4)` failed.") + call check(error, sum(abs(dct(x,4,3) - [-3.3871908980838743_rk, -2.1309424696909023_rk, & + 11.645661095452331_rk, 29.872472272322447_rk])) < eps, & + "`dct(x,n=4,t=3)` failed.") end subroutine test_modernized_dct @@ -55,15 +88,37 @@ subroutine test_modernized_idct(error) real(kind=rk) :: eps = 1.0e-10_rk real(kind=rk) :: x(4) = [1, 2, 3, 4] - call check(error, all(idct(dct(x))/(2.0_rk*(4.0_rk - 1.0_rk)) == [real(kind=rk) :: 1, 2, 3, 4]), & - "`idct(dct(x))/(2.0_rk*(4.0_rk-1.0_rk))` failed.") + ! inverse DCT-1 + call check(error, sum(abs(idct(dct(x,t=1),t=1)/(2*3) - x)) < eps, "`idct(dct(x,t=1),t=1)/(2*3)` failed.") + if (allocated(error)) return + call check(error, sum(abs(idct(dct(x,t=1), 2, 1)/(2*1) - [5.5_rk, 9.5_rk])) < eps,& + "`idct(dct(x,t=1), 2, 1)/(2*1)` failed.") + if (allocated(error)) return + call check(error, sum(abs(idct(dct(x,2,1), 4, 1)/(2*3) - [0.16666666666666666_rk, 0.33333333333333331_rk,& + 0.66666666666666663_rk, 0.83333333333333315_rk])) < eps,& + "`idct(dct(x,2,1), 4, 1)/(2*3)` failed.") + + ! inverse DCT-2 + x = [1, 2, 3, 4] + call check(error, sum(abs(idct(dct(x,t=2))/(4*4) - x)) < eps, & + "`idct(dct(x, t=2))/(4*4)` failed.") + if (allocated(error)) return + call check(error, sum(abs(idct(dct(x,t=2),n=2) - [22.156460020898692_rk, 57.843539979101308_rk])) < eps,& + "`idct(dct(x,t=2),n=2)` failed.") + if (allocated(error)) return + call check(error, sum(abs(idct(dct(x,n=2,t=2),n=4) - [6.7737481404944937_rk, 9.8352155994152106_rk,& + 14.164784400584789_rk, 17.226251859505506_rk])) < eps, "`idct(dct(x,n=2,t=2),n=4)` failed.") + + ! inverse DCT-3 + x = [1, 2, 3, 4] + call check(error, sum(abs(idct(dct(x,t=3),t=3)/(4*4) - x)) < eps, & + "`idct(dct(x, t=3), t=3)/(4*4)` failed.") if (allocated(error)) return - call check(error, all(idct(dct(x), 2)/(2.0_rk*(2.0_rk - 1.0_rk)) == [real(kind=rk) :: 5.5, 9.5]), & - "`idct(dct(x), 2)/(2.0_rk*(2.0_rk-1.0_rk))` failed.") + call check(error, sum(abs(idct(dct(x,t=3),n=2,t=3)/(4*2) - [1.4483415291679655_rk, 7.4608849947753271_rk])) < eps, & + "`idct(dct(x,t=3),n=2,t=3)/(4*2)` failed.") if (allocated(error)) return - call check(error, all(idct(dct(x, 2), 4)/(2.0_rk*(4.0_rk - 1.0_rk)) == & - [0.16666666666666666_rk, 0.33333333333333331_rk, 0.66666666666666663_rk, 0.83333333333333315_rk]), & - "`idct(dct(x, 2), 4)/(2.0_rk*(4.0_rk-1.0_rk))` failed.") + call check(error, sum(abs(idct(dct(x,n=2,t=3),n=4,t=3)/(4*4) - [0.5_rk, 0.70932417358418376_rk, 1.0_rk, & + 0.78858050747473762_rk])) < eps, "`idct(dct(x,n=2,t=3),n=4, t=3)/(4*4)` failed.") end subroutine test_modernized_idct diff --git a/test/test_fftpack_qct.f90 b/test/test_fftpack_qct.f90 deleted file mode 100644 index ed56b58..0000000 --- a/test/test_fftpack_qct.f90 +++ /dev/null @@ -1,77 +0,0 @@ -module test_fftpack_qct - - use fftpack, only: rk, dcosqi, dcosqf, dcosqb, qct, iqct - use testdrive, only: new_unittest, unittest_type, error_type, check - implicit none - private - - public :: collect_qct - -contains - - !> Collect all exported unit tests - subroutine collect_qct(testsuite) - !> Collection of tests - type(unittest_type), allocatable, intent(out) :: testsuite(:) - - testsuite = [ & - new_unittest("classic-qct-API", test_classic_qct), & - new_unittest("modernized-qct-API", test_modernized_qct), & - new_unittest("modernized-iqct-API", test_modernized_iqct) & - ] - - end subroutine collect_qct - - subroutine test_classic_qct(error) - type(error_type), allocatable, intent(out) :: error - real(kind=rk) :: w(3*4 + 15) - real(kind=rk) :: x(4) = [1, 2, 3, 4] - real(kind=rk) :: eps = 1.0e-10_rk - - 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 dcosqb(4, x, w) - call check(error, sum(abs(x/(4.0_rk*4.0_rk) - [real(kind=rk) :: 1, 2, 3, 4])) < eps, & - "`dcosqb` failed.") - - end subroutine test_classic_qct - - subroutine test_modernized_qct(error) - type(error_type), allocatable, intent(out) :: error - real(kind=rk) :: eps = 1.0e-10_rk - real(kind=rk) :: x(3) = [9, -9, 3] - - call check(error, sum(abs(qct(x, 2) - [-3.7279220613578570_rk, 21.727922061357859_rk])) < eps, & - "`qct(x, 2)` failed.") - if (allocated(error)) return - call check(error, sum(abs(qct(x, 3) - qct(x))) < eps, & - "`qct(x,3)` failed.") - if (allocated(error)) return - call check(error, sum(abs(qct(x, 4) - [-3.3871908980838743_rk, -2.1309424696909023_rk, & - 11.645661095452331_rk, 29.872472272322447_rk])) < eps, & - "`qct(x, 4)` failed.") - - end subroutine test_modernized_qct - - subroutine test_modernized_iqct(error) - type(error_type), allocatable, intent(out) :: error - real(kind=rk) :: eps = 1.0e-10_rk - real(kind=rk) :: x(4) = [1, 2, 3, 4] - - call check(error, sum(abs(iqct(qct(x))/(4.0_rk*4.0_rk) - [real(kind=rk) :: 1, 2, 3, 4])) < eps, & - "`iqct(qct(x)/(4.0_rk*4.0_rk)` failed.") - if (allocated(error)) return - call check(error, sum(abs(iqct(qct(x), 2)/(4.0_rk*2.0_rk) - [1.4483415291679655_rk, 7.4608849947753271_rk])) < eps, & - "`iqct(qct(x), 2)/(4.0_rk*2.0_rk)` failed.") - if (allocated(error)) return - call check(error, sum(abs(iqct(qct(x, 2), 4)/(4.0_rk*4.0_rk) - [0.5_rk, 0.70932417358418376_rk, 1.0_rk, & - 0.78858050747473762_rk])) < eps, & - "`iqct(qct(x, 2), 4)/(4.0_rk*4.0_rk)` failed.") - - end subroutine test_modernized_iqct - -end module test_fftpack_qct From cb99ff51d1eec04bd564f2406b13f304a527723d Mon Sep 17 00:00:00 2001 From: Chris Val Date: Thu, 2 Mar 2023 12:27:18 -0500 Subject: [PATCH 3/5] update & improve DCT docs --- doc/specs/fftpack.md | 747 ++++++++++++++++++++----------------------- 1 file changed, 345 insertions(+), 402 deletions(-) diff --git a/doc/specs/fftpack.md b/doc/specs/fftpack.md index 4a0e7fe..8274891 100644 --- a/doc/specs/fftpack.md +++ b/doc/specs/fftpack.md @@ -2,32 +2,30 @@ title: FFTPACK --- -# The `fftpack` Module - [TOC] -## Fourier transform of double complex periodic sequences -### `zffti` +# Discrete Fourier transform (DFT) of complex data +## `zffti` -#### Description +### Description Initializes the array `wsave` which is used in both `zfftf` and `zfftb`. The prime factorization of `n` together with a tabulation of the trigonometric functions are computed and stored in `wsave`. -#### Status +### Status Experimental. -#### Class +### Class Pure subroutine. -#### Syntax +### Syntax `call [[fftpack(module):zffti(interface)]](n, wsave)` -#### Argument +### Argument `n`: Shall be an `integer` scalar. This argument is `intent(in)`. @@ -40,11 +38,11 @@ The same work array can be used for both `zfftf` and `zfftb` as long as `n` remains unchanged. Different `wsave` arrays are required for different values of `n`. -##### Warning +#### Warning The contents of `wsave` must not be changed between calls of `zfftf` or `zfftb`. -#### Example +### Example ```fortran program demo_zffti @@ -55,9 +53,9 @@ program demo_zffti end program demo_zffti ``` -### `zfftf` +## `zfftf` -#### Description +### Description Computes the forward complex discrete fourier transform (the fourier analysis). Equivalently, `zfftf` computes the fourier coefficients of a complex periodic sequence. @@ -67,19 +65,19 @@ The transform is not normalized. To obtain a normalized transform the output mus The array `wsave` which is used by subroutine `zfftf` must be initialized by calling subroutine `zffti(n,wsave)`. -#### Status +### Status Experimental. -#### Class +### Class Pure subroutine. -#### Syntax +### Syntax `call [[fftpack(module):zfftf(interface)]](n, c, wsave)` -#### Argument +### Argument `n`: Shall be an `integer` scalar. This argument is `intent(in)`. @@ -106,11 +104,11 @@ This initialization does not have to be repeated so long as `n` remains unchange The same `wsave` array can be used by `zfftf` and `zfftb`. Contains initialization calculations which must not be destroyed between calls of subroutine `zfftf` or `zfftb`. -##### Warning +#### Warning The contents of `wsave` must not be changed between calls of `zfftf` or `zfftb`. -#### Example +### Example ```fortran program demo_zfftf @@ -123,9 +121,9 @@ program demo_zfftf end program demo_zfftf ``` -### `zfftb` +## `zfftb` -#### Description +### Description Unnormalized inverse of `zfftf`. @@ -137,19 +135,19 @@ The transform is not normalized. to obtain a normalized transform the output mus The array `wsave` which is used by subroutine `zfftf` must be initialized by calling subroutine `zffti(n,wsave)`. -#### Status +### Status Experimental. -#### Class +### Class Pure subroutine. -#### Syntax +### Syntax `call [[fftpack(module):zfftb(interface)]](n, c, wsave)` -#### Argument +### Argument `n`: Shall be an `integer` scalar. This argument is `intent(in)`. @@ -173,11 +171,11 @@ This argument is `intent(in)`. A `real` work array which must be dimensioned at least `4n+15` in the program that calls `zfftf`. The `wsave` array must be initialized by calling subroutine `zffti(n,wsave)` and a different `wsave` array must be used for each different value of `n`. This initialization does not have to be repeated so long as `n` remains unchanged thus subsequent transforms can be obtained faster than the first. The same `wsave` array can be used by `zfftf` and `zfftb`. Contains initialization calculations which must not be destroyed between calls of subroutine `zfftf` or `zfftb`. -##### Warning +#### Warning The contents of `wsave` must not be changed between calls of `zfftf` or `zfftb`. -#### Example +### Example ```fortran program demo_zfftb @@ -191,25 +189,25 @@ program demo_zfftb end program demo_zfftb ``` -### `fft` +## `fft` -#### Description +### Description Computes the forward complex discrete fourier transform (the fourier analysis). -#### Status +### Status Experimental. -#### Class +### Class Pure function. -#### Syntax +### Syntax `result = [[fftpack(module):fft(interface)]](x [, n])` -#### Argument +### Argument `x`: Shall be a `complex` and rank-1 array. This argument is `intent(in)`. @@ -218,15 +216,15 @@ This argument is `intent(in)`. This argument is `intent(in)` and `optional`. Defines the length of the Fourier transform. If `n` is not specified (the default) then `n = size(x)`. If `n <= size(x)`, `x` is truncated, if `n > size(x)`, `x` is zero-padded. -#### Return value +### Return value Returns a `complex` and rank-1 array, the Discrete Fourier Transform (DFT) of `x`. -#### Notes +### Notes Within numerical accuracy, `x == ifft(fft(x))/size(x)`. -#### Example +### Example ```fortran program demo_fft @@ -239,25 +237,25 @@ program demo_fft end program demo_fft ``` -### `ifft` +## `ifft` -#### Description +### Description Unnormalized inverse of `fft`. -#### Status +### Status Experimental. -#### Class +### Class Pure function. -#### Syntax +### Syntax `result = [[fftpack(module):ifft(interface)]](x [, n])` -#### Argument +### Argument `x`: Shall be a `complex` and rank-1 array. This argument is `intent(in)`. @@ -266,11 +264,11 @@ This argument is `intent(in)`. This argument is `intent(in)` and `optional`. Defines the length of the Fourier transform. If `n` is not specified (the default) then `n = size(x)`. If `n <= size(x)`, `x` is truncated, if `n > size(x)`, `x` is zero-padded. -#### Return value +### Return value Returns a `complex` and rank-1 array, the unnormalized inverse Discrete Fourier Transform (DFT) of `x`. -#### Example +### Example ```fortran program demo_ifft @@ -281,28 +279,28 @@ program demo_ifft end program demo_ifft ``` -## Fourier transform of double real periodic sequences -### `dffti` +# Discrete Fourier transform (DFT) of real data +## `dffti` -#### Description +### Description Initializes the array `wsave` which is used in both `dfftf` and `dfftb`. The prime factorization of `n` together with a tabulation of the trigonometric functions are computed and stored in `wsave`. -#### Status +### Status Experimental. -#### Class +### Class Pure subroutine. -#### Syntax +### Syntax `call [[fftpack(module):dffti(interface)]](n, wsave)` -#### Argument +### Argument `n`: Shall be an `integer` scalar. This argument is `intent(in)`. @@ -312,13 +310,13 @@ The length of the sequence to be transformed. This argument is `intent(out)`. A work array which must be dimensioned at least `2*n+15`. The same work array can be used for both `dfftf` and `dfftb` as long as `n` remains unchanged. -Different `wsave` arrays are required for different values of `n`. +Different `wsave` arrays are required for different values of `n`. -##### Warning +#### Warning The contents of `wsave` must not be changed between calls of `dfftf` or `dfftb`. -#### Example +### Example ```fortran program demo_dffti @@ -329,9 +327,9 @@ program demo_dffti end program demo_dffti ``` -### `dfftf` +## `dfftf` -#### Description +### Description Computes the fourier coefficients of a real perodic sequence (fourier analysis). The transform is defined below at output parameter `r`. @@ -340,19 +338,19 @@ The transform is not normalized. To obtain a normalized transform the output mus The array `wsave` which is used by subroutine `dfftf` must be initialized by calling subroutine `dffti(n,wsave)`. -#### Status +### Status Experimental. -#### Class +### Class Pure subroutine. -#### Syntax +### Syntax `call [[fftpack(module):dfftf(interface)]](n, r, wsave)` -#### Argument +### Argument `n`: Shall be an `integer` scalar. This argument is `intent(in)`. @@ -392,11 +390,11 @@ This initialization does not have to be repeated so long as `n` remains unchange The same `wsave` array can be used by `dfftf` and `dfftb`. Contains initialization calculations which must not be destroyed between calls of subroutine `dfftf` or `dfftb`. -##### Warning +#### Warning The contents of `wsave` must not be changed between calls of `dfftf` or `dfftb`. -#### Example +### Example ```fortran program demo_dfftf @@ -408,9 +406,9 @@ program demo_dfftf end program demo_dfftf ``` -### `dfftb` +## `dfftb` -#### Description +### Description Unnormalized inverse of `dfftf`. @@ -422,19 +420,19 @@ The transform is not normalized. To obtain a normalized transform the output mus The array `wsave` which is used by subroutine `dfftf` must be initialized by calling subroutine `dffti(n,wsave)`. -#### Status +### Status Experimental. -#### Class +### Class Pure subroutine. -#### Syntax +### Syntax `call [[fftpack(module):dfftb(interface)]](n, r, wsave)` -#### Argument +### Argument `n`: Shall be an `integer` scalar. This argument is `intent(in)`. @@ -468,11 +466,11 @@ This argument is `intent(in)`. A `real` work array which must be dimensioned at least `2n+15` in the program that calls `dfftf`. The `wsave` array must be initialized by calling subroutine `dffti(n,wsave)` and a different `wsave` array must be used for each different value of `n`. This initialization does not have to be repeated so long as `n` remains unchanged thus subsequent transforms can be obtained faster than the first. The same `wsave` array can be used by `dfftf` and `dfftb`. Contains initialization calculations which must not be destroyed between calls of subroutine `dfftf` or `dfftb`. -##### Warning +#### Warning The contents of `wsave` must not be changed between calls of `dfftf` or `dfftb`. -#### Example +### Example ```fortran program demo_dfftb @@ -485,25 +483,25 @@ program demo_dfftb end program demo_dfftb ``` -### `rfft` +## `rfft` -#### Description +### Description Discrete Fourier transform of a real sequence. -#### Status +### Status Experimental. -#### Class +### Class Pure function. -#### Syntax +### Syntax `result = [[fftpack(module):rfft(interface)]](x [, n])` -#### Argument +### Argument `x`: Shall be a `real` and rank-1 array. This argument is `intent(in)`. @@ -513,15 +511,15 @@ The data to transform. This argument is `intent(in)` and `optional`. Defines the length of the Fourier transform. If `n` is not specified (the default) then `n = size(x)`. If `n <= size(x)`, `x` is truncated, if `n > size(x)`, `x` is zero-padded. -#### Return value +### Return value Returns a `real` and rank-1 array, the Discrete Fourier Transform (DFT) of `x`. -#### Notes +### Notes Within numerical accuracy, `y == rfft(irfft(y))/size(y)`. -#### Example +### Example ```fortran program demo_rfft @@ -533,25 +531,25 @@ program demo_rfft end program demo_rfft ``` -### `irfft` +## `irfft` -#### Description +### Description Unnormalized inverse of `rfft`. -#### Status +### Status Experimental. -#### Class +### Class Pure function. -#### Syntax +### Syntax `result = [[fftpack(module):irfft(interface)]](x [, n])` -#### Argument +### Argument `x`: Shall be a `real` array. This argument is `intent(in)`. @@ -561,11 +559,11 @@ Transformed data to invert. This argument is `intent(in)` and `optional`. Defines the length of the Fourier transform. If `n` is not specified (the default) then `n = size(x)`. If `n <= size(x)`, `x` is truncated, if `n > size(x)`, `x` is zero-padded. -#### Return value +### Return value Returns a `real` and rank-1 array, the unnormalized inverse discrete Fourier transform. -#### Example +### Example ```fortran program demo_irfft @@ -576,28 +574,28 @@ program demo_irfft end program demo_irfft ``` -## Simplified fourier transform of double real periodic sequences +# Simplified discrete Fourier transform (DFT) of real data -### `dzffti` +## `dzffti` -#### Description +### Description Initializes the array `wsave` which is used in both `dzfftf` and `dzfftb`. The prime factorization of `n` together with a tabulation of the trigonometric functions are computed and stored in `wsave`. -#### Status +### Status Experimental -#### Class +### Class Prue function. -#### Syntax +### Syntax `call [[fftpack(module):dzffti(interface)]](n, wsave)` -#### Arguments +### Arguments `n`: Shall be an `integer` scalar. This argument is `intent(in)`. @@ -609,11 +607,11 @@ A work array which must be dimensioned at least `3*n+15`. The same work array can be used for both `dzfftf` and `dzfftb` as long as n remains unchanged. Different `wsave` arrays are required for different values of `n`. -##### Warning +#### Warning The contents of `wsave` must not be changed between calls of `dzfftf` or `dzfftb`. -#### Example +### Example ```fortran program demo_dzffti @@ -624,27 +622,27 @@ program demo_dzffti end program demo_dzffti ``` -### `dzfftf` +## `dzfftf` -#### Description +### Description Computes the fourier coefficients of a `real` perodic sequence (fourier analysis). The transform is defined below at output parameters `azero`, `a` and `b`. `dzfftf` is a simplified but **slower version** of `dfftf`. -#### Status +### Status Experimental -#### Class +### Class Pure subroutine. -#### Syntax +### Syntax `call [[fftpack(module):dzfftf(interface)]](n, r, azero, a, b, wsave)` -#### Arguments +### Arguments `n`: Shall be an `integer` scalar. This argument is `intent(in)`. @@ -685,7 +683,7 @@ In the program that calls `dzfftf`. The `wsave` array must be initialized by cal This initialization does not have to be repeated so long as `n` remains unchanged thus subsequent transforms can be obtained faster than the first. The same `wsave` array can be used by `dzfftf` and `dzfftb`. -#### Example +### Example ```fortran program demo_dzfftf @@ -698,27 +696,27 @@ program demo_dzfftf end program demo_dzfftf ``` -### `dzfftb` +## `dzfftb` -#### Description +### Description Computes a `real` perodic sequence from its fourier coefficients (fourier synthesis). The transform is defined below at output parameter `r`. `dzfftb` is a simplified but **slower version** of `dfftb`. -#### Status +### Status Experimental -#### Class +### Class Pure subroutine. -#### Syntax +### Syntax `call [[fftpack(module):dzfftb(interface)]](n, r, azero, a, b, wsave)` -#### Arguments +### Arguments `n`: Shall be an `integer` scalar. This argument is `intent(in)`. @@ -792,7 +790,7 @@ In the program that calls `dzfftf`. The `wsave` array must be initialized by cal This initialization does not have to be repeated so long as `n` remains unchanged thus subsequent transforms can be obtained faster than the first. The same `wsave` array can be used by `dzfftf` and `dzfftb`. -#### Example +### Example ```fortran program demo_dzfftb @@ -807,16 +805,16 @@ program demo_dzfftb end program demo_dzfftb ``` -## Cosine transform with odd wave numbers +# Discrete cosine transforms (DCT) + +## DCT type-1 (DCT-1) -### `dcosqi` +### Initialize DCT-1: `dcosti` #### Description -Initializes the array `wsave` which is used in both `dcosqf` and `dcosqb`. -The prime factorization of `n` together with -a tabulation of the trigonometric functions are computed and -stored in `wsave`. +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`. #### Status @@ -828,44 +826,49 @@ Pure subroutine. #### Syntax -`call [[fftpack(module):dcosqi(interface)]](n, wsave)` +`call [[fftpack(module):dcosti(interface)]](n , wsave)` #### Arguments -`n`: Shall be an `integer` scalar. +`n`: Shall be a `integer` scalar. This argument is `intent(in)`. -The length of the array to be transformed. -The method is most efficient when `n` is a product of small primes. +The length of the sequence to be transformed. +The method is most efficient when n-1 is a product of small primes. `wsave`: Shall be a `real` and rank-1 array. This argument is `intent(out)`. A work array which must be dimensioned at least `3*n+15`. -The same work array can be used for both `dcosqf` and `dcosqb` -as long as `n` remains unchanged. Different `wsave` arrays are required for different values of `n`. -The contents of `wsave` must not be changed between calls of `dcosqf` or `dcosqb`. +The contents of `wsave` must not be changed between calls of `dcost`. #### Example ```fortran -program demo_dcosqi - use fftpack, only: dcosqi +program demo_dcosti + use fftpack, only: dcosti real(kind=8) :: w(3*4 + 15) - call dcosqi(4, w) !! Initializes the array `w` which is used in both `dcosqf` and `dcosqb`. -end program demo_dcosqi + call dcosti(4, w) !! Initializes the array `w` which is used in subroutine `dcost`. +end program demo_dcosti ``` -### `dcosqf` +### Compute DCT-1: `dcost` -#### Decsription +#### Description -Computes the fast fourier transform of quarter wave data. -That is, `dcosqf` computes the coefficients in a cosine series representation with only odd wave numbers. +Computes the DCT-1 of the input real data. The transform is defined below at output parameter `x`. -`dcosqf` is the unnormalized inverse of `dcosqb` since a call of `dcosqf` followed by a call of `dcosqb` will multiply the input sequence `x` by `4*n`. +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\]). -The array `wsave` which is used by subroutine `dcosqf` must be initialized by calling subroutine `dcosqi(n,wsave)`. +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` +by `2*(n-1)`. + +The array `wsave` which is used by subroutine `dcost` must be initialized by calling subroutine `dcosti(n,wsave)`. #### Status @@ -877,65 +880,66 @@ Pure subroutine. #### Syntax -`call [[fftpack(module):dcosqf(interface)]](n, x, wsave)` +`call [[fftpack(module):dcost(interface)]](n, x, wsave)` #### Arguments -`n`: Shall be an `integer` scalar. +`n`: Shall be a `integer` scalar. This argument is `intent(in)`. -The length of the array `x` to be transformed. -The method is most efficient when `n` is a product of small primes. +The length of the sequence `x`. +`n` must be greater than `1`. +The method is most efficient when `n-1` is a product of small primes. `x`: Shall be a `real` and rank-1 array. -This argument is `intent(inout)`. -An array which contains the sequence to be transformed. +This argument is `intent(inout)`. +An array which contains the sequence to be transformed, and is overwritten +by the result. ``` for i=1,...,n - x(i) = x(1) plus the sum from k=2 to k=n of + x(i) = x(1)+(-1)**(i-1)*x(n) - 2*x(k)*cos((2*i-1)*(k-1)*pi/(2*n)) + + the sum from k=2 to k=n-1 - a call of dcosqf followed by a call of - cosqb will multiply the sequence x by 4*n. - therefore dcosqb is the unnormalized inverse - of dcosqf. + 2*x(k)*cos((k-1)*(i-1)*pi/(n-1)) + + a call of dcost followed by another call of + dcost will multiply the sequence x by 2*(n-1) + hence dcost is the unnormalized inverse + of itself. ``` `wsave`: Shall be a `real` and rank-1 array. This argument is `intent(in)`. -A work array which must be dimensioned at least `3*n+15` -in the program that calls `dcosqf`. -The `wsave` array must be initialized by calling subroutine `dcosqi(n,wsave)` and a different `wsave` array must be used for each different value of `n`. -This initialization does not have to be repeated so long as `n` remains unchanged thus subsequent transforms can be obtained faster than the first. - -##### Warning - -`wsave` contains initialization calculations which must not be destroyed between calls of `dcosqf` or `dcosqb`. +A work array which must be of length at least `3*n+15` in the program that calls `dcost`. +The `wsave` array must be initialized by calling subroutine `dcosti(n,wsave)` and a different `wsave` array must be used for each different value of `n`. +This initialization does not have to be repeated so long as `n` remains unchanged thus subsequent +transforms can be obtained faster than the first. +Contains initialization calculations which must not be destroyed between calls of `dcost`. #### Example ```fortran -program demo_dcosqf - use fftpack, only: dcosqi, dcosqf - real(kind=8) :: w(3*4 + 15) +program demo_dcost + use fftpack, only: dcosti, dcost real(kind=8) :: x(4) = [1, 2, 3, 4] - call dcosqi(4, w) - call dcosqf(4, x, w) !! `x`: [12.0, -9.10, 2.62, -1.51] -end program demo_dcosqf + real(kind=8) :: w(3*4 + 15) + call dcosti(4, w) + call dcost(4, x, w) !! `x`: [15.0, -4.0, 0.0, -1.0] + call dcost(4, x, w) !! `x`: [6.0, 12.0, 18.0, 24.0] +end program demo_dcost ``` -### `dcosqb` +## DCT of types 2, 3 (DCT-2, 3), a.k.a "Quarter" cosine transforms -#### Decsription +### Initialize DCT-2, 3: `dcosqi` -Computes the fast fourier transform of quarter wave data. -That is, `dcosqb` computes a sequence from its representation in terms of a cosine series with odd wave numbers. -The transform is defined below at output parameter `x`. - -`dcosqb` is the unnormalized inverse of `dcosqf` since a call of `dcosqb` followed by a call of `dcosqf` will multiply the input sequence `x` by `4*n`. +#### Description -The array `wsave` which is used by subroutine `dcosqb` must be initialized by calling subroutine `dcosqi(n,wsave)`. +Initializes the array `wsave` which is used in both `dcosqf` and `dcosqb`. +The prime factorization of `n` together with +a tabulation of the trigonometric functions are computed and +stored in `wsave`. #### Status @@ -947,201 +951,120 @@ Pure subroutine. #### Syntax -`call [[fftpack(module):dcosqf(interface)]](n, x, wsave)` +`call [[fftpack(module):dcosqi(interface)]](n, wsave)` #### Arguments `n`: Shall be an `integer` scalar. This argument is `intent(in)`. -The length of the array `x` to be transformed. +The length of the array to be transformed. The method is most efficient when `n` is a product of small primes. -`x`: Shall be a `real` and rank-1 array. -This argument is `intent(inout)`. -An array which contains the sequence to be transformed. -``` -for i=1,...,n - - x(i)= the sum from k=1 to k=n of - - 4*x(k)*cos((2*k-1)*(i-1)*pi/(2*n)) - - a call of dcosqb followed by a call of - dcosqf will multiply the sequence x by 4*n. - therefore dcosqf is the unnormalized inverse - of dcosqb. -``` - `wsave`: Shall be a `real` and rank-1 array. -This argument is `intent(in)`. -A work array which must be dimensioned at least `3*n+15` -in the program that calls `dcosqb`. -The `wsave` array must be initialized by calling subroutine `dcosqi(n,wsave)` and a different `wsave` array must be used for each different value of `n`. -This initialization does not have to be repeated so long as `n` remains unchanged thus subsequent transforms can be obtained faster than the first. - -##### Warning - -`wsave` contains initialization calculations which must not be destroyed between calls of `dcosqf` or `dcosqb`. +This argument is `intent(out)`. +A work array which must be dimensioned at least `3*n+15`. +The same work array can be used for both `dcosqf` and `dcosqb` +as long as `n` remains unchanged. +Different `wsave` arrays are required for different values of `n`. +The contents of `wsave` must not be changed between calls of `dcosqf` or `dcosqb`. #### Example ```fortran -program demo_dcosqb - use fftpack, only: dcosqi, dcosqf, dcosqb +program demo_dcosqi + use fftpack, only: dcosqi real(kind=8) :: w(3*4 + 15) - real(kind=8) :: x(4) = [4, 3, 5, 10] - call dcosqi(4, w) - call dcosqf(4, x, w) - call dcosqb(4, x, w) !! `x`: [1.0, 2.0, 3.0, 4.0] * 4 * n, n = 4, which is unnormalized. -end program demo_dcosqb + call dcosqi(4, w) !! Initializes the array `w` which is used in both `dcosqf` and `dcosqb`. +end program demo_dcosqi ``` -### `qct` +### Compute DCT-3: `dcosqf` #### Description -Forward transform of quarter wave data. - -#### Status - -Experimental. - -#### Class - -Pure function. - -#### Syntax - -`result = [[fftpack(module):qct(interface)]](x [, n])` - -#### Argument - -`x`: Shall be a `real` and rank-1 array. -This argument is `intent(in)`. -The data to transform. - -`n`: Shall be an `integer` scalar. -This argument is `intent(in)` and `optional`. -Defines the length of the Fourier transform. If `n` is not specified (the default) then `n = size(x)`. If `n <= size(x)`, `x` is truncated, if `n > size(x)`, `x` is zero-padded. - -#### Return value - -Returns a `real` and rank-1 array, the Quarter-Cosine Transform (QCT) of `x`. - -#### Notes - -Within numerical accuracy, `x == iqct(qct(x))/(4*size(x))`. - -#### Example - -```fortran -program demo_qct - use fftpack, only: qct - real(kind=8) :: x(4) = [1, 2, 3, 4] - print *, qct(x,3) !! [7.4, -5.0, 0.53]. - print *, qct(x) !! [12.0, -9.10, 2.62, -1.51]. - print *, qct(x,5) !! [14.4, -6.11, -5.0, 4.4, -2.65]. -end program demo_qct -``` - -### `iqct` +Computes the DCT-3 of the input real data. +The transform is defined below at output parameter `x`. -#### Description +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`. -Unnormalized inverse of `qct`. +The array `wsave` which is used by subroutine `dcosqf` must be initialized by calling subroutine `dcosqi(n,wsave)`. #### Status -Experimental. +Experimental #### Class -Pure function. +Pure subroutine. #### Syntax -`result = [[fftpack(module):iqct(interface)]](x [, n])` - -#### Argument +`call [[fftpack(module):dcosqf(interface)]](n, x, wsave)` -`x`: Shall be a `real` array. -This argument is `intent(in)`. -Transformed data to invert. +#### Arguments `n`: Shall be an `integer` scalar. -This argument is `intent(in)` and `optional`. -Defines the length of the Fourier transform. If `n` is not specified (the default) then `n = size(x)`. If `n <= size(x)`, `x` is truncated, if `n > size(x)`, `x` is zero-padded. - -#### Return value - -Returns a `real` and rank-1 array, the unnormalized inverse Quarter-Cosine Transform. - -#### Example +This argument is `intent(in)`. +The length of the array `x` to be transformed. +The method is most efficient when `n` is a product of small primes. -```fortran -program demo_iqct - use fftpack, only: qct, iqct - real(kind=8) :: x(4) = [1, 2, 3, 4] - print *, iqct(qct(x))/(4.0*4.0) !! [1.0, 2.0, 3.0, 4.0] - print *, iqct(qct(x), 3)/(4.0*3.0) !! [1.84, 2.71, 5.47] -end program demo_iqct +`x`: Shall be a `real` and rank-1 array. +This argument is `intent(inout)`. +An array which contains the sequence to be transformed, and is overwritten by +the result. ``` +for i=1,...,n -## Cosine transform of a real even sequence - -### `dcosti` - -#### 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`. - -#### Status - -Experimental - -#### Class - -Pure subroutine. - -#### Syntax + x(i) = x(1) plus the sum from k=2 to k=n of -`call [[fftpack(module):dcosti(interface)]](n , wsave)` + 2*x(k)*cos((2*i-1)*(k-1)*pi/(2*n)) -#### Arguments + a call of dcosqf followed by a call of + cosqb will multiply the sequence x by 4*n. + therefore dcosqb is the unnormalized inverse + of dcosqf. +``` -`n`: Shall be a `integer` scalar. +`wsave`: Shall be a `real` and rank-1 array. This argument is `intent(in)`. -The length of the sequence to be transformed. -The method is most efficient when n-1 is a product of small primes. +A work array which must be dimensioned at least `3*n+15` +in the program that calls `dcosqf`. +The `wsave` array must be initialized by calling subroutine `dcosqi(n,wsave)` and a different `wsave` array must be used for each different value of `n`. +This initialization does not have to be repeated so long as `n` remains unchanged thus subsequent transforms can be obtained faster than the first. -`wsave`: Shall be a `real` and rank-1 array. -This argument is `intent(out)`. -A work array which must be dimensioned at least `3*n+15`. -Different `wsave` arrays are required for different values of `n`. -The contents of `wsave` must not be changed between calls of `dcost`. +**Warning**: `wsave` contains initialization calculations which must not be +destroyed between calls of `dcosqf` or `dcosqb` of the same `n`. #### Example ```fortran -program demo_dcosti - use fftpack, only: dcosti +program demo_dcosqf + use fftpack, only: dcosqi, dcosqf real(kind=8) :: w(3*4 + 15) - call dcosti(4, w) !! Initializes the array `w` which is used in subroutine `dcost`. -end program demo_dcosti + real(kind=8) :: x(4) = [1, 2, 3, 4] + call dcosqi(4, w) + call dcosqf(4, x, w) !! `x`: [12.0, -9.10, 2.62, -1.51] +end program demo_dcosqf ``` -### `dcost` +### Compute DCT-2: `dcosqb` #### Description -Computes the discrete fourier cosine transform of an even sequence `x(i)`. +Computes the DCT-2 of the input real data. The transform is defined below at output parameter `x`. -`dcost` is the unnormalized inverse of itself since a call of `dcost` followed by another call of `dcost` will multiply the input sequence `x` by `2*(n-1)`. -The transform is defined below at output parameter `x`. +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. -The array `wsave` which is used by subroutine `dcost` must be initialized by calling subroutine `dcosti(n,wsave)`. +Also, `dcosqb` (DCT-2) is the unnormalized inverse of `dcosqf` (DCT-3), since a +call of `dcosqb` followed by a call of `dcosqf` will multiply the input sequence +`x` by `4*n`. + +The array `wsave` which is used by subroutine `dcosqb` must be initialized by calling subroutine `dcosqi(n,wsave)`. #### Status @@ -1153,124 +1076,134 @@ Pure subroutine. #### Syntax -`call [[fftpack(module):dcost(interface)]](n, x, wsave)` +`call [[fftpack(module):dcosqf(interface)]](n, x, wsave)` #### Arguments -`n`: Shall be a `integer` scalar. +`n`: Shall be an `integer` scalar. This argument is `intent(in)`. -The length of the sequence `x`. -`n` must be greater than `1`. -The method is most efficient when `n-1` is a product of small primes. +The length of the array `x` to be transformed. +The method is most efficient when `n` is a product of small primes. `x`: Shall be a `real` and rank-1 array. -This argument is `intent(inout)`. -An array which contains the sequence to be transformed. +This argument is `intent(inout)`. +An array which contains the sequence to be transformed, and is overwritten by +the result. ``` for i=1,...,n - x(i) = x(1)+(-1)**(i-1)*x(n) - - + the sum from k=2 to k=n-1 + x(i)= the sum from k=1 to k=n of - 2*x(k)*cos((k-1)*(i-1)*pi/(n-1)) + 4*x(k)*cos((2*k-1)*(i-1)*pi/(2*n)) - a call of dcost followed by another call of - dcost will multiply the sequence x by 2*(n-1) - hence dcost is the unnormalized inverse - of itself. + a call of dcosqb followed by a call of + dcosqf will multiply the sequence x by 4*n. + therefore dcosqf is the unnormalized inverse + of dcosqb. ``` `wsave`: Shall be a `real` and rank-1 array. This argument is `intent(in)`. -A work array which must be dimensioned at least `3*n+15` in the program that calls `dcost`. -The `wsave` array must be initialized by calling subroutine `dcosti(n,wsave)` and a different `wsave` array must be used for each different value of `n`. -This initialization does not have to be repeated so long as `n` remains unchanged thus subsequent -transforms can be obtained faster than the first. -Contains initialization calculations which must not be destroyed between calls of `dcost`. +A work array which must be dimensioned at least `3*n+15` +in the program that calls `dcosqb`. +The `wsave` array must be initialized by calling subroutine `dcosqi(n,wsave)` and a different `wsave` array must be used for each different value of `n`. +This initialization does not have to be repeated so long as `n` remains unchanged thus subsequent transforms can be obtained faster than the first. + +**Warning**: `wsave` contains initialization calculations which must not be +destroyed between calls of `dcosqf` or `dcosqb` of the same `n`. #### Example ```fortran -program demo_dcost - use fftpack, only: dcosti, dcost - real(kind=8) :: x(4) = [1, 2, 3, 4] +program demo_dcosqb + use fftpack, only: dcosqi, dcosqf, dcosqb real(kind=8) :: w(3*4 + 15) - call dcosti(4, w) - call dcost(4, x, w) !! Computes the discrete fourier cosine (forward) transform of an even sequence, `x`(unnormalized): [15.0, -4.0, 0.0, -1.0] - call dcost(4, x, w) !! Computes the discrete fourier cosine (backward) transform of an even sequence, `x`(unnormalized): [6.0, 12.0, 18.0, 24.0] -end program demo_dcost + real(kind=8) :: x(4) = [4, 3, 5, 10] + call dcosqi(4, w) + call dcosqf(4, x, w) + call dcosqb(4, x, w) !! [64.0, 48.0, 80.0, 160.0] +end program demo_dcosqb ``` -### `dct` +## Simplified DCT of types 1, 2, 3: `dct` -#### Description +### Description -Discrete fourier cosine (forward) transform of an even sequence. +Discrete cosine transforms (DCT) of types 1, 2, 3. +This is a more flexible interface for the DCT of types 1, 2 and 3, albeit +slightly slower than the in-place DCT procedures. -#### Status +### Status Experimental. -#### Class +### Class Pure function. -#### Syntax +### Syntax -`result = [[fftpack(module):dct(interface)]](x [, n])` +`result = [[fftpack(module):dct(interface)]](x [, n, t])` -#### Argument +### Argument `x`: Shall be a `real` and rank-1 array. -This argument is `intent(in)`. +This argument is `intent(in)`. The data to transform. `n`: Shall be an `integer` scalar. -This argument is `intent(in)` and `optional`. -Defines the length of the Fourier transform. If `n` is not specified (the default) then `n = size(x)`. If `n <= size(x)`, `x` is truncated, if `n > size(x)`, `x` is zero-padded. +This argument is `intent(in)` and `optional`. +Defines the length of the DCT. If `n` is not specified (the default) then `n = size(x)`. If `n <= size(x)`, `x` is truncated, if `n > size(x)`, `x` is zero-padded. + +`t`: Shall be an `integer` scalar, equal to `1`, `2` or `3`. +This argument is `intent(in)` and `optional`. +Defines the type of DCT to be performed. The default type is `2`. #### Return value -Returns a `real` and rank-1 array, the Discrete-Cosine Transform (DCT) of `x`. +Returns a `real` and rank-1 array, the DCT type-`t` of the input data `x`. #### Notes -Within numerical accuracy, `y == dct(idct(y))/2*(size(y) - 1)`. +Within numerical accuracy, +- `x == idct(dct(x, t=1), t=1) / (2*(size(x) - 1))` +- `x == idct(dct(x, t=2), t=2) / (4*size(x))` +- `x == idct(dct(x, t=3), t=3) / (4*size(x))` -#### Example +### Example ```fortran program demo_dct use fftpack, only: dct real(kind=8) :: x(4) = [1, 2, 3, 4] - print *, dct(x,3) !! [8.0, -2.0, 0.0]. - print *, dct(x) !! [15.0, -4.0, 0.0, -1.0]. - print *, dct(x,5) !! [19.0, -1.8, -5.0, 3.8, -5.0]. - print *, dct(dct(x))/(2*(4 - 1)) !! (normalized): [1.0, 2.0, 3.0, 4.0] + print *, dct(x,3,1) !! [8.0, -2.0, 0.0]. + print *, dct(x,t=1) !! [15.0, -4.0, 0.0, -1.0]. + print *, dct(x,5,2) !! [14.36, -6.11, -5.0, 4.40, -2.65]. + print *, dct(dct(x,t=1),t=1)/(2*(4 - 1)) !! (normalized): [1.0, 2.0, 3.0, 4.0] end program demo_dct ``` -### `idct` +## Simplified inverse DCT of types 1, 2, 3: `idct` -#### Description +### Description -Unnormalized inverse of `dct`. -In fact, `idct` and `dct` have the same effect, `idct` = `dct`. +Unnormalized inverse discrete cosine transform (IDCT) of types 1, 2 and 3. +This is a more flexible interface for the IDCT of types 1, 2 and 3, +albeit slightly slower than the in-place DCT procedures. -#### Status +### Status Experimental. -#### Class +### Class Pure function. -#### Syntax +### Syntax -`result = [[fftpack(module):idct(interface)]](x [, n])` +`result = [[fftpack(module):idct(interface)]](x [, n, t])` -#### Argument +### Argument `x`: Shall be a `real` array. This argument is `intent(in)`. @@ -1280,53 +1213,63 @@ Transformed data to invert. This argument is `intent(in)` and `optional`. Defines the length of the Fourier transform. If `n` is not specified (the default) then `n = size(x)`. If `n <= size(x)`, `x` is truncated, if `n > size(x)`, `x` is zero-padded. +`t`: Shall be an `integer` scalar, equal to `1` or `2`. +This argument is `intent(in)` and `optional`. +Defines the type of the IDCT to be performed. The default type is `2`. + #### Return value -Returns a `real` and rank-1 array, the inverse Discrete-Cosine Transform (iDCT) of `x`. +Returns a `real` and rank-1 array, the IDCT type-`t` of the input data `x`. -#### Example +#### Notes + +Within numerical accuracy, +- `x == idct(dct(x, t=1), t=1) / (2*(size(x) - 1))` +- `x == idct(dct(x, t=2), t=2) / (4*size(x))` +- `x == idct(dct(x, t=3), t=3) / (4*size(x))` + +### Example ```fortran program demo_idct use fftpack, only: dct, idct real(kind=8) :: x(4) = [1, 2, 3, 4] - print *, idct(dct(x))/(2*(4-1)) !! (normalized): [1.0, 2.0, 3.0, 4.0] - print *, idct(idct(x))/(2*(4-1)) !! (normalized): [1.0, 2.0, 3.0, 4.0] - print *, idct(dct(x), 3) !! (unnormalized): [7.0, 15.0, 23.0] + print *, idct(dct(x,t=1),t=1)/(2*(4-1)) !! (normalized): [1.0, 2.0, 3.0, 4.0] + print *, idct(dct(x,t=2),t=2)/(4*4) !! (normalized): [1.0, 2.0, 3.0, 4.0] + print *, idct(dct(x), n=3) !! (unnormalized): [22.06, 32.5, 65.65] end program demo_idct ``` +# Utility functions -## Utility functions +## `fftshift` -### `fftshift` - -#### Description +### Description Rearranges the Fourier transform by moving the zero-frequency component to the center of the array. -#### Status +### Status Experimental. -#### Class +### Class Pure function. -#### Syntax +### Syntax `result = [[fftpack(module):fftshift(interface)]](x)` -#### Argument +### Argument `x`: Shall be a `complex/real` and rank-1 array. This argument is `intent(in)`. -#### Return value +### Return value Returns the `complex/real` and rank-1 Fourier transform by moving the zero-frequency component to the center of the array. -#### Example +### Example ```fortran program demo_fftshift @@ -1340,34 +1283,34 @@ program demo_fftshift end program demo_fftshift ``` -### `ifftshift` +## `ifftshift` -#### Description +### Description Rearranges the Fourier transform with zero frequency shifting back to the original transform output. In other words, `ifftshift` is the result of undoing `fftshift`. -#### Status +### Status Experimental. -#### Class +### Class Pure function. -#### Syntax +### Syntax `result = [[fftpack(module):ifftshift(interface)]](x)` -#### Argument +### Argument `x`: Shall be a `complex/real` and rank-1 array. This argument is `intent(in)`. -#### Return value +### Return value Returns the `complex/real` and rank-1 Fourier transform with zero frequency shifting back to the original transform output. -#### Example +### Example ```fortran program demo_ifftshift @@ -1381,34 +1324,34 @@ program demo_ifftshift end program demo_ifftshift ``` -### `fftfreq` +## `fftfreq` -#### Description +### Description Returns the integer frequency (or wavenumber) values that correspond to the coefficients calculated by the complex discrete Fourier transform, in the standard order (zero frequency first). -#### Status +### Status Experimental. -#### Class +### Class Pure function. -#### Syntax +### Syntax `result = [[fftpack(module):fftfreq(interface)]](n)` -#### Argument +### Argument `n`: Shall be an `integer`, equal to the length of the corresponding complex discrete Fourier transform. This argument is `intent(in)`. -#### Return value +### Return value Returns the `integer` and rank-1 array of the transform's frequency values in the standard order (zero frequency first). -#### Example +### Example ```fortran program demo_fftfreq @@ -1418,34 +1361,34 @@ program demo_fftfreq end program demo_fftfreq ``` -### `rfftfreq` +## `rfftfreq` -#### Description +### Description Returns the integer frequency (or wavenumber) values that correspond to the coefficients calculated by the real discrete Fourier transform, in the standard order (zero frequency first). -#### Status +### Status Experimental. -#### Class +### Class Pure function. -#### Syntax +### Syntax `result = [[fftpack(module):rfftfreq(interface)]](n)` -#### Argument +### Argument `n`: Shall be an `integer`, equal to the length of the corresponding real discrete Fourier transform. This argument is `intent(in)`. -#### Return value +### Return value Returns the `integer` and rank-1 array of the transform's frequency values in the standard order (zero frequency first). -#### Example +### Example ```fortran program demo_rfftfreq From fa10e9cfdb3286a54c7330107ed72081fe0f3ce2 Mon Sep 17 00:00:00 2001 From: Chris Val Date: Fri, 3 Mar 2023 15:02:18 -0500 Subject: [PATCH 4/5] review 1 changes --- doc/specs/fftpack.md | 30 ++++++++++++------------- src/fftpack.f90 | 8 +++---- src/fftpack_dct.f90 | 24 ++++++++++---------- test/test_fftpack_dct.f90 | 46 +++++++++++++++++++++------------------ 4 files changed, 56 insertions(+), 52 deletions(-) diff --git a/doc/specs/fftpack.md b/doc/specs/fftpack.md index 8274891..1fd471f 100644 --- a/doc/specs/fftpack.md +++ b/doc/specs/fftpack.md @@ -1143,7 +1143,7 @@ Pure function. ### Syntax -`result = [[fftpack(module):dct(interface)]](x [, n, t])` +`result = [[fftpack(module):dct(interface)]](x [, n, type])` ### Argument @@ -1155,7 +1155,7 @@ The data to transform. This argument is `intent(in)` and `optional`. Defines the length of the DCT. If `n` is not specified (the default) then `n = size(x)`. If `n <= size(x)`, `x` is truncated, if `n > size(x)`, `x` is zero-padded. -`t`: Shall be an `integer` scalar, equal to `1`, `2` or `3`. +`type`: Shall be an `integer` scalar, equal to `1`, `2` or `3`. This argument is `intent(in)` and `optional`. Defines the type of DCT to be performed. The default type is `2`. @@ -1166,9 +1166,9 @@ Returns a `real` and rank-1 array, the DCT type-`t` of the input data `x`. #### Notes Within numerical accuracy, -- `x == idct(dct(x, t=1), t=1) / (2*(size(x) - 1))` -- `x == idct(dct(x, t=2), t=2) / (4*size(x))` -- `x == idct(dct(x, t=3), t=3) / (4*size(x))` +- `x == idct(dct(x, type=1), type=1) / (2*(size(x) - 1))` +- `x == idct(dct(x, type=2), type=2) / (4*size(x))` +- `x == idct(dct(x, type=3), type=3) / (4*size(x))` ### Example @@ -1177,9 +1177,9 @@ program demo_dct use fftpack, only: dct real(kind=8) :: x(4) = [1, 2, 3, 4] print *, dct(x,3,1) !! [8.0, -2.0, 0.0]. - print *, dct(x,t=1) !! [15.0, -4.0, 0.0, -1.0]. + print *, dct(x,type=1) !! [15.0, -4.0, 0.0, -1.0]. print *, dct(x,5,2) !! [14.36, -6.11, -5.0, 4.40, -2.65]. - print *, dct(dct(x,t=1),t=1)/(2*(4 - 1)) !! (normalized): [1.0, 2.0, 3.0, 4.0] + print *, dct(dct(x,type=1),type=1)/(2*(4 - 1)) !! (normalized): [1.0, 2.0, 3.0, 4.0] end program demo_dct ``` @@ -1201,7 +1201,7 @@ Pure function. ### Syntax -`result = [[fftpack(module):idct(interface)]](x [, n, t])` +`result = [[fftpack(module):idct(interface)]](x [, n, type])` ### Argument @@ -1213,7 +1213,7 @@ Transformed data to invert. This argument is `intent(in)` and `optional`. Defines the length of the Fourier transform. If `n` is not specified (the default) then `n = size(x)`. If `n <= size(x)`, `x` is truncated, if `n > size(x)`, `x` is zero-padded. -`t`: Shall be an `integer` scalar, equal to `1` or `2`. +`type`: Shall be an `integer` scalar, equal to `1` or `2`. This argument is `intent(in)` and `optional`. Defines the type of the IDCT to be performed. The default type is `2`. @@ -1224,9 +1224,9 @@ Returns a `real` and rank-1 array, the IDCT type-`t` of the input data `x`. #### Notes Within numerical accuracy, -- `x == idct(dct(x, t=1), t=1) / (2*(size(x) - 1))` -- `x == idct(dct(x, t=2), t=2) / (4*size(x))` -- `x == idct(dct(x, t=3), t=3) / (4*size(x))` +- `x == idct(dct(x, type=1), type=1) / (2*(size(x) - 1))` +- `x == idct(dct(x, type=2), type=2) / (4*size(x))` +- `x == idct(dct(x, type=3), type=3) / (4*size(x))` ### Example @@ -1234,9 +1234,9 @@ Within numerical accuracy, program demo_idct use fftpack, only: dct, idct real(kind=8) :: x(4) = [1, 2, 3, 4] - print *, idct(dct(x,t=1),t=1)/(2*(4-1)) !! (normalized): [1.0, 2.0, 3.0, 4.0] - print *, idct(dct(x,t=2),t=2)/(4*4) !! (normalized): [1.0, 2.0, 3.0, 4.0] - print *, idct(dct(x), n=3) !! (unnormalized): [22.06, 32.5, 65.65] + print *, idct(dct(x,type=1),type=1)/(2*(4-1)) !! (normalized): [1.0, 2.0, 3.0, 4.0] + print *, idct(dct(x,type=2),type=2)/(4*4) !! (normalized): [1.0, 2.0, 3.0, 4.0] + print *, idct(dct(x), n=3) !! (unnormalized): [22.06, 32.5, 65.65] end program demo_idct ``` diff --git a/src/fftpack.f90 b/src/fftpack.f90 index 3fff9d6..d02479a 100644 --- a/src/fftpack.f90 +++ b/src/fftpack.f90 @@ -247,10 +247,10 @@ end function irfft_rk !> Dsicrete cosine transforms. !> ([Specification](../page/specs/fftpack.html#dct)) interface dct - pure module function dct_rk(x, n, t) result(result) + pure module function dct_rk(x, n, type) result(result) real(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n - integer, intent(in), optional :: t + integer, intent(in), optional :: type real(kind=rk), allocatable :: result(:) end function dct_rk end interface dct @@ -260,10 +260,10 @@ end function dct_rk !> Inverse discrete cosine transforms. !> ([Specification](../page/specs/fftpack.html#idct)) interface idct - pure module function idct_rk(x, n, t) result(result) + pure module function idct_rk(x, n, type) result(result) real(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n - integer, intent(in), optional :: t + integer, intent(in), optional :: type real(kind=rk), allocatable :: result(:) end function idct_rk end interface idct diff --git a/src/fftpack_dct.f90 b/src/fftpack_dct.f90 index e8569ca..e765215 100644 --- a/src/fftpack_dct.f90 +++ b/src/fftpack_dct.f90 @@ -3,10 +3,10 @@ contains !> Discrete cosine transforms of types 1, 2, 3. - pure module function dct_rk(x, n, t) result(result) + pure module function dct_rk(x, n, type) result(result) real(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n - integer, intent(in), optional :: t + integer, intent(in), optional :: type real(kind=rk), allocatable :: result(:) integer :: lenseq, lensav, i @@ -25,7 +25,7 @@ pure module function dct_rk(x, n, t) result(result) end if ! Default to DCT-2 - if (.not.present(t)) then + if (.not.present(type)) then lensav = 3*lenseq + 15 allocate (wsave(lensav)) call dcosqi(lenseq, wsave) @@ -33,17 +33,17 @@ pure module function dct_rk(x, n, t) result(result) return end if - if (t == 1) then ! DCT-1 + if (type == 1) then ! DCT-1 lensav = 3*lenseq + 15 allocate (wsave(lensav)) call dcosti(lenseq, wsave) call dcost(lenseq, result, wsave) - else if (t == 2) then ! DCT-2 + else if (type == 2) then ! DCT-2 lensav = 3*lenseq + 15 allocate (wsave(lensav)) call dcosqi(lenseq, wsave) call dcosqb(lenseq, result, wsave) - else if (t == 3) then ! DCT-3 + else if (type == 3) then ! DCT-3 lensav = 3*lenseq + 15 allocate (wsave(lensav)) call dcosqi(lenseq, wsave) @@ -52,10 +52,10 @@ pure module function dct_rk(x, n, t) result(result) end function dct_rk !> Inverse discrete cosine transforms of types 1, 2, 3. - pure module function idct_rk(x, n, t) result(result) + pure module function idct_rk(x, n, type) result(result) real(kind=rk), intent(in) :: x(:) integer, intent(in), optional :: n - integer, intent(in), optional :: t + integer, intent(in), optional :: type real(kind=rk), allocatable :: result(:) integer :: lenseq, lensav, i @@ -74,7 +74,7 @@ pure module function idct_rk(x, n, t) result(result) end if ! Default to t=2; inverse DCT-2 is DCT-3 - if (.not.present(t)) then + if (.not.present(type)) then lensav = 3*lenseq + 15 allocate (wsave(lensav)) call dcosqi(lenseq, wsave) @@ -82,17 +82,17 @@ pure module function idct_rk(x, n, t) result(result) return end if - if (t == 1) then ! inverse DCT-1 is DCT-1 + if (type == 1) then ! inverse DCT-1 is DCT-1 lensav = 3*lenseq + 15 allocate (wsave(lensav)) call dcosti(lenseq, wsave) call dcost(lenseq, result, wsave) - else if (t == 2) then ! inverse DCT-2 is DCT-3 + else if (type == 2) then ! inverse DCT-2 is DCT-3 lensav = 3*lenseq + 15 allocate (wsave(lensav)) call dcosqi(lenseq, wsave) call dcosqf(lenseq, result, wsave) - else if (t == 3) then ! inverse DCT-3 is DCT-2 + else if (type == 3) then ! inverse DCT-3 is DCT-2 lensav = 3*lenseq + 15 allocate (wsave(lensav)) call dcosqi(lenseq, wsave) diff --git a/test/test_fftpack_dct.f90 b/test/test_fftpack_dct.f90 index 0256cd5..f680058 100644 --- a/test/test_fftpack_dct.f90 +++ b/test/test_fftpack_dct.f90 @@ -56,7 +56,7 @@ subroutine test_modernized_dct(error) ! DCT-1 call check(error, sum(abs(dct(x,2,1) - [0.0_rk, 18.0_rk])) < eps, "`dct(x,2,1)` failed.") if (allocated(error)) return - call check(error, sum(abs(dct(x,3,1) - dct(x,t=1))) < eps, "`dct(x,3,1)` failed.") + call check(error, sum(abs(dct(x,3,1) - dct(x,type=1))) < eps, "`dct(x,3,1)` failed.") if (allocated(error)) return call check(error, sum(abs(dct(x,4,1) - [real(kind=rk) :: -3, -3.0000000000000036_rk, 15, 33])) < eps,& "`dct(x,4,1)` failed.") @@ -74,12 +74,12 @@ subroutine test_modernized_dct(error) call check(error, sum(abs(dct(x,2,3) - [-3.7279220613578570_rk, 21.727922061357859_rk])) < eps, & "`dct(x,2,3)` failed.") if (allocated(error)) return - call check(error, sum(abs(dct(x,3,3) - dct(x,t=3))) < eps, & + call check(error, sum(abs(dct(x,3,3) - dct(x,type=3))) < eps, & "`dct(x,3,3)` failed.") if (allocated(error)) return call check(error, sum(abs(dct(x,4,3) - [-3.3871908980838743_rk, -2.1309424696909023_rk, & 11.645661095452331_rk, 29.872472272322447_rk])) < eps, & - "`dct(x,n=4,t=3)` failed.") + "`dct(x,4,3)` failed.") end subroutine test_modernized_dct @@ -89,36 +89,40 @@ subroutine test_modernized_idct(error) real(kind=rk) :: x(4) = [1, 2, 3, 4] ! inverse DCT-1 - call check(error, sum(abs(idct(dct(x,t=1),t=1)/(2*3) - x)) < eps, "`idct(dct(x,t=1),t=1)/(2*3)` failed.") + call check(error, sum(abs(idct(dct(x,type=1),type=1)/(2*3) - x)) < eps, & + "`idct(dct(x,type=1),type=1)/(2*3)` failed.") if (allocated(error)) return - call check(error, sum(abs(idct(dct(x,t=1), 2, 1)/(2*1) - [5.5_rk, 9.5_rk])) < eps,& - "`idct(dct(x,t=1), 2, 1)/(2*1)` failed.") + call check(error, sum(abs(idct(dct(x,type=1), 2, 1)/(2*1) - [5.5_rk, 9.5_rk])) < eps, & + "`idct(dct(x,type=1), 2, 1)/(2*1)` failed.") if (allocated(error)) return - call check(error, sum(abs(idct(dct(x,2,1), 4, 1)/(2*3) - [0.16666666666666666_rk, 0.33333333333333331_rk,& - 0.66666666666666663_rk, 0.83333333333333315_rk])) < eps,& + call check(error, sum(abs(idct(dct(x,2,1), 4, 1)/(2*3) - [0.16666666666666666_rk, 0.33333333333333331_rk, & + 0.66666666666666663_rk, 0.83333333333333315_rk])) < eps, & "`idct(dct(x,2,1), 4, 1)/(2*3)` failed.") - + ! inverse DCT-2 x = [1, 2, 3, 4] - call check(error, sum(abs(idct(dct(x,t=2))/(4*4) - x)) < eps, & - "`idct(dct(x, t=2))/(4*4)` failed.") + call check(error, sum(abs(idct(dct(x,type=2))/(4*4) - x)) < eps, & + "`idct(dct(x, type=2))/(4*4)` failed.") if (allocated(error)) return - call check(error, sum(abs(idct(dct(x,t=2),n=2) - [22.156460020898692_rk, 57.843539979101308_rk])) < eps,& - "`idct(dct(x,t=2),n=2)` failed.") + call check(error, sum(abs(idct(dct(x,type=2),n=2) - [22.156460020898692_rk, 57.843539979101308_rk])) < eps, & + "`idct(dct(x, type=2),n=2)` failed.") if (allocated(error)) return - call check(error, sum(abs(idct(dct(x,n=2,t=2),n=4) - [6.7737481404944937_rk, 9.8352155994152106_rk,& - 14.164784400584789_rk, 17.226251859505506_rk])) < eps, "`idct(dct(x,n=2,t=2),n=4)` failed.") + call check(error, sum(abs(idct(dct(x,n=2,type=2),n=4) - [6.7737481404944937_rk, 9.8352155994152106_rk, & + 14.164784400584789_rk, 17.226251859505506_rk])) < eps, & + "`idct(dct(x, n=2, type=2), n=4)` failed.") ! inverse DCT-3 x = [1, 2, 3, 4] - call check(error, sum(abs(idct(dct(x,t=3),t=3)/(4*4) - x)) < eps, & - "`idct(dct(x, t=3), t=3)/(4*4)` failed.") + call check(error, sum(abs(idct(dct(x,type=3),type=3)/(4*4) - x)) < eps, & + "`idct(dct(x, type=3), type=3)/(4*4)` failed.") if (allocated(error)) return - call check(error, sum(abs(idct(dct(x,t=3),n=2,t=3)/(4*2) - [1.4483415291679655_rk, 7.4608849947753271_rk])) < eps, & - "`idct(dct(x,t=3),n=2,t=3)/(4*2)` failed.") + call check(error, sum(abs(idct(dct(x,type=3),n=2,type=3)/(4*2) - & + [1.4483415291679655_rk, 7.4608849947753271_rk])) < eps, & + "`idct(dct(x, type=3), n=2, type=3)/(4*2)` failed.") if (allocated(error)) return - call check(error, sum(abs(idct(dct(x,n=2,t=3),n=4,t=3)/(4*4) - [0.5_rk, 0.70932417358418376_rk, 1.0_rk, & - 0.78858050747473762_rk])) < eps, "`idct(dct(x,n=2,t=3),n=4, t=3)/(4*4)` failed.") + call check(error, sum(abs(idct(dct(x,n=2,type=3),n=4,type=3)/(4*4) - & + [0.5_rk, 0.70932417358418376_rk, 1.0_rk, 0.78858050747473762_rk])) < eps, & + "`idct(dct(x, n=2, type=3), n=4, type=3)/(4*4)` failed.") end subroutine test_modernized_idct From 62396ea7c52203d9be8994feda49a894a5a220ee Mon Sep 17 00:00:00 2001 From: Chris Val Date: Fri, 3 Mar 2023 15:10:32 -0500 Subject: [PATCH 5/5] docs reference link --- doc/specs/fftpack.md | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/doc/specs/fftpack.md b/doc/specs/fftpack.md index 1fd471f..be8b73f 100644 --- a/doc/specs/fftpack.md +++ b/doc/specs/fftpack.md @@ -1159,11 +1159,11 @@ Defines the length of the DCT. If `n` is not specified (the default) then `n = s This argument is `intent(in)` and `optional`. Defines the type of DCT to be performed. The default type is `2`. -#### Return value +### Return value Returns a `real` and rank-1 array, the DCT type-`t` of the input data `x`. -#### Notes +### Notes Within numerical accuracy, - `x == idct(dct(x, type=1), type=1) / (2*(size(x) - 1))` @@ -1217,11 +1217,11 @@ Defines the length of the Fourier transform. If `n` is not specified (the defaul This argument is `intent(in)` and `optional`. Defines the type of the IDCT to be performed. The default type is `2`. -#### Return value +### Return value Returns a `real` and rank-1 array, the IDCT type-`t` of the input data `x`. -#### Notes +### Notes Within numerical accuracy, - `x == idct(dct(x, type=1), type=1) / (2*(size(x) - 1))` @@ -1240,6 +1240,10 @@ program demo_idct end program demo_idct ``` +## References + +[1] Wikipedia, "Discrete cosine transform", [https://en.wikipedia.org/wiki/Discrete_cosine_transform](https://en.wikipedia.org/wiki/Discrete_cosine_transform) + # Utility functions ## `fftshift`