Mercurial > hg > sv-dependency-builds

annotate src/fftw-3.3.5/genfft/conv.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 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 open Complex
Chris@42 23 open Util
Chris@42 24
Chris@42 25 let polyphase m a ph i = a (m * i + ph)
Chris@42 26
Chris@42 27 let rec divmod n i =
Chris@42 28 if (i < 0) then
Chris@42 29 let (a, b) = divmod n (i + n)
Chris@42 30 in (a - 1, b)
Chris@42 31 else (i / n, i mod n)
Chris@42 32
Chris@42 33 let unpolyphase m a i = let (x, y) = divmod m i in a y x
Chris@42 34
Chris@42 35 let lift2 f a b i = f (a i) (b i)
Chris@42 36
Chris@42 37 (* convolution of signals A and B *)
Chris@42 38 let rec conv na a nb b =
Chris@42 39 let rec naive na a nb b i =
Chris@42 40 sigma 0 na (fun j -> (a j) @* (b (i - j)))
Chris@42 41
Chris@42 42 and recur na a nb b =
Chris@42 43 if (na <= 1 || nb <= 1) then
Chris@42 44 naive na a nb b
Chris@42 45 else
Chris@42 46 let p = polyphase 2 in
Chris@42 47 let ee = conv (na - na / 2) (p a 0) (nb - nb / 2) (p b 0)
Chris@42 48 and eo = conv (na - na / 2) (p a 0) (nb / 2) (p b 1)
Chris@42 49 and oe = conv (na / 2) (p a 1) (nb - nb / 2) (p b 0)
Chris@42 50 and oo = conv (na / 2) (p a 1) (nb / 2) (p b 1) in
Chris@42 51 unpolyphase 2 (function
Chris@42 52 0 -> fun i -> (ee i) @+ (oo (i - 1))
Chris@42 53 | 1 -> fun i -> (eo i) @+ (oe i)
Chris@42 54 | _ -> failwith "recur")
Chris@42 55
Chris@42 56
Chris@42 57 (* Karatsuba variant 1: (a+bx)(c+dx) = (ac+bdxx)+((a+b)(c+d)-ac-bd)x *)
Chris@42 58 and karatsuba1 na a nb b =
Chris@42 59 let p = polyphase 2 in
Chris@42 60 let ae = p a 0 and nae = na - na / 2
Chris@42 61 and ao = p a 1 and nao = na / 2
Chris@42 62 and be = p b 0 and nbe = nb - nb / 2
Chris@42 63 and bo = p b 1 and nbo = nb / 2 in
Chris@42 64 let ae = infinite nae ae and ao = infinite nao ao
Chris@42 65 and be = infinite nbe be and bo = infinite nbo bo in
Chris@42 66 let aeo = lift2 (@+) ae ao and naeo = nae
Chris@42 67 and beo = lift2 (@+) be bo and nbeo = nbe in
Chris@42 68 let ee = conv nae ae nbe be
Chris@42 69 and oo = conv nao ao nbo bo
Chris@42 70 and eoeo = conv naeo aeo nbeo beo in
Chris@42 71
Chris@42 72 let q = function
Chris@42 73 0 -> fun i -> (ee i) @+ (oo (i - 1))
Chris@42 74 | 1 -> fun i -> (eoeo i) @- ((ee i) @+ (oo i))
Chris@42 75 | _ -> failwith "karatsuba1" in
Chris@42 76 unpolyphase 2 q
Chris@42 77
Chris@42 78 (* Karatsuba variant 2:
Chris@42 79 (a+bx)(c+dx) = ((a+b)c-b(c-dxx))+x((a+b)c-a(c-d)) *)
Chris@42 80 and karatsuba2 na a nb b =
Chris@42 81 let p = polyphase 2 in
Chris@42 82 let ae = p a 0 and nae = na - na / 2
Chris@42 83 and ao = p a 1 and nao = na / 2
Chris@42 84 and be = p b 0 and nbe = nb - nb / 2
Chris@42 85 and bo = p b 1 and nbo = nb / 2 in
Chris@42 86 let ae = infinite nae ae and ao = infinite nao ao
Chris@42 87 and be = infinite nbe be and bo = infinite nbo bo in
Chris@42 88
Chris@42 89 let c1 = conv nae (lift2 (@+) ae ao) nbe be
Chris@42 90 and c2 = conv nao ao (nbo + 1) (fun i -> be i @- bo (i - 1))
Chris@42 91 and c3 = conv nae ae nbe (lift2 (@-) be bo) in
Chris@42 92
Chris@42 93 let q = function
Chris@42 94 0 -> lift2 (@-) c1 c2
Chris@42 95 | 1 -> lift2 (@-) c1 c3
Chris@42 96 | _ -> failwith "karatsuba2" in
Chris@42 97 unpolyphase 2 q
Chris@42 98
Chris@42 99 and karatsuba na a nb b =
Chris@42 100 let m = na + nb - 1 in
Chris@42 101 if (m < !Magic.karatsuba_min) then
Chris@42 102 recur na a nb b
Chris@42 103 else
Chris@42 104 match !Magic.karatsuba_variant with
Chris@42 105 1 -> karatsuba1 na a nb b
Chris@42 106 | 2 -> karatsuba2 na a nb b
Chris@42 107 | _ -> failwith "unknown karatsuba variant"
Chris@42 108
Chris@42 109 and via_circular na a nb b =
Chris@42 110 let m = na + nb - 1 in
Chris@42 111 if (m < !Magic.circular_min) then
Chris@42 112 karatsuba na a nb b
Chris@42 113 else
Chris@42 114 let rec find_min n = if n >= m then n else find_min (2 * n) in
Chris@42 115 circular (find_min 1) a b
Chris@42 116
Chris@42 117 in
Chris@42 118 let a = infinite na a and b = infinite nb b in
Chris@42 119 let res = array (na + nb - 1) (via_circular na a nb b) in
Chris@42 120 infinite (na + nb - 1) res
Chris@42 121
Chris@42 122 and circular n a b =
Chris@42 123 let via_dft n a b =
Chris@42 124 let fa = Fft.dft (-1) n a
Chris@42 125 and fb = Fft.dft (-1) n b
Chris@42 126 and scale = inverse_int n in
Chris@42 127 let fab i = ((fa i) @* (fb i)) @* scale in
Chris@42 128 Fft.dft 1 n fab
Chris@42 129
Chris@42 130 in via_dft n a b