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
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
|
// 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 sinq functions by
// Paul N Swarztrauber, placed in the public domain at
// http://www.netlib.org/fftpack/.
package fftpack
import "math"
// Sinqi initializes the array work which is used in both Sinqf and
// Sinqb. 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 which must be dimensioned at least 3*n.
// The same work array can be used for both Sinqf and Sinqb
// 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 Sinqf or Sinqb.
//
// ifac An integer work array of length at least 15.
func Sinqi(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)
}
// Sinqf computes the Fast Fourier Transform of quarter wave data.
// That is, Sinqf computes the coefficients in a sine series
// representation with only odd wave numbers. The transform is
// defined below at output parameter x.
//
// Sinqb is the unnormalized inverse of Sinqf since a call of Sinqf
// followed by a call of Sinqb will multiply the input sequence x
// by 4*n.
//
// The array work which is used by subroutine Sinqf must be
// initialized by calling subroutine Sinqi(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 Sinqf. The work array must be
// initialized by calling subroutine Sinqi(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] = (-1)^(i)*x[n-1]
// + the sum from k=0 to k=n-2 of
// 2*x[k]*sin((2*i+1)*k*pi/(2*n))
//
// A call of Sinqf followed by a call of
// Sinqb will multiply the sequence x by 4*n.
// Therefore Sinqb is the unnormalized inverse
// of Sinqf.
//
// work Contains initialization calculations which must not
// be destroyed between calls of Sinqf or Sinqb.
func Sinqf(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 == 1 {
return
}
for k := 0; k < n/2; k++ {
kc := n - k - 1
x[k], x[kc] = x[kc], x[k]
}
Cosqf(n, x, work, ifac)
for k := 1; k < n; k += 2 {
x[k] = -x[k]
}
}
// Sinqb computes the Fast Fourier Transform of quarter wave data.
// That is, Sinqb computes a sequence from its representation in
// terms of a sine series with odd wave numbers. The transform is
// defined below at output parameter x.
//
// Sinqf is the unnormalized inverse of Sinqb since a call of Sinqb
// followed by a call of Sinqf will multiply the input sequence x
// by 4*n.
//
// The array work which is used by subroutine Sinqb must be
// initialized by calling subroutine Sinqi(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 Sinqb. The work array must be
// initialized by calling subroutine Sinqi(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]*sin((2*k+1)*i*pi/(2*n))
//
// A call of Sinqb followed by a call of
// Sinqf will multiply the sequence x by 4*n.
// Therefore Sinqf is the unnormalized inverse
// of Sinqb.
//
// work Contains initialization calculations which must not
// be destroyed between calls of Sinqb or Sinqf.
func Sinqb(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")
}
switch n {
case 1:
x[0] *= 4
fallthrough
case 0:
return
default:
for k := 1; k < n; k += 2 {
x[k] = -x[k]
}
Cosqb(n, x, work, ifac)
for k := 0; k < n/2; k++ {
kc := n - k - 1
x[k], x[kc] = x[kc], x[k]
}
}
}
|
