Mercurial > hg > sv-dependency-builds

annotate src/fftw-3.3.8/genfft/trig.ml @ 84:08ae793730bd

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