Mercurial > hg > sv-dependency-builds

annotate src/fftw-3.3.5/genfft/trig.ml @ 56:af97cad61ff0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add updated build of PortAudio for OSX
author Chris Cannam <cannam@all-day-breakfast.com>
date Tue, 03 Jan 2017 15:10:52 +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