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 | <http://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 | #include "fastfouriertransform.hpp" |
22 | #include "utilities.hpp" |
23 | #include <ql/math/fastfouriertransform.hpp> |
24 | #include <ql/math/array.hpp> |
25 | #include <complex> |
26 | #include <vector> |
27 | #include <functional> |
28 | |
29 | using namespace QuantLib; |
30 | using namespace boost::unit_test_framework; |
31 | |
32 | void FastFourierTransformTest::testSimple() { |
33 | BOOST_TEST_MESSAGE("Testing complex direct FFT..." ); |
34 | typedef std::complex<Real> cx; |
35 | cx a[] = { cx(0,0), cx(1,1), cx(3,3), cx(4,4), |
36 | cx(4,4), cx(3,3), cx(1,1), cx(0,0) }; |
37 | cx b[8]; |
38 | FastFourierTransform fft(3); |
39 | fft.transform(inBegin: a, inEnd: a+8, out: b); |
40 | cx expected[] = { cx(16,16), cx(-4.8284,-11.6569), |
41 | cx(0,0), cx(-0.3431,0.8284), |
42 | cx(0,0), cx(0.8284, -0.3431), |
43 | cx(0,0), cx(-11.6569,-4.8284) }; |
44 | for (size_t i = 0; i<8; i++) { |
45 | if ((std::fabs(x: b[i].real() - expected[i].real()) > 1.0e-2) || |
46 | (std::fabs(x: b[i].imag() - expected[i].imag()) > 1.0e-2)) |
47 | BOOST_ERROR("Convolution(" << i << ")\n" |
48 | << std::setprecision(4) << std::scientific |
49 | << " calculated: " << b[i] << "\n" |
50 | << " expected: " << expected[i]); |
51 | } |
52 | } |
53 | |
54 | void FastFourierTransformTest::testInverse() { |
55 | BOOST_TEST_MESSAGE("Testing convolution via inverse FFT..." ); |
56 | Array x(3); |
57 | x[0] = 1; |
58 | x[1] = 2; |
59 | x[2] = 3; |
60 | |
61 | size_t order = FastFourierTransform::min_order(inputSize: x.size())+1; |
62 | FastFourierTransform fft(order); |
63 | size_t nFrq = fft.output_size(); |
64 | std::vector< std::complex<Real> > ft (nFrq); |
65 | std::vector< Real > tmp (nFrq); |
66 | std::complex<Real> z = std::complex<Real>(); |
67 | |
68 | fft.inverse_transform(inBegin: x.begin(), inEnd: x.end(), out: ft.begin()); |
69 | for (Size i=0; i<nFrq; ++i) { |
70 | tmp[i] = std::norm(z: ft[i]); |
71 | ft[i] = z; |
72 | } |
73 | fft.inverse_transform(inBegin: tmp.begin(), inEnd: tmp.end(), out: ft.begin()); |
74 | |
75 | // 0 |
76 | Real calculated = ft[0].real() / nFrq; |
77 | Real expected = x[0]*x[0] + x[1]*x[1] + x[2]*x[2]; |
78 | if (fabs (x: calculated - expected) > 1.0e-10) |
79 | BOOST_ERROR("Convolution(0)\n" |
80 | << std::setprecision(16) << std::scientific |
81 | << " calculated: " << calculated << "\n" |
82 | << " expected: " << expected); |
83 | |
84 | // 1 |
85 | calculated = ft[1].real() / nFrq; |
86 | expected = x[0]*x[1] + x[1]*x[2]; |
87 | if (fabs (x: calculated - expected) > 1.0e-10) |
88 | BOOST_ERROR("Convolution(1)\n" |
89 | << std::setprecision(16) << std::scientific |
90 | << " calculated: " << calculated << "\n" |
91 | << " expected: " << expected); |
92 | |
93 | // 2 |
94 | calculated = ft[2].real() / nFrq; |
95 | expected = x[0]*x[2]; |
96 | if (fabs (x: calculated - expected) > 1.0e-10) |
97 | BOOST_ERROR("Convolution(1)\n" |
98 | << std::setprecision(16) << std::scientific |
99 | << " calculated: " << calculated << "\n" |
100 | << " expected: " << expected); |
101 | |
102 | } |
103 | |
104 | |
105 | |
106 | test_suite* FastFourierTransformTest::suite() { |
107 | auto* suite = BOOST_TEST_SUITE("fast fourier transform tests" ); |
108 | suite->add(QUANTLIB_TEST_CASE(&FastFourierTransformTest::testSimple)); |
109 | suite->add(QUANTLIB_TEST_CASE(&FastFourierTransformTest::testInverse)); |
110 | return suite; |
111 | } |
112 | |
113 | |