Mercurial > hg > sv-dependency-builds

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