Hi all,
I am trying to calculate derivatives of chebyshev coefficients using FFTW
for cosine(x). However I am not getting back the right values when I
perform inverse chebyshev transform. Can anyone check the code below and
tell me the possible errors. Thank you.
program test
implicit none
#include "fftw3.f"
INTEGER, PARAMETER :: N=4
INTEGER(8) :: plan
REAL(8), DIMENSION(N) :: in, out, out1, in1
REAL(8), DIMENSION(N) :: dout
REAL :: pi, y
INTEGER :: j, p
pi = 3.14
write(*,*) 'Input array:'
do j = 1,N,1
y =-1*cos(pi*(j)/N)
in(j) = cos(2*pi*y/N)
enddo
call dfftw_plan_r2r_1d ( plan, N, in, out,
^ FFTW_REDFT00, FFTW_ESTIMATE )
call dfftw_execute ( plan )
call dfftw_destroy_plan ( plan )
out = out/6
do j = 1,N,1
write(*,*) ' out(',j,') = ',out(j)
enddo
dout(N) = 0
dout(N-1) = 2*(N)*out(N)
do j = N-2,1,-1
dout(j) = dout(j+2) + 2*(j+1)*out(j+1)
enddo
call dfftw_plan_r2r_1d(plan, N+1, dout, out1,
^ FFTW_REDFT00, FFTW_ESTIMATE)
call dfftw_execute ( plan )
call dfftw_destroy_plan ( plan )
write(*,*) 'Output array after inverse FFT:'
do j = 1,N,1
write(*,*) ' out1(',j,') = ',out1(j)
enddo
do j = 1,N,1
y =-1*cos(pi*(j)/N)
in1(j) = -1*sin(2*pi*y/N)
write(*,*) ' in1(',j,') = ',in1(j)
enddo
end
Regards,
Sasidhar
DCT derivative
Started by ●November 21, 2008
Reply by ●November 21, 20082008-11-21
On Nov 21, 11:42�am, "sasidhar" <sasidha...@gmail.com> wrote:> Hi all, > > I am trying to calculate derivatives of chebyshev coefficients using FFTW > for cosine(x). However I am not getting back the right values when I > perform inverse chebyshev transform. Can anyone check the code below and > tell me the possible errors. Thank you. > > � � � program test > � � � implicit none > #include "fftw3.f" > > � � � INTEGER, PARAMETER �:: �N=4 > � � � INTEGER(8) :: �plan > � � � REAL(8), DIMENSION(N) :: �in, out, out1, in1 > � � � REAL(8), DIMENSION(N) :: dout > � � � REAL :: �pi, y > � � � INTEGER :: j, p > > � � � pi = 3.14 > � � � write(*,*) 'Input array:' > � � � do j = 1,N,1 > � � � � y =-1*cos(pi*(j)/N) > � � � � in(j) = cos(2*pi*y/N) > � � � enddo > > � � � call dfftw_plan_r2r_1d ( plan, N, in, out, > � � �^ � � � � � � � � � � � � FFTW_REDFT00, FFTW_ESTIMATE ) > � � � call dfftw_execute ( plan ) > � � � call dfftw_destroy_plan ( plan ) > � � � out = out/6 > � � � do j = 1,N,1 > � � � � write(*,*) ' � �out(',j,') = ',out(j) > � � � enddo > > � � � dout(N) = 0 > � � � dout(N-1) = 2*(N)*out(N) > > � � � do j = N-2,1,-1 > � � � � dout(j) = dout(j+2) + 2*(j+1)*out(j+1) > � � � enddo > > � � � call dfftw_plan_r2r_1d(plan, N+1, dout, out1, > � � �^ � � � � � � � � � � � � FFTW_REDFT00, FFTW_ESTIMATE) > � � � call dfftw_execute ( plan ) > � � � call dfftw_destroy_plan ( plan ) > > � � � write(*,*) 'Output array after inverse FFT:' > � � � do j = 1,N,1 > � � � � write(*,*) ' � �out1(',j,') = ',out1(j) > � � � enddo > > � � � do j = 1,N,1 > � � � � y =-1*cos(pi*(j)/N) > � � � � in1(j) = -1*sin(2*pi*y/N) > � � � � write(*,*) ' � �in1(',j,') = ',in1(j) > � � � enddo > � � � end > > �Regards, > �SasidharHello Sasidhar, I don't know about your code, but derivatives of coefficients (constants?) are zero. Or are you trying to find derivatives of Chebyshev polynomials? These do have well defined forms in terms Chebyshev polynomials (although both 1st and 2nd kind are needed though). Or you can just stay with polynomials. Maybe if you can give a clearer description of what you are really trying to do, then we may be better able to help. Are you trying to do differentiation via FFTs? You know aliasing can be an issue here! Clay
Reply by ●November 21, 20082008-11-21
clay@claysturner.com wrote:> On Nov 21, 11:42 am, "sasidhar" <sasidha...@gmail.com> wrote: >> Hi all, >> >> I am trying to calculate derivatives of chebyshev coefficients using FFTW >> for cosine(x). However I am not getting back the right values when I >> perform inverse chebyshev transform. Can anyone check the code below and >> tell me the possible errors. Thank you. >> >> program test >> implicit none >> #include "fftw3.f" >> >> INTEGER, PARAMETER :: N=4 >> INTEGER(8) :: plan >> REAL(8), DIMENSION(N) :: in, out, out1, in1 >> REAL(8), DIMENSION(N) :: dout >> REAL :: pi, y >> INTEGER :: j, p >> >> pi = 3.14 >> write(*,*) 'Input array:' >> do j = 1,N,1 >> y =-1*cos(pi*(j)/N) >> in(j) = cos(2*pi*y/N) >> enddo >> >> call dfftw_plan_r2r_1d ( plan, N, in, out, >> ^ FFTW_REDFT00, FFTW_ESTIMATE ) >> call dfftw_execute ( plan ) >> call dfftw_destroy_plan ( plan ) >> out = out/6 >> do j = 1,N,1 >> write(*,*) ' out(',j,') = ',out(j) >> enddo >> >> dout(N) = 0 >> dout(N-1) = 2*(N)*out(N) >> >> do j = N-2,1,-1 >> dout(j) = dout(j+2) + 2*(j+1)*out(j+1) >> enddo >> >> call dfftw_plan_r2r_1d(plan, N+1, dout, out1, >> ^ FFTW_REDFT00, FFTW_ESTIMATE) >> call dfftw_execute ( plan ) >> call dfftw_destroy_plan ( plan ) >> >> write(*,*) 'Output array after inverse FFT:' >> do j = 1,N,1 >> write(*,*) ' out1(',j,') = ',out1(j) >> enddo >> >> do j = 1,N,1 >> y =-1*cos(pi*(j)/N) >> in1(j) = -1*sin(2*pi*y/N) >> write(*,*) ' in1(',j,') = ',in1(j) >> enddo >> end >> >> Regards, >> Sasidhar > > Hello Sasidhar, > > I don't know about your code, but derivatives of coefficients > (constants?) are zero. Or are you trying to find derivatives of > Chebyshev polynomials? These do have well defined forms in terms > Chebyshev polynomials (although both 1st and 2nd kind are needed > though). Or you can just stay with polynomials. > > Maybe if you can give a clearer description of what you are really > trying to do, then we may be better able to help. > > Are you trying to do differentiation via FFTs? You know aliasing can > be an issue here!Clay, What is a Chebychev transform? Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Reply by ●November 21, 20082008-11-21
On Nov 21, 2:52�pm, Jerry Avins <j...@ieee.org> wrote:> c...@claysturner.com wrote: > > On Nov 21, 11:42 am, "sasidhar" <sasidha...@gmail.com> wrote: > >> Hi all, > > >> I am trying to calculate derivatives of chebyshev coefficients using FFTW > >> for cosine(x). However I am not getting back the right values when I > >> perform inverse chebyshev transform. Can anyone check the code below and > >> tell me the possible errors. Thank you. > > >> � � � program test > >> � � � implicit none > >> #include "fftw3.f" > > >> � � � INTEGER, PARAMETER �:: �N=4 > >> � � � INTEGER(8) :: �plan > >> � � � REAL(8), DIMENSION(N) :: �in, out, out1, in1 > >> � � � REAL(8), DIMENSION(N) :: dout > >> � � � REAL :: �pi, y > >> � � � INTEGER :: j, p > > >> � � � pi = 3.14 > >> � � � write(*,*) 'Input array:' > >> � � � do j = 1,N,1 > >> � � � � y =-1*cos(pi*(j)/N) > >> � � � � in(j) = cos(2*pi*y/N) > >> � � � enddo > > >> � � � call dfftw_plan_r2r_1d ( plan, N, in, out, > >> � � �^ � � � � � � � � � � � � FFTW_REDFT00, FFTW_ESTIMATE ) > >> � � � call dfftw_execute ( plan ) > >> � � � call dfftw_destroy_plan ( plan ) > >> � � � out = out/6 > >> � � � do j = 1,N,1 > >> � � � � write(*,*) ' � �out(',j,') = ',out(j) > >> � � � enddo > > >> � � � dout(N) = 0 > >> � � � dout(N-1) = 2*(N)*out(N) > > >> � � � do j = N-2,1,-1 > >> � � � � dout(j) = dout(j+2) + 2*(j+1)*out(j+1) > >> � � � enddo > > >> � � � call dfftw_plan_r2r_1d(plan, N+1, dout, out1, > >> � � �^ � � � � � � � � � � � � FFTW_REDFT00, FFTW_ESTIMATE) > >> � � � call dfftw_execute ( plan ) > >> � � � call dfftw_destroy_plan ( plan ) > > >> � � � write(*,*) 'Output array after inverse FFT:' > >> � � � do j = 1,N,1 > >> � � � � write(*,*) ' � �out1(',j,') = ',out1(j) > >> � � � enddo > > >> � � � do j = 1,N,1 > >> � � � � y =-1*cos(pi*(j)/N) > >> � � � � in1(j) = -1*sin(2*pi*y/N) > >> � � � � write(*,*) ' � �in1(',j,') = ',in1(j) > >> � � � enddo > >> � � � end > > >> �Regards, > >> �Sasidhar > > > Hello Sasidhar, > > > I don't know about your code, but derivatives of coefficients > > (constants?) are zero. Or are you trying to find derivatives of > > Chebyshev polynomials? These do have well defined forms in terms > > Chebyshev polynomials (although both 1st and 2nd kind are needed > > though). Or you can just stay with polynomials. > > > Maybe if you can give a clearer description of what you are really > > trying to do, then we may be better able to help. > > > Are you trying to do differentiation via FFTs? You know aliasing can > > be an issue here! > > Clay, What is a Chebychev transform? > > Jerry > -- > Engineering is the art of making what you want from things you can get. > �����������������������������������������������������������������������- Hide quoted text - > > - Show quoted text -Hey Jerry, Basically it is a form of expansion of a function in terms of Chebyshev polys where they are sampled on a nonuniformly spaced set of nodes. The nodes are found by cos(k*pi/N) where k goes from 0,1,2,3,.., N. Using a set from -1 to +1. Scaling and translation may be used for other domains. The transform is also self inverting. On the "chebyshev nodes", the discrete chebyshev functions (formed by sampling the continuous ones) are orthogonal with respect to the chebyshev weighting function. Chebys have a bunch of neat properties such as a composition of chebys is also a chebys! But by doing a cosine transformation on the domain, the cheby transform can be found via a FFT. And since it is self inverting, the OP wants to take a sample function and go around the loop and see if all of his steps are working properly. He needs to do a discrete cheby transform by the direct method so as to compare that halfway step to what he gets from the FFT method. Then he can try the inverse to see if he gets back to where he started as a sanity check. IHTH, Clay Clay
Reply by ●November 21, 20082008-11-21
>Basically it is a form of expansion of a function in terms of >Chebyshev polys where they are sampled on a nonuniformly spaced set of >nodes. The nodes are found by cos(k*pi/N) where k goes from >0,1,2,3,.., N. Using a set from -1 to +1. Scaling and translation may >be used for other domains. The transform is also self inverting. >On the "chebyshev nodes", the discrete chebyshev functions (formed by >sampling the continuous ones) are orthogonal with respect to the >chebyshev weighting function. > >Chebys have a bunch of neat properties such as a composition of chebys >is also a chebys! But by doing a cosine transformation on the domain, >the cheby transform can be found via a FFT. And since it is self >inverting, the OP wants to take a sample function and go around the >loop and see if all of his steps are working properly. > >He needs to do a discrete cheby transform by the direct method so as >to compare that halfway step to what he gets from the FFT method. Then >he can try the inverse to see if he gets back to where he started as a >sanity check.Hello Clay, Once the chebyshev coefficients are calculated, they can be manipulated in a way so that the inverse of the transformations gives a derivative of the original function to which chebyshev transformation has been carried upon. This is performed using the recurrence relation that I have shown above. "dout(N) = 0 dout(N-1) = 2*(N)*out(N) do j = N-2,1,-1 dout(j) = dout(j+2) + 2*(j+1)*out(j+1) enddo" Here dout array gives the manipulated coefficients and out array gives the chebyshev coefficients. However when I try to calculate the inverse, I am not getting the derivative function of the original. Regards, Sasidhar
Reply by ●November 22, 20082008-11-22
On Nov 21, 5:00�pm, "sasidhar" <sasidha...@gmail.com> wrote:> Here dout array gives the manipulated coefficients and out array gives the > chebyshev coefficients. However when I try to calculate the inverse, I am > not getting the derivative function of the original.For one thing, you may be missing a scale factor. FFTW computes the unnormalized transform, so you need to multiply somewhere by a scale factor yourself to get the inverse. See: http://www.fftw.org/doc/1d-Real_002deven-DFTs-_0028DCTs_0029.html More generally, you need to pay close attention to FFTW's normalization convention compared to the reference you are using for the Chebyshev polynomials; compare the formulas in the manual above to those in your references. Regards, Steven G. Johnson
Reply by ●November 22, 20082008-11-22
>For one thing, you may be missing a scale factor. FFTW computes the >unnormalized transform, so you need to multiply somewhere by a scale >factor yourself to get the inverse. See: > >http://www.fftw.org/doc/1d-Real_002deven-DFTs-_0028DCTs_0029.html > >More generally, you need to pay close attention to FFTW's >normalization convention compared to the reference you are using for >the Chebyshev polynomials; compare the formulas in the manual above to >those in your references. > >Regards, >Steven G. Johnson >Hello Steven, Thanks for the reply. I guess you are talking about normalization carried over with foactor of 2(n-1). I have taken care of that in the code right after doing the DCT transforms. call dfftw_plan_r2r_1d ( plan, N, in, out, ^ FFTW_REDFT00, FFTW_ESTIMATE ) call dfftw_execute ( plan ) call dfftw_destroy_plan ( plan ) out = out/6 I have already tested the inverse chebyshev transformation for the input function and it works fine. But if I perform the recurrence relation to calculate the differentiation coefficients and then use inverse chebyshev, I am not able to obtain the derivative function. Regards, Sasidhar
Reply by ●November 24, 20082008-11-24
On Nov 22, 6:31�pm, "sasidhar" <sasidha...@gmail.com> wrote:> I have already tested the inverse chebyshev transformation for the input > function and it works fine. But if I perform the recurrence relation to > calculate the differentiation coefficients and then use inverse chebyshev, > I am not able to obtain the derivative function.I haven't looked at your recurrence relation, but you should realize that there are several different normalization conventions for DCTs besides just overall scale factors. e.g. sometimes the x_0 element is scaled differently to make the transform unitary, whereas FFTW's normalization is chosen to enforce the correspondence with a DFT of real-even data. You should carefully compare FFTW's definition in the manual with the definition in the reference you are using for the recurrence relation, and adjust accordingly. Regards, Steven G. Johnson
Reply by ●May 16, 20092009-05-16
> >>For one thing, you may be missing a scale factor. FFTW computes the >>unnormalized transform, so you need to multiply somewhere by a scale >>factor yourself to get the inverse. See: >> >>http://www.fftw.org/doc/1d-Real_002deven-DFTs-_0028DCTs_0029.html >> >>More generally, you need to pay close attention to FFTW's >>normalization convention compared to the reference you are using for >>the Chebyshev polynomials; compare the formulas in the manual above to >>those in your references. >> >>Regards, >>Steven G. Johnson >> > >Hello Steven, > >Thanks for the reply. I guess you are talking about normalizationcarried>over with foactor of 2(n-1). I have taken care of that in the code right >after doing the DCT transforms. > > call dfftw_plan_r2r_1d ( plan, N, in, out, > ^ FFTW_REDFT00, FFTW_ESTIMATE ) > call dfftw_execute ( plan ) > call dfftw_destroy_plan ( plan ) > out = out/6 > >I have already tested the inverse chebyshev transformation for the input >function and it works fine. But if I perform the recurrence relation to >calculate the differentiation coefficients and then use inversechebyshev,>I am not able to obtain the derivative function. > >Regards, >Sasidhar > >Sasidhar, Did you ever get your spectral derivative code to work? I'm trying to do the same thing, using the DCT from FFTW and following the recurrence in Boyd's Spectral Methods book (pp. 298), but no luck for me either. Here's my fortran (pretty much the same, but using different DCT routines): call dfftw_plan_r2r_1d(plan1, n, in, out, FFTW_REDFT10, FFTW_ESTIMATE) call dfftw_execute_r2r(plan1, in, out) call dfftw_destroy_plan(plan1) d(n) = 0.0 d(n-1) = 0.0 do i=n-2,2,-1 d(i) = d(i+2) + 2*(i+1)*out(i+1) end do d(1) = (d(3) + 2*2*out(2))/2.0 d = d/(2.0*dble(n)) call dfftw_plan_r2r_1d(plan2, n, d, out, FFTW_REDFT01, FFTW_ESTIMATE) call dfftw_execute_r2r(plan2, d, out) call dfftw_destroy_plan(plan2) The test function I've been using is f(x) = x*exp(x) on [-1,1], and the results from the above give me something that looks sort of like the derivative flipped upside-down. Anyone out there successfully accomplished calculating derivatives with the discrete cosine transform that could lend a clue or two? Cheers, Josh
Reply by ●May 30, 20092009-05-30
Finally fixed it!
Here's a working Fortran subroutine for using the FFTW DCT to calculate
derivatives. The error in the above two posted codes is an off-by-one on
the multiplier in the recurrence. Also there is a lurking sign error that
I'm not quite able to trace from the FFTW definitions all the way through
to the recurrence.
subroutine dct_diff(dydx, y, n)
! use the discrete cosine transform from fftw to calculate the
derivative
integer, intent(in) :: n ! number of data points
double precision, intent(in), dimension(n) :: y ! function values
double precision, intent(out), dimension(n) :: dydx ! derivative
values
double precision, dimension(n) :: beta ! derivative coefficients
double precision, dimension(n) :: alpha ! function coefficients
integer :: plan1, plan2, i
integer :: n_logical ! the logical size of the transform, size of
! the corresponding real symmetric DFT
! the logical size depends on the type of transform, check the docs:
! http://www.fftw.org/doc/1d-Real_002deven-DFTs-_0028DCTs_0029.html
n_logical = 2*(n)
! forward DCT:
call dfftw_plan_r2r_1d(plan1, n, y, alpha, FFTW_REDFT10,
FFTW_ESTIMATE)
call dfftw_execute_r2r(plan1, y, alpha)
call dfftw_destroy_plan(plan1)
! recurrence for the derivative coefficients:
beta(n) = 0D0
beta(n-1) = 2D0 * dble(n-1) * alpha(n)
do i = n-1, 2, -1
beta(i-1) = beta(i+1) + 2D0 * dble(i-1) * alpha(i)
end do
beta = - beta ! this makes it work, but not sure why!
! inverse DCT:
call dfftw_plan_r2r_1d(plan2, n, beta, dydx, FFTW_REDFT01,
FFTW_ESTIMATE)
call dfftw_execute_r2r(plan2, beta, dydx)
call dfftw_destroy_plan(plan2)
! FFTW computes the un-normalized transforms, normalize by logical
size
dydx = dydx / dble(n_logical)
end subroutine dct_diff
