Sign in

username:

password:



Not a member?

Search compdsp



Search tips

comp.dsp by Keywords

Adaptive Filter | ADPCM | ADSP | ADSP-2181 | Aliasing | AMR | Anti-Aliasing | ARMA | Autocorrelation | AutoCovariance | Beamforming | Bessel | Blackfin | Butterworth | C6713 | CCS | Chebyshev | CIC Filter | Circular Convolution | Code Composer Studio | Comb Filter | Compression | Convolution | Cross Correlation | DCT | Decimation | Deconvolution | Demodulation | DM642 | DSP Boards | DSP/BIOS | DTMF | Echo Cancellation | Equalization | Equalizer | ETSI | EZLITE (Ez-kit Lite) | FFT | FFTW | FIR Filter | Fixed Point | FSK | G.711 | G.723 | G.729 | Gaussian Noise | Goertzel | GPIO | Hilbert Transform | IFFT | IIR Filter | Interpolation | Invariance | JTAG | Kalman | Laplace Transform | Levinson | LPC | McBSP | MIPS | Modulation | MPEG | Multirate | Notch Filter | Nyquist | OFDM | Oversampling | Pink Noise | Pitch | PLL | Polyphase | QAM | QDMA | Quantization | Quantizer | Radar | Random Noise | Reed Solomon | Remez | Resampling | RTDX | Sampling | Sharc | TI C6711 | Undersampling | Viterbi | Wavelets | White Noise | Wiener Filter | Windowing | XDS510PP | Z Transform

Sponsor

Industry's highest performing at the lowest power DSPs now as low as $5.00*
Start development today!
*volume pricing for 10ku

Discussion Groups

Free Online Books

See Also

Embedded SystemsFPGAElectronics

Discussion Groups | Comp.DSP | DCT derivative

There are 10 messages in this thread.

You are currently looking at messages 0 to 10.


DCT derivative - sasidhar - 2008-11-21 11:42:00

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
______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

Re: DCT derivative - 2008-11-21 14:10:00



On Nov 21, 11:42=A0am, "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.
>
> =A0 =A0 =A0 program test
> =A0 =A0 =A0 implicit none
> #include "fftw3.f"
>
> =A0 =A0 =A0 INTEGER, PARAMETER =A0:: =A0N=3D4
> =A0 =A0 =A0 INTEGER(8) :: =A0plan
> =A0 =A0 =A0 REAL(8), DIMENSION(N) :: =A0in, out, out1, in1
> =A0 =A0 =A0 REAL(8), DIMENSION(N) :: dout
> =A0 =A0 =A0 REAL :: =A0pi, y
> =A0 =A0 =A0 INTEGER :: j, p
>
> =A0 =A0 =A0 pi =3D 3.14
> =A0 =A0 =A0 write(*,*) 'Input array:'
> =A0 =A0 =A0 do j =3D 1,N,1
> =A0 =A0 =A0 =A0 y =3D-1*cos(pi*(j)/N)
> =A0 =A0 =A0 =A0 in(j) =3D cos(2*pi*y/N)
> =A0 =A0 =A0 enddo
>
> =A0 =A0 =A0 call dfftw_plan_r2r_1d ( plan, N, in, out,
> =A0 =A0 =A0^ =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 FFTW_REDFT00=
, FFTW_ESTIMATE )
> =A0 =A0 =A0 call dfftw_execute ( plan )
> =A0 =A0 =A0 call dfftw_destroy_plan ( plan )
> =A0 =A0 =A0 out =3D out/6
> =A0 =A0 =A0 do j =3D 1,N,1
> =A0 =A0 =A0 =A0 write(*,*) ' =A0 =A0out(',j,') =3D ',out(j)
> =A0 =A0 =A0 enddo
>
> =A0 =A0 =A0 dout(N) =3D 0
> =A0 =A0 =A0 dout(N-1) =3D 2*(N)*out(N)
>
> =A0 =A0 =A0 do j =3D N-2,1,-1
> =A0 =A0 =A0 =A0 dout(j) =3D dout(j+2) + 2*(j+1)*out(j+1)
> =A0 =A0 =A0 enddo
>
> =A0 =A0 =A0 call dfftw_plan_r2r_1d(plan, N+1, dout, out1,
> =A0 =A0 =A0^ =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 FFTW_REDFT00=
, FFTW_ESTIMATE)
> =A0 =A0 =A0 call dfftw_execute ( plan )
> =A0 =A0 =A0 call dfftw_destroy_plan ( plan )
>
> =A0 =A0 =A0 write(*,*) 'Output array after inverse FFT:'
> =A0 =A0 =A0 do j =3D 1,N,1
> =A0 =A0 =A0 =A0 write(*,*) ' =A0 =A0out1(',j,') =3D ',out1(j)
> =A0 =A0 =A0 enddo
>
> =A0 =A0 =A0 do j =3D 1,N,1
> =A0 =A0 =A0 =A0 y =3D-1*cos(pi*(j)/N)
> =A0 =A0 =A0 =A0 in1(j) =3D -1*sin(2*pi*y/N)
> =A0 =A0 =A0 =A0 write(*,*) ' =A0 =A0in1(',j,') =3D ',in1(j)
> =A0 =A0 =A0 enddo
> =A0 =A0 =A0 end
>
> =A0Regards,
> =A0Sasidhar

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
______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

Re: DCT derivative - Jerry Avins - 2008-11-21 14:52:00

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.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

Re: DCT derivative - 2008-11-21 15:48:00

On Nov 21, 2:52=A0pm, 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 F=
FTW
> >> for cosine(x). However I am not getting back the right values when I
> >> perform inverse chebyshev transform. Can anyone check the code below a=
nd
> >> tell me the possible errors. Thank you.
>
> >> =A0 =A0 =A0 program test
> >> =A0 =A0 =A0 implicit none
> >> #include "fftw3.f"
>
> >> =A0 =A0 =A0 INTEGER, PARAMETER =A0:: =A0N=3D4
> >> =A0 =A0 =A0 INTEGER(8) :: =A0plan
> >> =A0 =A0 =A0 REAL(8), DIMENSION(N) :: =A0in, out, out1, in1
> >> =A0 =A0 =A0 REAL(8), DIMENSION(N) :: dout
> >> =A0 =A0 =A0 REAL :: =A0pi, y
> >> =A0 =A0 =A0 INTEGER :: j, p
>
> >> =A0 =A0 =A0 pi =3D 3.14
> >> =A0 =A0 =A0 write(*,*) 'Input array:'
> >> =A0 =A0 =A0 do j =3D 1,N,1
> >> =A0 =A0 =A0 =A0 y =3D-1*cos(pi*(j)/N)
> >> =A0 =A0 =A0 =A0 in(j) =3D cos(2*pi*y/N)
> >> =A0 =A0 =A0 enddo
>
> >> =A0 =A0 =A0 call dfftw_plan_r2r_1d ( plan, N, in, out,
> >> =A0 =A0 =A0^ =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 FFTW_REDF=
T00, FFTW_ESTIMATE )
> >> =A0 =A0 =A0 call dfftw_execute ( plan )
> >> =A0 =A0 =A0 call dfftw_destroy_plan ( plan )
> >> =A0 =A0 =A0 out =3D out/6
> >> =A0 =A0 =A0 do j =3D 1,N,1
> >> =A0 =A0 =A0 =A0 write(*,*) ' =A0 =A0out(',j,') =3D ',out(j)
> >> =A0 =A0 =A0 enddo
>
> >> =A0 =A0 =A0 dout(N) =3D 0
> >> =A0 =A0 =A0 dout(N-1) =3D 2*(N)*out(N)
>
> >> =A0 =A0 =A0 do j =3D N-2,1,-1
> >> =A0 =A0 =A0 =A0 dout(j) =3D dout(j+2) + 2*(j+1)*out(j+1)
> >> =A0 =A0 =A0 enddo
>
> >> =A0 =A0 =A0 call dfftw_plan_r2r_1d(plan, N+1, dout, out1,
> >> =A0 =A0 =A0^ =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 FFTW_REDF=
T00, FFTW_ESTIMATE)
> >> =A0 =A0 =A0 call dfftw_execute ( plan )
> >> =A0 =A0 =A0 call dfftw_destroy_plan ( plan )
>
> >> =A0 =A0 =A0 write(*,*) 'Output array after inverse FFT:'
> >> =A0 =A0 =A0 do j =3D 1,N,1
> >> =A0 =A0 =A0 =A0 write(*,*) ' =A0 =A0out1(',j,') =3D ',out1(j)
> >> =A0 =A0 =A0 enddo
>
> >> =A0 =A0 =A0 do j =3D 1,N,1
> >> =A0 =A0 =A0 =A0 y =3D-1*cos(pi*(j)/N)
> >> =A0 =A0 =A0 =A0 in1(j) =3D -1*sin(2*pi*y/N)
> >> =A0 =A0 =A0 =A0 write(*,*) ' =A0 =A0in1(',j,') =3D ',in1(j)
> >> =A0 =A0 =A0 enddo
> >> =A0 =A0 =A0 end
>
> >> =A0Regards,
> >> =A0Sasidhar
>
> > 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.
> =AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=
=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=
=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF- Hide qu=
oted 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
______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

Re: DCT derivative - sasidhar - 2008-11-21 17:00:00

>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
______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

Re: DCT derivative - Steven G. Johnson - 2008-11-22 12:18:00

On Nov 21, 5:00=A0pm, "sasidhar" <sasidha...@gmail.com> wrote:
> Here dout array gives the manipulated coefficients and out array gives th=
e
> 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
______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

Re: DCT derivative - sasidhar - 2008-11-22 18:31:00

>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
______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

Re: DCT derivative - Steven G. Johnson - 2008-11-24 21:27:00

On Nov 22, 6:31=A0pm, "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
______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

Re: DCT derivative - jstults - 2009-05-16 14:49:00

>
>>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
>
>

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
______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

Re: DCT derivative - jstults - 2009-05-30 21:57:00

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
______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.