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