github.com/gopherd/gonum@v0.0.4/dsp/fourier/internal/fftpack/cosq.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 cosq functions by
// Paul N Swarztrauber, placed in the public domain at
// http://www.netlib.org/fftpack/.

package fftpack

import "math"

// Cosqi initializes the array work which is used in both Cosqf
// and Cosqb. 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 parameters
//
// work    A work array which must be dimensioned at least 3*n.
//         The same work array can be used for both Cosqf and Cosqb
//         as long as n remains unchanged. Different work arrays
//         are required for different values of n. The contents of
//         work must not be changed between calls of Cosqf or Cosqb.
//
// ifac    An integer work array of length at least 15.
func Cosqi(n int, work []float64, ifac []int) {
	if len(work) < 3*n {
		panic("fourier: short work")
	}
	if len(ifac) < 15 {
		panic("fourier: short ifac")
	}
	dt := 0.5 * math.Pi / float64(n)
	for k := range work[:n] {
		work[k] = math.Cos(float64(k+1) * dt)
	}
	Rffti(n, work[n:], ifac)
}

// Cosqf computes the Fast Fourier Transform of quarter wave data.
// That is, Cosqf computes the coefficients in a cosine series
// representation with only odd wave numbers. The transform is
// defined below at output parameter x.
//
// Cosqb is the unnormalized inverse of Cosqf since a call of Cosqf
// followed by a call of Cosqb will multiply the input sequence x
// by 4*n.
//
// The array work which is used by subroutine Cosqf must be
// initialized by calling subroutine Cosqi(n,work).
//
// Input parameters
//
// n       The length of the array x to be transformed. The method
//         is most efficient when n is a product of small primes.
//
// x       An array which contains the sequence to be transformed.
//
// work    A work array which must be dimensioned at least 3*n
//         in the program that calls Cosqf. The work array must be
//         initialized by calling subroutine Cosqi(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=0, ..., n-1
//           x[i] = x[i] + the sum from k=0 to k=n-2 of
//               2*x[k]*cos((2*i+1)*k*pi/(2*n))
//
//         A call of Cosqf followed by a call of
//         Cosqb will multiply the sequence x by 4*n.
//         Therefore Cosqb is the unnormalized inverse
//         of Cosqf.
//
// work    Contains initialization calculations which must not
//         be destroyed between calls of Cosqf or Cosqb.
func Cosqf(n int, x, work []float64, ifac []int) {
	if len(x) < n {
		panic("fourier: short sequence")
	}
	if len(work) < 3*n {
		panic("fourier: short work")
	}
	if len(ifac) < 15 {
		panic("fourier: short ifac")
	}
	if n < 2 {
		return
	}
	if n == 2 {
		tsqx := math.Sqrt2 * x[1]
		x[1] = x[0] - tsqx
		x[0] += tsqx
		return
	}
	cosqf1(n, x, work, work[n:], ifac)
}

func cosqf1(n int, x, w, xh []float64, ifac []int) {
	for k := 1; k < (n+1)/2; k++ {
		kc := n - k
		xh[k] = x[k] + x[kc]
		xh[kc] = x[k] - x[kc]
	}
	if n%2 == 0 {
		xh[(n+1)/2] = 2 * x[(n+1)/2]
	}
	for k := 1; k < (n+1)/2; k++ {
		kc := n - k
		x[k] = w[k-1]*xh[kc] + w[kc-1]*xh[k]
		x[kc] = w[k-1]*xh[k] - w[kc-1]*xh[kc]
	}
	if n%2 == 0 {
		x[(n+1)/2] = w[(n-1)/2] * xh[(n+1)/2]
	}
	Rfftf(n, x, xh, ifac)
	for i := 2; i < n; i += 2 {
		x[i-1], x[i] = x[i-1]-x[i], x[i-1]+x[i]
	}
}

// Cosqb computes the Fast Fourier Transform of quarter wave data.
// That is, Cosqb computes a sequence from its representation in
// terms of a cosine series with odd wave numbers. The transform
// is defined below at output parameter x.
//
// Cosqf is the unnormalized inverse of Cosqb since a call of Cosqb
// followed by a call of Cosqf will multiply the input sequence x
// by 4*n.
//
// The array work which is used by subroutine Cosqb must be
// initialized by calling subroutine Cosqi(n,work).
//
//
// Input parameters
//
// n       The length of the array x to be transformed. The method
//         is most efficient when n is a product of small primes.
//
// x       An array which contains the sequence to be transformed.
//
// work    A work array which must be dimensioned at least 3*n
//         in the program that calls Cosqb. The work array must be
//         initialized by calling subroutine Cosqi(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=0, ..., n-1
//           x[i]= the sum from k=0 to k=n-1 of
//             4*x[k]*cos((2*k+1)*i*pi/(2*n))
//
//         A call of Cosqb followed by a call of
//         Cosqf will multiply the sequence x by 4*n.
//         Therefore Cosqf is the unnormalized inverse
//         of Cosqb.
//
// work    Contains initialization calculations which must not
//         be destroyed between calls of Cosqb or Cosqf.
func Cosqb(n int, x, work []float64, ifac []int) {
	if len(x) < n {
		panic("fourier: short sequence")
	}
	if len(work) < 3*n {
		panic("fourier: short work")
	}
	if len(ifac) < 15 {
		panic("fourier: short ifac")
	}

	if n < 2 {
		x[0] *= 4
		return
	}
	if n == 2 {
		x[0], x[1] = 4*(x[0]+x[1]), 2*math.Sqrt2*(x[0]-x[1])
		return
	}
	cosqb1(n, x, work, work[n:], ifac)
}

func cosqb1(n int, x, w, xh []float64, ifac []int) {
	for i := 2; i < n; i += 2 {
		x[i-1], x[i] = x[i-1]+x[i], x[i]-x[i-1]
	}
	x[0] *= 2
	if n%2 == 0 {
		x[n-1] *= 2
	}
	Rfftb(n, x, xh, ifac)
	for k := 1; k < (n+1)/2; k++ {
		kc := n - k
		xh[k] = w[k-1]*x[kc] + w[kc-1]*x[k]
		xh[kc] = w[k-1]*x[k] - w[kc-1]*x[kc]
	}
	if n%2 == 0 {
		x[(n+1)/2] *= 2 * w[(n-1)/2]
	}
	for k := 1; k < (n+1)/2; k++ {
		x[k] = xh[k] + xh[n-k]
		x[n-k] = xh[k] - xh[n-k]
	}
	x[0] *= 2
}