Mercurial > hg > sv-dependency-builds

annotate src/fftw-3.3.5/genfft/trig.ml @ 83:ae30d91d2ffe

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Replace these with versions built using an older toolset (so as to avoid ABI compatibilities when linking on Ubuntu 14.04 for packaging purposes)
author Chris Cannam
date Fri, 07 Feb 2020 11:51:13 +0000
parents 2cd0e3b3e1fd
children
rev   line source
Chris@42 1 (*
Chris@42 2 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
Chris@42 3 * Copyright (c) 2003, 2007-14 Matteo Frigo
Chris@42 4 * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology
Chris@42 5 *
Chris@42 6 * This program is free software; you can redistribute it and/or modify
Chris@42 7 * it under the terms of the GNU General Public License as published by
Chris@42 8 * the Free Software Foundation; either version 2 of the License, or
Chris@42 9 * (at your option) any later version.
Chris@42 10 *
Chris@42 11 * This program is distributed in the hope that it will be useful,
Chris@42 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
Chris@42 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
Chris@42 14 * GNU General Public License for more details.
Chris@42 15 *
Chris@42 16 * You should have received a copy of the GNU General Public License
Chris@42 17 * along with this program; if not, write to the Free Software
Chris@42 18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
Chris@42 19 *
Chris@42 20 *)
Chris@42 21
Chris@42 22 (* trigonometric transforms *)
Chris@42 23 open Util
Chris@42 24
Chris@42 25 (* DFT of real input *)
Chris@42 26 let rdft sign n input =
Chris@42 27 Fft.dft sign n (Complex.real @@ input)
Chris@42 28
Chris@42 29 (* DFT of hermitian input *)
Chris@42 30 let hdft sign n input =
Chris@42 31 Fft.dft sign n (Complex.hermitian n input)
Chris@42 32
Chris@42 33 (* DFT real transform of vectors of two real numbers,
Chris@42 34 multiplication by (NaN I), and summation *)
Chris@42 35 let dft_via_rdft sign n input =
Chris@42 36 let f = rdft sign n input
Chris@42 37 in fun i ->
Chris@42 38 Complex.plus
Chris@42 39 [Complex.real (f i);
Chris@42 40 Complex.times (Complex.nan Expr.I) (Complex.imag (f i))]
Chris@42 41
Chris@42 42 (* Discrete Hartley Transform *)
Chris@42 43 let dht sign n input =
Chris@42 44 let f = Fft.dft sign n (Complex.real @@ input) in
Chris@42 45 (fun i ->
Chris@42 46 Complex.plus [Complex.real (f i); Complex.imag (f i)])
Chris@42 47
Chris@42 48 let trigI n input =
Chris@42 49 let twon = 2 * n in
Chris@42 50 let input' = Complex.hermitian twon input
Chris@42 51 in
Chris@42 52 Fft.dft 1 twon input'
Chris@42 53
Chris@42 54 let interleave_zero input = fun i ->
Chris@42 55 if (i mod 2) == 0
Chris@42 56 then Complex.zero
Chris@42 57 else
Chris@42 58 input ((i - 1) / 2)
Chris@42 59
Chris@42 60 let trigII n input =
Chris@42 61 let fourn = 4 * n in
Chris@42 62 let input' = Complex.hermitian fourn (interleave_zero input)
Chris@42 63 in
Chris@42 64 Fft.dft 1 fourn input'
Chris@42 65
Chris@42 66 let trigIII n input =
Chris@42 67 let fourn = 4 * n in
Chris@42 68 let twon = 2 * n in
Chris@42 69 let input' = Complex.hermitian fourn
Chris@42 70 (fun i ->
Chris@42 71 if (i == 0) then
Chris@42 72 Complex.real (input 0)
Chris@42 73 else if (i == twon) then
Chris@42 74 Complex.uminus (Complex.real (input 0))
Chris@42 75 else
Chris@42 76 Complex.antihermitian twon input i)
Chris@42 77 in
Chris@42 78 let dft = Fft.dft 1 fourn input'
Chris@42 79 in fun k -> dft (2 * k + 1)
Chris@42 80
Chris@42 81 let zero_extend n input = fun i ->
Chris@42 82 if (i >= 0 && i < n)
Chris@42 83 then input i
Chris@42 84 else Complex.zero
Chris@42 85
Chris@42 86 let trigIV n input =
Chris@42 87 let fourn = 4 * n
Chris@42 88 and eightn = 8 * n in
Chris@42 89 let input' = Complex.hermitian eightn
Chris@42 90 (zero_extend fourn (Complex.antihermitian fourn
Chris@42 91 (interleave_zero input)))
Chris@42 92 in
Chris@42 93 let dft = Fft.dft 1 eightn input'
Chris@42 94 in fun k -> dft (2 * k + 1)
Chris@42 95
Chris@42 96 let make_dct scale nshift trig =
Chris@42 97 fun n input ->
Chris@42 98 trig (n - nshift) (Complex.real @@ (Complex.times scale) @@
Chris@42 99 (zero_extend n input))
Chris@42 100 (*
Chris@42 101 * DCT-I: y[k] = sum x[j] cos(pi * j * k / n)
Chris@42 102 *)
Chris@42 103 let dctI = make_dct Complex.one 1 trigI
Chris@42 104
Chris@42 105 (*
Chris@42 106 * DCT-II: y[k] = sum x[j] cos(pi * (j + 1/2) * k / n)
Chris@42 107 *)
Chris@42 108 let dctII = make_dct Complex.one 0 trigII
Chris@42 109
Chris@42 110 (*
Chris@42 111 * DCT-III: y[k] = sum x[j] cos(pi * j * (k + 1/2) / n)
Chris@42 112 *)
Chris@42 113 let dctIII = make_dct Complex.half 0 trigIII
Chris@42 114
Chris@42 115 (*
Chris@42 116 * DCT-IV y[k] = sum x[j] cos(pi * (j + 1/2) * (k + 1/2) / n)
Chris@42 117 *)
Chris@42 118 let dctIV = make_dct Complex.half 0 trigIV
Chris@42 119
Chris@42 120 let shift s input = fun i -> input (i - s)
Chris@42 121
Chris@42 122 (* DST-x input := TRIG-x (input / i) *)
Chris@42 123 let make_dst scale nshift kshift jshift trig =
Chris@42 124 fun n input ->
Chris@42 125 Complex.real @@
Chris@42 126 (shift (- jshift)
Chris@42 127 (trig (n + nshift) (Complex.uminus @@
Chris@42 128 (Complex.times Complex.i) @@
Chris@42 129 (Complex.times scale) @@
Chris@42 130 Complex.real @@
Chris@42 131 (shift kshift (zero_extend n input)))))
Chris@42 132
Chris@42 133 (*
Chris@42 134 * DST-I: y[k] = sum x[j] sin(pi * j * k / n)
Chris@42 135 *)
Chris@42 136 let dstI = make_dst Complex.one 1 1 1 trigI
Chris@42 137
Chris@42 138 (*
Chris@42 139 * DST-II: y[k] = sum x[j] sin(pi * (j + 1/2) * k / n)
Chris@42 140 *)
Chris@42 141 let dstII = make_dst Complex.one 0 0 1 trigII
Chris@42 142
Chris@42 143 (*
Chris@42 144 * DST-III: y[k] = sum x[j] sin(pi * j * (k + 1/2) / n)
Chris@42 145 *)
Chris@42 146 let dstIII = make_dst Complex.half 0 1 0 trigIII
Chris@42 147
Chris@42 148 (*
Chris@42 149 * DST-IV y[k] = sum x[j] sin(pi * (j + 1/2) * (k + 1/2) / n)
Chris@42 150 *)
Chris@42 151 let dstIV = make_dst Complex.half 0 0 0 trigIV
Chris@42 152