1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
|
// 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 cost functions by
// Paul N Swarztrauber, placed in the public domain at
// http://www.netlib.org/fftpack/.
package fftpack
import "math"
// Costi initializes the array work which is used in subroutine
// Cost. 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.
// Different work arrays are required for different values
// of n. The contents of work must not be changed between
// calls of Cost.
//
// ifac An integer work array of length at least 15.
func Costi(n int, work []float64, ifac []int) {
if len(work) < 3*n {
panic("fourier: short work")
}
if len(ifac) < 15 {
panic("fourier: short ifac")
}
if n < 4 {
return
}
dt := math.Pi / float64(n-1)
for k := 1; k < n/2; k++ {
fk := float64(k)
work[k] = 2 * math.Sin(fk*dt)
work[n-k-1] = 2 * math.Cos(fk*dt)
}
Rffti(n-1, work[n:], ifac)
}
// Cost computes the Discrete Fourier Cosine Transform of an even
// sequence x(i). The transform is defined below at output parameter x.
//
// Cost is the unnormalized inverse of itself since a call of Cost
// followed by another call of Cost will multiply the input sequence
// x by 2*(n-1). The transform is defined below at output parameter x
//
// The array work which is used by subroutine Cost must be
// initialized by calling subroutine Costi(n,work).
//
// Input parameters:
//
// n The length of the sequence x. n must be greater than 1.
// The method is most efficient when n-1 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 Cost. The work array must be
// initialized by calling subroutine Costi(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) = x(1)+(-1)**(i-1)*x(n)
// + the sum from k=2 to k=n-1
// 2*x(k)*cos((k-1)*(i-1)*pi/(n-1))
//
// A call of Cost followed by another call of
// Cost will multiply the sequence x by 2*(n-1).
// Hence Cost is the unnormalized inverse
// of itself.
//
// work Contains initialization calculations which must not be
// destroyed between calls of Cost.
//
// ifac An integer work array of length at least 15.
func Cost(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
}
switch n {
case 2:
x[0], x[1] = x[0]+x[1], x[0]-x[1]
case 3:
x1p3 := x[0] + x[2]
tx2 := 2 * x[1]
x[1] = x[0] - x[2]
x[0] = x1p3 + tx2
x[2] = x1p3 - tx2
default:
c1 := x[0] - x[n-1]
x[0] += x[n-1]
for k := 1; k < n/2; k++ {
kc := n - k - 1
t1 := x[k] + x[kc]
t2 := x[k] - x[kc]
c1 += work[kc] * t2
t2 *= work[k]
x[k] = t1 - t2
x[kc] = t1 + t2
}
if n%2 != 0 {
x[n/2] *= 2
}
Rfftf(n-1, x, work[n:], ifac)
xim2 := x[1]
x[1] = c1
for i := 3; i < n; i += 2 {
xi := x[i]
x[i] = x[i-2] - x[i-1]
x[i-1] = xim2
xim2 = xi
}
if n%2 != 0 {
x[n-1] = xim2
}
}
}
|
