blob: 8f38d4b674790006ff216a09321eb681f0cbc268 [file] [log] [blame]
 // Copyright ©2018 The Gonum Authors. All rights reserved. // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. // This is a translation of the FFTPACK sint functions by // Paul N Swarztrauber, placed in the public domain at // http://www.netlib.org/fftpack/. package fftpack import "math" // Sinti initializes the array work which is used in subroutine Sint. // The prime factorization of n together with a tabulation of the // trigonometric functions are computed and stored in work. // // Input parameter // // n The length of the sequence to be transformed. The method // is most efficient when n+1 is a product of small primes. // // Output parameter // // work A work array with at least ceil(2.5*n) locations. // Different work arrays are required for different values // of n. The contents of work must not be changed between // calls of Sint. // // ifac An integer work array of length at least 15. func Sinti(n int, work []float64, ifac []int) { if len(work) < 5*(n+1)/2 { panic("fourier: short work") } if len(ifac) < 15 { panic("fourier: short ifac") } if n <= 1 { return } dt := math.Pi / float64(n+1) for k := 0; k < n/2; k++ { work[k] = 2 * math.Sin(float64(k+1)*dt) } Rffti(n+1, work[n/2:], ifac) } // Sint computes the Discrete Fourier Sine Transform of an odd // sequence x(i). The transform is defined below at output parameter x. // // Sint is the unnormalized inverse of itself since a call of Sint // followed by another call of Sint will multiply the input sequence // x by 2*(n+1). // // The array work which is used by subroutine Sint must be // initialized by calling subroutine Sinti(n,work). // // Input parameters // // n The length of the sequence to be transformed. The method // is most efficient when n+1 is the product of small primes. // // x An array which contains the sequence to be transformed. // // // work A work array with dimension at least ceil(2.5*n) // in the program that calls Sint. The work array must be // initialized by calling subroutine Sinti(n,work) and a // different work array must be used for each different // value of n. This initialization does not have to be // repeated so long as n remains unchanged thus subsequent // transforms can be obtained faster than the first. // // ifac An integer work array of length at least 15. // // Output parameters // // x for i=1,...,n // x(i)= the sum from k=1 to k=n // 2*x(k)*sin(k*i*pi/(n+1)) // // A call of Sint followed by another call of // Sint will multiply the sequence x by 2*(n+1). // Hence Sint is the unnormalized inverse // of itself. // // work Contains initialization calculations which must not be // destroyed between calls of Sint. // ifac Contains initialization calculations which must not be // destroyed between calls of Sint. func Sint(n int, x, work []float64, ifac []int) { if len(x) < n { panic("fourier: short sequence") } if len(work) < 5*(n+1)/2 { panic("fourier: short work") } if len(ifac) < 15 { panic("fourier: short ifac") } if n == 0 { return } sint1(n, x, work, work[n/2:], work[n/2+n+1:], ifac) } func sint1(n int, war, was, xh, x []float64, ifac []int) { const sqrt3 = 1.73205080756888 for i := 0; i < n; i++ { xh[i] = war[i] war[i] = x[i] } switch n { case 1: xh[0] *= 2 case 2: xh[0], xh[1] = sqrt3*(xh[0]+xh[1]), sqrt3*(xh[0]-xh[1]) default: x[0] = 0 for k := 0; k < n/2; k++ { kc := n - k - 1 t1 := xh[k] - xh[kc] t2 := was[k] * (xh[k] + xh[kc]) x[k+1] = t1 + t2 x[kc+1] = t2 - t1 } if n%2 != 0 { x[n/2+1] = 4 * xh[n/2] } rfftf1(n+1, x, xh, war, ifac) xh[0] = 0.5 * x[0] for i := 2; i < n; i += 2 { xh[i-1] = -x[i] xh[i] = xh[i-2] + x[i-1] } if n%2 == 0 { xh[n-1] = -x[n] } } for i := 0; i < n; i++ { x[i] = war[i] war[i] = xh[i] } }