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