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