github.com/gopherd/gonum@v0.0.4/dsp/fourier/internal/fftpack/sint.go

// 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]
	}
}