dctw.c File Reference

Discrete Cosine Transform based on FFTW. More...

#include "libftsh.h"

Defines
#define	ENTEREXIT   0
Functions
void	dctw (REAL *fcn, Ftrigtw_Plan *main_plan, int dir)
void	rfftw_to_NR_format (REAL *data_nr, REAL *data_rfftw, int numpts, int dir)
void	dctw_init (Ftrigtw_Plan *main_plan, int flags, int numpts)

---

Detailed Description

Discrete Cosine Transform based on FFTW.

Summary:

• dctw_init initializes; do once.
• dctw does the transform.
• rfftw_to_NR_format is auxiliary routine for dctw

---

Function Documentation

void dctw	(	REAL *	fcn,
		Ftrigtw_Plan *	main_plan,
		int	dir
	)

Computes the disrete cosine transform of fcn and puts the coefficients back into fcn

This is the DCT at half integer points.

INPUTS:

• fcn -- array containing the data
• main_plan -- a structure holding precomputed trig functions, the plan used by fftw,and the number of points initialized by dctw_init
• dir -- the direction to transform. 1 is expansion and -1 evaluation.

OUTPUT: fcn.

• In forward mode (dir==1), puts the coefficients in fcn. \[f(k) = norm \sum_{j=0}^{numpts-1} f(j) \cos\left(\pi\frac{k(j+1/2)}{numpts}\right) \] set for $L^2$ normalization.

       numpts-1              k  (j+0.5)
  fcn[k]= SUM   fcn[j] cos( PI ----------------  ) *norm
         j=0                         numpts   
         set for L^2 normalization.
  

• In reverse mode inverts this process.

NOTES:

• comment better
• requires even numpts. checked by dctw_init
• precompute sqrt(2.0)?
• slow and awkward

Here is the call graph for this function:

Here is the caller graph for this function:

void dctw_init	(	Ftrigtw_Plan *	main_plan,
		int	flags,
		int	numpts
	)

This routine initializes for dctw

We initialize trigonometic things used by the lct and dct portions,
and call the initailizer for the fftw real transform (forward).

INPUTS: 
- main_plan -- a structure that will hold the things we generate
  Memory for its fields is allocated in this routine.
- flags -- flags passed directly to the rfftw initializer
- numpts -- the number of points on which we will transform
  MUST be even

OUTPUTS: main_plan is filled in:

NOTES:
- comment better

keep these always as double so we don't lose precision

Here is the caller graph for this function:

void rfftw_to_NR_format	(	REAL *	data_nr,
		REAL *	data_rfftw,
		int	numpts,
		int	dir
	)

Translates between the formats for the real discrete fourier transform used by the rfftw package and the Numerical Recipies routines. This helps allow the NR DCT routine to call the rfftw.

INPUTS:

• data_nr -- The data in the nr format: real[0],real[n/2],real[1],imag[1]... ...real[n/2-1],imag[n/2-1]
• data_rfftw--The data in the rfftw format: real[0],real[1],...real[n/2], imag[(n+1)/2-1],....,imag[1]
• numpts -- the length n of the vectors; MUST be even
• dir
  • 1 for rfftw -> nr
  • -1 for nr -> rfftw

OUTPUTS:

• if dir==1, data_nr is written onto
• if dir==-1, data_rfftw is written onto

NOTES:

• we must have REAL==fftw_real
• there is a sign convention difference on the imaginary parts realfft_pre(,,,,1); versus rfftw_one(rw_plan_forward,,);

Here is the caller graph for this function: