This class serves as a wrapper for the FFTW Fourier transform
libraries for use with complex vectors. It serves two major purposes. First, it
encapsulates the allocation of the forward and reverse plan structures used by
the FFTW functions to keep track of information that only needs to be
calculated once for transforms of a given size. Second, it provides routines
for automatically converting from the std::complex<T>
classes used by VecMat to the fftw_complex structure used
by FFTW. It also provides routines for operating directly on vectors
of fftw_complex values.
Use of this class requires you to include the header file "vm_fftw.h", which automatically includes the "fftw.h" and "rfftw.h" headers.
explicit FFTWPlan(size_t len = 0)Construct a plan of size len. If len is zero then the
resize() member function can be used to set a valid plan size, or
the plan size will be set automatically when fourier() or
ifourier() are called. Note that any size can be used, but the
transforms will be very slow if len has large prime factors.
FFTWPlan(const FFTWPlan& plan)This constructs an exact copy of plan. Note that two plans of
the same size will function identically, so this constructor is only provided
for the sake of completeness. For example, it might be necessary if you want to
include an FFTWPlan object as a data member of another class.
FFTWPlan& operator=(const FFTWPlan& plan)The assignment operator simply resizes the plan it is called for to the size
of plan. Like the copy constructor, this is included only for
the sake of completeness.
size_t size() constThis member function returns the current size of the plan.
void resize(size_t len)len. Resizing to zero can be done
to free up any memory being used by the plan once it is no longer needed.
void fourier(const Vec_CPLX_REAL& array, Vec_CPLX_REAL&
arrayft)Performs the forward Fourier transform of array, and stores the
results in arrayft. Pointer casting is used to allow the function
to send pointers of type FFTW_CPLX *, that point to the vector's
data. This will only work if the CPLX_REAL class uses the same
storage scheme as FFTW_CPLX. This means that it must only have two
data members, and they must be stored with the real value first, and the
imaginary value second. The vectors array and arrayft
must be the same size. If the plan is not set to this size, it will be resized.
If the vectors are size zero, the function simply returns, and the plan size is
not altered.
void fourier(const Vec_FFTW_CPLX& array, Vec_FFTW_CPLX&
arrayft)Performs the forward Fourier transform of array, and stores the
results in arrayft. The data in array is directly
operated on. The vectors array and arrayft must be
the same size. If the plan is not set to this size, it will be resized. If the
vectors are size zero, the function simply returns, and the plan size is not
altered.
void ifourier(const Vec_CPLX_REAL& arrayft, Vec_CPLX_REAL&
array)Performs the inverse Fourier transform of arrayft, and stores
the results in array. Pointer casting is used to allow the function
to send pointers of type FFTW_CPLX *, that point to the vector's
data. This will only work if the CPLX_REAL class uses the same
storage scheme as FFTW_CPLX. This means that it must only have two
data members, and they must be stored with the real value first, and the
imaginary value second. The vectors array and arrayft
must be the same size. If the plan is not set to this size, it will be resized.
If the vectors are size zero, the function simply returns, and the plan size is
not altered.
void ifourier(const Vec_FFTW_CPLX& arrayft, Vec_FFTW_CPLX&
array)Performs the inverse Fourier transform of arrayft, and stores
the results in array. The data in array is directly
operated on. The vectors array and arrayft must be
the same size. If the plan is not set to this size, it will be resized. If
the vectors are size zero, the function simply returns, and the plan size is
not altered.
This class serves as a wrapper for the FFTW Fourier transform libraries for use with real vectors It encapsulates the allocation of the forward and reverse plan structures used by the FFTW functions to keep track of information that only needs to be calculated once for transforms of a given size. Note that the FFTW libraries have a special format for storing the Fourier transforms of real arrays. The real part of the positive frequencies are stored in the indices where the positive frequencies would normally be stored, and the imaginary part of the positive frequencies are stored in the indices where the negative frequencies would normally be stored. This allows the result to be stored in a real array of the same size as the input array. See the FFTW documentation for more details.
Use of this class requires you to include the header file "vm_fftw.h".
explicite RFFTWPlan(size_t len = 0)Construct a plan of size len. If len is zero then
the resize() member function can be used to set a valid plan size,
or the plan size will be set automatically when fourier() or
ifourier() are called. Note that any size can be used, but the
transforms will be very slow if len has large prime factors.
RFFTWPlan(const RFFTWPlan& plan)This constructs an exact copy of plan. Note that two plans of
the same size will function identically, so this constructor is only provided
for the sake of completeness. For example, it might be necessary if you want to
include an RFFTWPlan object as a data member of another class.
RFFTWPlan& operator=(const RFFTWPlan& plan)The assignment operator simply resizes the plan it is called for to the
size of plan. Like the copy constructor, this is included only for
the sake of completeness.
size_t size() constThis member function returns the current size of the plan.
void resize(size_t len)Reinitializes the plan to size len. Resizing to zero can be
done to free up any memory being used by the plan once it is no longer
needed.
void fourier(const Vec_REAL& array, Vec_REAL&
arrayft)Performs the forward Fourier transform of array, and stores the
results in arrayft. This function operates directly on the data of
the vectors, so no copying is done. The vectors array and
arrayft must be the same size. If the plan is not set to this
size, it will be resized. If the vectors are size zero, the function simply
returns, and the plan size is not altered.
void ifourier(const Vec_REAL& arrayft, Vec_REAL&
array)Performs the inverse Fourier transform of arrayft, and stores
the results in array. This function operates directly on the data
of the vectors, so no copying is done. The vectors array and
arrayft must be the same size. If the plan is not set to this
size, it will be resized. If the vectors are size zero, the function simply
returns, and the plan size is not altered. Note that due to how the
FFTW transforms for real data work, the inverse transform will corrupt
the contents of arrayft.
Vec_REAL power(const Vec_REAL& array, size_t res)This function calculates the power-spectrum of array. The
parameter res sets the window size to be used for the calculation.
The power is estimated by averaging over 50% overlapping periodograms of length
res. A Welch windowing function is used to eliminate aperiodicity
artifacts. The resulting power is returned in a real vector of length
res. Note that only the first res / 2 + 1 frequency
bins will be non-zero. This returned vector can therefore be immediately
inverse Fourier transformed to give the auto-correlation function. As with
fourier() and ifourier(), if the size of the plan is
not equal to res, then it will be resized. If res is
zero then an empty vector will be returned, and the plan will not be
resized.
Vec_REAL cross_spectrum(const Vec_REAL& array1, const
Vec_REAL& array2, size_t res)This function calculates the cross-spectrum of array1 and
array2. The parameter res sets the window size to be
used for the calculation. The cross-spectrum is estimated by averaging over 50%
overlapping windows of length res. A Welch windowing function is
used to eliminate aperiodicity artifacts. The result is returned in a real
vector of size res. This returned vector can then be immediately
inverse Fourier transformed to give the cross-correlation function. As with
fourier() and ifourier(), if the size of the plan is
not equal to res, then it will be resized. If res is
zero, then an empty vector will be returned, and the plan will not be resized.
Both arrays must be non-empty and the same size.
Vec_REAL convolve(const Vec_REAL& signal, const Vec_REAL&
response, size_t pad = 0)This function convolves signal with response. The
two vectors do not need to be the same length, but response must
not be longer than signal. If necessary, response
will be padded with zeros to give the proper length. This padding is done by
the convention that response[0] is the response function at time
zero, and that the last half of response represents negative time.
This means that the zero padding is done in the middle, not at the end. The
optional parameter pad can be used to pad signal to
remove aperiodicity artifacts. If signal is not periodic, then a
padding value equal of at least half the size of response should
be used. For best results, pick a value of pad such that the
padded length has only small prime factors. If the plan size is not equal to
the size of the padded signal, then it will be resized.