zfftiInitializes 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.
Experimental.
Pure subroutine.
call zffti(n, wsave)
n: Shall be an integer scalar.
This argument is intent(in).
The length of the sequence to be transformed.
wsave: Shall be a real array.
This argument is intent(out).
A work array which must be dimensioned at least 4*n+15.
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.
The contents of wsave must not be changed between calls of zfftf or zfftb.
program demo_zffti
use fftpack, only: zffti
complex(kind=8) :: x(4) = [1.0, 2.0, 3.0, 4.0]
real(kind=8) :: w(31)
call zffti(4,w)
end program demo_zffti
zfftfComputes the forward complex discrete fourier transform (the fourier analysis).
Equivalently, zfftf computes the fourier coefficients of a complex periodic sequence.
The transform is defined below at output parameter c.
The transform is not normalized. To obtain a normalized transform the output must be divided by n. Otherwise a call of zfftf followed by a call of zfftb will multiply the sequence by n.
The array wsave which is used by subroutine zfftf must be initialized by calling subroutine zffti(n,wsave).
Experimental.
Pure subroutine.
call zfftf(n, c, wsave)
n: Shall be an integer scalar.
This argument is intent(in).
The length of the complex sequence c. The method is more efficient when n is the product of small primes.
c: Shall be a complex and rank-1 array.
This argument is intent(inout).
A complex array of length n which contains the sequence.
for j=1,...,n
c(j)=the sum from k=1,...,n of
c(k)*exp(-i*(j-1)*(k-1)*2*pi/n)
where i=sqrt(-1)
wsave: Shall be a real array.
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.
The contents of wsave must not be changed between calls of zfftf or zfftb.
program demo_zfftf
use fftpack, only: zffti, zfftf
complex(kind=8) :: x(4)
real(kind=8) :: w(31)
x = [real(kind=8) :: 1.0, 2.0, 3.0, 4.0]
call zffti(4,w)
call zfftf(4,x,w) !! `x` returns [(10.0,0.0), (-2.0,2.0), (-2.0,0.0), (-2.0,-2.0)].
end program demo_zfftf
zfftbUnnormalized inverse of zfftf.
Computes the backward complex discrete fourier transform (the fourier synthesis).
Equivalently, zfftb computes a complex periodic sequence from its fourier coefficients.
The transform is defined below at output parameter c.
The transform is not normalized. to obtain a normalized transform the output must be divided by n. Otherwise a call of zfftf followed by a call of zfftb will multiply the sequence by n.
The array wsave which is used by subroutine zfftf must be initialized by calling subroutine zffti(n,wsave).
Experimental.
Pure subroutine.
call zfftb(n, c, wsave)
n: Shall be an integer scalar.
This argument is intent(in).
The length of the complex sequence c. The method is more efficient when n is the product of small primes.
c: Shall be a complex array.
This argument is intent(inout).
A complex array of length n which contains the sequence.
for j=1,...,n
c(j)=the sum from k=1,...,n of
c(k)*exp(-i*(j-1)*(k-1)*2*pi/n)
where i=sqrt(-1)
wsave: Shall be a real array.
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.
The contents of wsave must not be changed between calls of zfftf or zfftb.
program demo_zfftb
use fftpack, only: zffti, zfftf, zfftb
complex(kind=8) :: x(4)
real(kind=8) :: w(31)
x = [real(kind=8) :: 1.0, 2.0, 3.0, 4.0]
call zffti(4,w)
call zfftf(4,x,w) !! `x` returns [(10.0,0.0), (-2.0,2.0), (-2.0,0.0), (-2.0,-2.0)].
call zfftb(4,x,w) !! `x` returns [(4.0,0.0), (8.0,0.0), (12.0,0.0), (16.0,0.0)].
end program demo_zfftb
fftComputes the forward complex discrete fourier transform (the fourier analysis).
Experimental.
Pure function.
result = fft(x [, n])
x: Shall be a complex and rank-1 array.
This argument is intent(in).
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.
Returns a complex and rank-1 array, the Discrete Fourier Transform (DFT) of x.
Within numerical accuracy, x == ifft(fft(x))/size(x).
program demo_fft
use fftpack, only: fft
complex(kind=8) :: x(4)
x = [real(kind=8) :: 1.0, 2.0, 3.0, 4.0]
print *, fft(x) !! [(10.0,0.0), (-2.0,2.0), (-2.0,0.0), (-2.0,-2.0)].
print *, fft(x,3) !! [(6.0,0.0), (-1.5,0.86), (-1.5,0.86)].
print *, fft(x,5) !! [(10.0,0.0), (-4.0,1.3), (1.5,-2.1), (1.5,2.1), (-4.0,1.3)].
end program demo_fft
ifftUnnormalized inverse of fft.
Experimental.
Pure function.
result = ifft(x [, n])
x: Shall be a complex and rank-1 array.
This argument is intent(in).
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.
Returns a complex and rank-1 array, the unnormalized inverse Discrete Fourier Transform (DFT) of x.
program demo_ifft
use fftpack, only: fft, ifft
complex(kind=8) :: x(4) = [1.0, 2.0, 3.0, 4.0]
print *, ifft(fft(x))/4.0 !! [(1.0,0.0), (2.0,0.0), (3.0,0.0), (4.0,0.0)]
print *, ifft(fft(x), 3) !! [(6.0,2.0), (10.3,-1.0), (13.73,-1.0)]
end program demo_ifft
dfftiInitializes 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.
Experimental.
Pure subroutine.
call dffti(n, wsave)
n: Shall be an integer scalar.
This argument is intent(in).
The length of the sequence to be transformed.
wsave: Shall be a real array.
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.
The contents of wsave must not be changed between calls of dfftf or dfftb.
program demo_dffti
use fftpack, only: dffti
real(kind=8) :: x(4) = [1.0, 2.0, 3.0, 4.0]
real(kind=8) :: w(23)
call dffti(4,w)
end program demo_dffti
dfftfComputes the fourier coefficients of a real perodic sequence (fourier analysis).
The transform is defined below at output parameter r.
The transform is not normalized. To obtain a normalized transform the output must be divided by n. Otherwise a call of dfftf followed by a call of dfftb will multiply the sequence by n.
The array wsave which is used by subroutine dfftf must be initialized by calling subroutine dffti(n,wsave).
Experimental.
Pure subroutine.
call dfftf(n, r, wsave)
n: Shall be an integer scalar.
This argument is intent(in).
The length of the real sequence r. The method is more efficient when n is the product of small primes.
n may change so long as different work arrays are provided.
r: Shall be a real array.
This argument is intent(inout).
A real array of length n which contains the sequence.
r(1) = the sum from i=1 to i=n of r(i)
if n is even set l =n/2 , if n is odd set l = (n+1)/2
then for k = 2,...,l
r(2*k-2) = the sum from i = 1 to i = n of
r(i)*cos((k-1)*(i-1)*2*pi/n)
r(2*k-1) = the sum from i = 1 to i = n of
-r(i)*sin((k-1)*(i-1)*2*pi/n)
if n is even
r(n) = the sum from i = 1 to i = n of
(-1)**(i-1)*r(i)
wsave: Shall be a real array.
This argument is intent(in).
A real work array which must be dimensioned at least 4n+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.
The contents of wsave must not be changed between calls of dfftf or dfftb.
program demo_dfftf
use fftpack, only: dffti, dfftf
real(kind=8) :: x(4) = [1, 2, 3, 4]
real(kind=8) :: w(23)
call dffti(4,w)
call dfftf(4,x,w) !! `x` returns [10.0, -2.0, 2.0, -2.0].
end program demo_dfftf
dfftbUnnormalized inverse of dfftf.
Computes the backward real discrete fourier transform (the fourier synthesis).
Equivalently, dfftb computes a real periodic sequence from its fourier coefficients.
The transform is defined below at output parameter c.
The transform is not normalized. To obtain a normalized transform the output must be divided by n. Otherwise a call of dfftf followed by a call of dfftb will multiply the sequence by n.
The array wsave which is used by subroutine dfftf must be initialized by calling subroutine dffti(n,wsave).
Experimental.
Pure subroutine.
call dfftb(n, r, wsave)
n: Shall be an integer scalar.
This argument is intent(in).
The length of the real sequence r. The method is more efficient when n is the product of small primes.
r: Shall be a real array.
This argument is intent(inout).
A real array of length n which contains the sequence.
for n even and for i = 1,...,n
r(i) = r(1)+(-1)**(i-1)*r(n)
plus the sum from k=2 to k=n/2 of
2.*r(2*k-2)*cos((k-1)*(i-1)*2*pi/n)
-2.*r(2*k-1)*sin((k-1)*(i-1)*2*pi/n)
for n odd and for i = 1,...,n
r(i) = r(1) plus the sum from k=2 to k=(n+1)/2 of
2.*r(2*k-2)*cos((k-1)*(i-1)*2*pi/n)
-2.*r(2*k-1)*sin((k-1)*(i-1)*2*pi/n)
wsave: Shall be a real array.
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.
The contents of wsave must not be changed between calls of dfftf or dfftb.
program demo_dfftb
use fftpack, only: dffti, dfftf, dfftb
real(kind=8) :: x(4) = [1, 2, 3, 4]
real(kind=8) :: w(31)
call dffti(4,w)
call dfftf(4,x,w) !! `x` returns [10.0, -2.0, 2.0, -2.0].
call dfftb(4,x,w) !! `x` returns [4.0, 8.0, 12.0, 16.0], which is not normalized.
end program demo_dfftb
rfftDiscrete Fourier transform of a real sequence.
Experimental.
Pure function.
result = rfft(x [, n])
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.
Returns a real and rank-1 array, the Discrete Fourier Transform (DFT) of x.
Within numerical accuracy, y == rfft(irfft(y))/size(y).
program demo_rfft
use fftpack, only: rfft
real(kind=8) :: x(4) = [1, 2, 3, 4]
print *, rfft(x,3) !! [6.0, -1.5, 0.87].
print *, rfft(x) !! [10.0, -2.0, 2.0, -2.0].
print *, rfft(x,5) !! [10.0, -4.0, -1.3, 1.5, -2.1].
end program demo_rfft
irfftUnnormalized inverse of rfft.
Experimental.
Pure function.
result = irfft(x [, n])
x: Shall be a real array.
This argument is intent(in).
Transformed data to invert.
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.
Returns a real and rank-1 array, the unnormalized inverse discrete Fourier transform.
program demo_irfft
use fftpack, only: rfft, irfft
real(kind=8) :: x(4) = [1, 2, 3, 4]
print *, irfft(rfft(x))/4.0 !! [1.0, 2.0, 3.0, 4.0]
print *, irfft(rfft(x), 3) !! [6.0, 8.53, 15.46]
end program demo_irfft
dzfftiInitializes 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.
Experimental
Prue function.
call dzffti(n, wsave)
n: Shall be an integer scalar.
This argument is intent(in).
The length of the sequence to be transformed.
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 dzfftf and dzfftb 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 dzfftf or dzfftb.
program demo_dzffti
use fftpack, only: dzffti
real(kind=8) :: x(4) = [1, 2, 3, 4]
real(kind=8) :: w(3*4 + 15)
call dzffti(4, w) !! Initializes the array `w` which is used in both `dzfftf` and `dzfftb`.
end program demo_dzffti
dzfftfComputes 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.
Experimental
Pure subroutine.
call dzfftf(n, r, azero, a, b, wsave)
n: Shall be an integer scalar.
This argument is intent(in).
The length of the array r to be transformed.
The method is most efficient when n is the product of small primes.
r: Shall be a real and rank-1 array.
This argument is intent(in).
A real array of length n which contains the sequence to be transformed. r is not destroyed.
azero: Shall be a real scalar.
This argument is intent(out).
The sum from i=1 to i=n of r(i)/n.
a, b: Shall be a real and rank-1 array.
This argument is intent(out).
for n even b(n/2)=0. and a(n/2) is the sum from i=1 to i=n of (-1)**(i-1)*r(i)/n
for n even define kmax=n/2-1
for n odd define kmax=(n-1)/2
then for k=1,...,kmax
a(k) equals the sum from i=1 to i=n of
2./n*r(i)*cos(k*(i-1)*2*pi/n)
b(k) equals the sum from i=1 to i=n of
2./n*r(i)*sin(k*(i-1)*2*pi/n)
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 dzfftf. The wsave array must be initialized by calling subroutine dzffti(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 dzfftf and dzfftb.
program demo_dzfftf
use fftpack, only: dzffti, dzfftf
real(kind=8) :: x(4) = [1, 2, 3, 4]
real(kind=8) :: w(3*4 + 15)
real(kind=8) :: azero, a(4/2), b(4/2)
call dzffti(4, w)
call dzfftf(4, x, azero, a, b, w) !! `azero`: 2.5; `a`: [-1.0, -0.5]; `b`: [-1.0, -0.0]
end program demo_dzfftf
dzfftbComputes 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.
Experimental
Pure subroutine.
call dzfftb(n, r, azero, a, b, wsave)
n: Shall be an integer scalar.
This argument is intent(in).
The length of the output array r.
The method is most efficient when n is the product of small primes.
r: Shall be a real and rank-1 array.
This argument is intent(out).
if n is even define kmax=n/2
if n is odd define kmax=(n-1)/2
then for i=1,...,n
r(i)=azero plus the sum from k=1 to k=kmax of
a(k)*cos(k*(i-1)*2*pi/n)+b(k)*sin(k*(i-1)*2*pi/n)
Complex notation:
for j=1,...,n
r(j) equals the sum from k=-kmax to k=kmax of
c(k)*exp(i*k*(j-1)*2*pi/n)
where
c(k) = .5*cmplx(a(k),-b(k)) for k=1,...,kmax
c(-k) = conjg(c(k))
c(0) = azero
and i=sqrt(-1)
Amplitude - phase notation:
for i=1,...,n
r(i) equals azero plus the sum from k=1 to k=kmax of
alpha(k)*cos(k*(i-1)*2*pi/n+beta(k))
where
alpha(k) = sqrt(a(k)*a(k)+b(k)*b(k))
cos(beta(k))=a(k)/alpha(k)
sin(beta(k))=-b(k)/alpha(k)
azero: Shall be a real scalar.
This argument is intent(in).
The constant fourier coefficient.
a, b: Shall be a real and rank-1 array.
This argument is intent(in).
Arrays which contain the remaining fourier coefficients these arrays are not destroyed.
The length of these arrays depends on whether n is even or odd.
if n is even n/2 locations are required
if n is odd (n-1)/2 locations are required
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 dzfftf. The wsave array must be initialized by calling subroutine dzffti(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 dzfftf and dzfftb.
program demo_dzfftb
use fftpack, only: dzffti, dzfftf, dzfftb
real(kind=8) :: x(4) = [1, 2, 3, 4]
real(kind=8) :: w(3*4 + 15)
real(kind=8) :: azero, a(4/2), b(4/2)
call dzffti(4, w)
call dzfftf(4, x, azero, a, b, w) !! `azero`: 2.5; `a`: [-1.0, -0.5]; `b`: [-1.0, -0.0]
x = 0.0
call dzfftb(4, x, azero, a, b, w) !! `x`: [1.0, 2.0, 3.0, 4.0]
end program demo_dzfftb
dcosti or dct_t1iInitializes 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.
Experimental
Pure subroutine.
call dcosti(n , wsave)
n: Shall be a integer scalar.
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.
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.
program demo_dcosti
use fftpack, only: dcosti
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
dcost or dct_t1Computes 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 ) is
equivalent to the DFT of data abcdedcb (size ).
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).
Experimental
Pure subroutine.
call dcost(n, x, wsave)
n: Shall be a 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.
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
x(i) = x(1)+(-1)**(i-1)*x(n)
+ the sum from k=2 to k=n-1
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 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.
program demo_dcost
use fftpack, only: dcosti, dcost
real(kind=8) :: x(4) = [1, 2, 3, 4]
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
dcosqi or dct_t23iInitializes 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.
The two procedures are completely equivalent and expect the same arguments. It is a matter of personal preference which one you choose to use.
Experimental
Pure subroutine.
call dcosqi(n, wsave)
n: Shall be an 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.
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.
program demo_dcosqi
use fftpack, only: dcosqi
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
dcosqf or dct_t3Computes 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.
The array wsave which is used by subroutine dcosqf must be initialized by calling subroutine dcosqi(n,wsave).
Experimental
Pure subroutine.
call dcosqf(n, x, wsave)
n: Shall be an 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.
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
x(i) = x(1) plus the sum from k=2 to k=n of
2*x(k)*cos((2*i-1)*(k-1)*pi/(2*n))
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.
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 of the same n.
program demo_dcosqf
use fftpack, only: dcosqi, dcosqf
real(kind=8) :: w(3*4 + 15)
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
dcosqb or dct_t2Computes 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.
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).
Experimental
Pure subroutine.
call dcosqb(n, x, wsave)
n: Shall be an 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.
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
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 of the same n.
program demo_dcosqb
use fftpack, only: dcosqi, dcosqf, dcosqb
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) !! [64.0, 48.0, 80.0, 160.0]
end program demo_dcosqb
dctDiscrete 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.
Experimental.
Pure function.
result = dct(x [, n, type])
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 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.
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.
Returns a real and rank-1 array, the DCT type-t of the input data x.
Within numerical accuracy,
- 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))
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,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,type=1),type=1)/(2*(4 - 1)) !! (normalized): [1.0, 2.0, 3.0, 4.0]
end program demo_dct
idctUnnormalized 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.
Experimental.
Pure function.
result = idct(x [, n, type])
x: Shall be a real array.
This argument is intent(in).
Transformed data to invert.
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.
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.
Returns a real and rank-1 array, the IDCT type-t of the input data x.
Within numerical accuracy,
- 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))
program demo_idct
use fftpack, only: dct, idct
real(kind=8) :: x(4) = [1, 2, 3, 4]
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
[1] Wikipedia, "Discrete cosine transform", https://en.wikipedia.org/wiki/Discrete_cosine_transform
fftshiftRearranges the Fourier transform by moving the zero-frequency component to the center of the array.
Experimental.
Pure function.
result = fftshift(x)
x: Shall be a complex/real and rank-1 array.
This argument is intent(in).
Returns the complex/real and rank-1 Fourier transform by moving the zero-frequency component to the center of the array.
program demo_fftshift
use fftpack, only: fftshift
complex(kind=8) :: c(5) = [1, 2, 3, 4, 5]
real(kind=8) :: x(5) = [1, 2, 3, 4, 5]
print *, fftshift(c(1:4)) !! [(3.0,0.0), (4.0,0.0), (1.0,0.0), (2.0,0.0)]
print *, fftshift(c) !! [(4.0,0.0), (5.0,0.0), (1.0,0.0), (2.0,0.0), (3.0,0.0)]
print *, fftshift(x(1:4)) !! [3.0, 4.0, 1.0, 2.0]
print *, fftshift(x) !! [4.0, 5.0, 1.0, 2.0, 3.0]
end program demo_fftshift
ifftshiftRearranges the Fourier transform with zero frequency shifting back to the original transform output. In other words, ifftshift is the result of undoing fftshift.
Experimental.
Pure function.
result = ifftshift(x)
x: Shall be a complex/real and rank-1 array.
This argument is intent(in).
Returns the complex/real and rank-1 Fourier transform with zero frequency shifting back to the original transform output.
program demo_ifftshift
use fftpack, only: fftshift, ifftshift
complex(kind=8) :: c(5) = [1, 2, 3, 4, 5]
real(kind=8) :: x(5) = [1, 2, 3, 4, 5]
print *, ifftshift(fftshift(c(1:4))) !! [(1.0,0.0), (2.0,0.0), (3.0,0.0), (4.0,0.0)]
print *, ifftshift(fftshift(c) ) !! [(1.0,0.0), (2.0,0.0), (3.0,0.0), (4.0,0.0), (5.0,0.0)]
print *, ifftshift(fftshift(x(1:4))) !! [1.0, 2.0, 3.0, 4.0]
print *, ifftshift(fftshift(x)) !! [1.0, 2.0, 3.0, 4.0, 5.0]
end program demo_ifftshift
fftfreqReturns 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).
Experimental.
Pure function.
result = fftfreq(n)
n: Shall be an integer, equal to the length of the corresponding complex discrete Fourier transform.
This argument is intent(in).
Returns the integer and rank-1 array of the transform's frequency values in the standard order (zero frequency first).
program demo_fftfreq
use fftpack, only: fftfreq
print *, fftfreq(4) ! [0, 1, -2, -1]
print *, fftfreq(5) ! [0, 1, 2, -2, -1]
end program demo_fftfreq
rfftfreqReturns 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).
Experimental.
Pure function.
result = rfftfreq(n)
n: Shall be an integer, equal to the length of the corresponding real discrete Fourier transform.
This argument is intent(in).
Returns the integer and rank-1 array of the transform's frequency values in the standard order (zero frequency first).
program demo_rfftfreq
use fftpack, only: rfftfreq
print *, rfftfreq(4) ! [0, 1, 1, -2]
print *, rfftfreq(5) ! [0, 1, 1, 2, 2]
end program demo_rfftfreq