blob: 3c24f8d4c2abaa8488f867f7f0e48f4224298293 [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 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 }