/src/quantlib/ql/math/fastfouriertransform.hpp

Source
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2006 Joseph Wang
 Copyright (C) 2009 Liquidnet Holdings, Inc.

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <https://www.quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

/*! \file fastfouriertransform.hpp
    \brief Fast Fourier Transform
*/

// Based on public domain code by Christopher Diggins

#ifndef quantlib_fast_fourier_transform_hpp
#define quantlib_fast_fourier_transform_hpp

#include <ql/errors.hpp>
#include <ql/types.hpp>
#include <vector>
#include <iterator>

namespace QuantLib {

    //! FFT implementation
    class FastFourierTransform {
      public:

        //! the minimum order required for the given input size
        static std::size_t min_order(std::size_t inputSize) {
            return static_cast<std::size_t>(
                std::ceil(std::log(static_cast<Real>(inputSize)) / M_LN2));
        }

        FastFourierTransform(std::size_t order)
        : cs_(order), sn_(order) {
            std::size_t m = static_cast<std::size_t>(1) << order;
            cs_[order - 1] = std::cos (2 * M_PI / m);
            sn_[order - 1] = std::sin (2 * M_PI / m);
            for (std::size_t i = order - 1; i > 0; --i) {
                cs_ [i - 1] = cs_[i]*cs_[i] - sn_[i]*sn_[i];
                sn_ [i - 1] = 2*sn_[i]*cs_[i];
            }
        }

        //! The required size for the output vector
        std::size_t output_size() const {
            return (static_cast<std::size_t>(1) << cs_.size());
        }

        //! FFT transform.
        /*! The output sequence must be allocated by the user */
        template<typename InputIterator, typename RandomAccessIterator>
        void transform(InputIterator inBegin, InputIterator inEnd,
                       RandomAccessIterator out) const {
            transform_impl(inBegin, inEnd, out, false);
        }

        //! Inverse FFT transform.
        /*! The output sequence must be allocated by the user. */
        template<typename InputIterator, typename RandomAccessIterator>
        void inverse_transform(InputIterator inBegin, InputIterator inEnd,
                               RandomAccessIterator out) const {
            transform_impl(inBegin, inEnd, out, true);
        }

      private:
        std::vector<double> cs_, sn_;

        template<typename InputIterator, typename RandomAccessIterator>
        void transform_impl(InputIterator inBegin, InputIterator inEnd,
                            RandomAccessIterator out,
                            bool inverse) const {
            typedef
                typename std::iterator_traits<RandomAccessIterator>::value_type
                                                                       complex;
            const std::size_t order = cs_.size();
            const auto N = std::size_t(static_cast<std::size_t>(1) << order);
            std::size_t i = 0;
            for (; inBegin != inEnd; ++i, ++inBegin) {
                *(out + bit_reverse(i, order)) = *inBegin;
            }
            QL_REQUIRE (i <= N, "FFT order is too small");
            for (std::size_t s = 1; s <= order; ++s) {
                std::size_t m = static_cast<std::size_t>(1) << s;
                complex w(1.0);
                complex wm(cs_[s-1], inverse ? sn_[s-1] : -sn_[s-1]);
                for (std::size_t j = 0; j < m/2; ++j) {
                    for (std::size_t k = j; k < N; k += m) {
                        complex t = w * (*(out + k + m/2));
                        complex u = *(out + k);
                        *(out + k) = u + t;
                        *(out + k + m/2) = u - t;
                    }
                    w *= wm;
                }
            }
        }

        static std::size_t bit_reverse(std::size_t x, std::size_t order) {
            std::size_t n = 0;
            for (std::size_t i = 0; i < order; ++i) {
                n <<= 1;
                n |= (x & 1);
                x >>= 1;
            }
            return n;
        }
    };

}

#endif

Line	Count	Source
1		/* -- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -- */
2
3		/*
4		Copyright (C) 2006 Joseph Wang
5		Copyright (C) 2009 Liquidnet Holdings, Inc.
6
7		This file is part of QuantLib, a free-software/open-source library
8		for financial quantitative analysts and developers - http://quantlib.org/
9
10		QuantLib is free software: you can redistribute it and/or modify it
11		under the terms of the QuantLib license. You should have received a
12		copy of the license along with this program; if not, please email
13		<quantlib-dev@lists.sf.net>. The license is also available online at
14		<https://www.quantlib.org/license.shtml>.
15
16		This program is distributed in the hope that it will be useful, but WITHOUT
17		ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18		FOR A PARTICULAR PURPOSE. See the license for more details.
19		*/
20
21		/*! \file fastfouriertransform.hpp
22		\brief Fast Fourier Transform
23		*/
24
25		// Based on public domain code by Christopher Diggins
26
27		#ifndef quantlib_fast_fourier_transform_hpp
28		#define quantlib_fast_fourier_transform_hpp
29
30		#include <ql/errors.hpp>
31		#include <ql/types.hpp>
32		#include <vector>
33		#include <iterator>
34
35		namespace QuantLib {
36
37		//! FFT implementation
38		class FastFourierTransform {
39		public:
40
41		//! the minimum order required for the given input size
42	0	static std::size_t min_order(std::size_t inputSize) {
43	0	return static_cast<std::size_t>(
44	0	std::ceil(std::log(static_cast<Real>(inputSize)) / M_LN2));
45	0	}
46
47		FastFourierTransform(std::size_t order)
48	0	: cs_(order), sn_(order) {
49	0	std::size_t m = static_cast<std::size_t>(1) << order;
50	0	cs_[order - 1] = std::cos (2 * M_PI / m);
51	0	sn_[order - 1] = std::sin (2 * M_PI / m);
52	0	for (std::size_t i = order - 1; i > 0; --i) {
53	0	cs_ [i - 1] = cs_[i]cs_[i] - sn_[i]sn_[i];
54	0	sn_ [i - 1] = 2sn_[i]cs_[i];
55	0	}
56	0	}
57
58		//! The required size for the output vector
59	0	std::size_t output_size() const {
60	0	return (static_cast<std::size_t>(1) << cs_.size());
61	0	}
62
63		//! FFT transform.
64		/! The output sequence must be allocated by the user /
65		template<typename InputIterator, typename RandomAccessIterator>
66		void transform(InputIterator inBegin, InputIterator inEnd,
67	0	RandomAccessIterator out) const {
68	0	transform_impl(inBegin, inEnd, out, false);
69	0	} Unexecuted instantiation: void QuantLib::FastFourierTransform::transform<std::__1::__wrap_iter<std::__1::complex<double>>, std::__1::__wrap_iter<std::__1::complex<double>> >(std::__1::__wrap_iter<std::__1::complex<double>>, std::__1::__wrap_iter<std::__1::complex<double>>, std::__1::__wrap_iter<std::__1::complex<double>>) const Unexecuted instantiation: void QuantLib::FastFourierTransform::transform<std::__1::__wrap_iter<double>, std::__1::__wrap_iter<std::__1::complex<double>> >(std::__1::__wrap_iter<double>, std::__1::__wrap_iter<double>, std::__1::__wrap_iter<std::__1::complex<double>>) const Unexecuted instantiation: void QuantLib::FastFourierTransform::transform<double, std::__1::__wrap_iter<std::__1::complex<double>> >(double, double, std::__1::__wrap_iter<std::__1::complex<double>*>) const
70
71		//! Inverse FFT transform.
72		/! The output sequence must be allocated by the user. /
73		template<typename InputIterator, typename RandomAccessIterator>
74		void inverse_transform(InputIterator inBegin, InputIterator inEnd,
75		RandomAccessIterator out) const {
76		transform_impl(inBegin, inEnd, out, true);
77		}
78
79		private:
80		std::vector<double> cs_, sn_;
81
82		template<typename InputIterator, typename RandomAccessIterator>
83		void transform_impl(InputIterator inBegin, InputIterator inEnd,
84		RandomAccessIterator out,
85	0	bool inverse) const {
86	0	typedef
87	0	typename std::iterator_traits<RandomAccessIterator>::value_type
88	0	complex;
89	0	const std::size_t order = cs_.size();
90	0	const auto N = std::size_t(static_cast<std::size_t>(1) << order);
91	0	std::size_t i = 0;
92	0	for (; inBegin != inEnd; ++i, ++inBegin) {
93	0	(out + bit_reverse(i, order)) = inBegin;
94	0	}
95	0	QL_REQUIRE (i <= N, "FFT order is too small");
96	0	for (std::size_t s = 1; s <= order; ++s) {
97	0	std::size_t m = static_cast<std::size_t>(1) << s;
98	0	complex w(1.0);
99	0	complex wm(cs_[s-1], inverse ? sn_[s-1] : -sn_[s-1]);
100	0	for (std::size_t j = 0; j < m/2; ++j) {
101	0	for (std::size_t k = j; k < N; k += m) {
102	0	complex t = w * (*(out + k + m/2));
103	0	complex u = *(out + k);
104	0	*(out + k) = u + t;
105	0	*(out + k + m/2) = u - t;
106	0	}
107	0	w *= wm;
108	0	}
109	0	}
110	0	} Unexecuted instantiation: void QuantLib::FastFourierTransform::transform_impl<std::__1::__wrap_iter<std::__1::complex<double>>, std::__1::__wrap_iter<std::__1::complex<double>> >(std::__1::__wrap_iter<std::__1::complex<double>>, std::__1::__wrap_iter<std::__1::complex<double>>, std::__1::__wrap_iter<std::__1::complex<double>>, bool) const Unexecuted instantiation: void QuantLib::FastFourierTransform::transform_impl<std::__1::__wrap_iter<double>, std::__1::__wrap_iter<std::__1::complex<double>> >(std::__1::__wrap_iter<double>, std::__1::__wrap_iter<double>, std::__1::__wrap_iter<std::__1::complex<double>>, bool) const Unexecuted instantiation: void QuantLib::FastFourierTransform::transform_impl<double, std::__1::__wrap_iter<std::__1::complex<double>> >(double, double, std::__1::__wrap_iter<std::__1::complex<double>*>, bool) const
111
112	0	static std::size_t bit_reverse(std::size_t x, std::size_t order) {
113	0	std::size_t n = 0;
114	0	for (std::size_t i = 0; i < order; ++i) {
115	0	n <<= 1;
116	0	n \|= (x & 1);
117	0	x >>= 1;
118	0	}
119	0	return n;
120	0	}
121		};
122
123		}
124
125		#endif

Coverage Report

Created: 2026-01-25 06:59