From b5003e0c532bbc0fefefc710036378ce4c399324 Mon Sep 17 00:00:00 2001 From: azamora Date: Wed, 25 Mar 2026 16:34:18 -0700 Subject: [PATCH 1/4] adding dependency --- go.mod | 1 + go.sum | 4 ++++ 2 files changed, 5 insertions(+) diff --git a/go.mod b/go.mod index dcd7e7e..25ee6f5 100644 --- a/go.mod +++ b/go.mod @@ -22,6 +22,7 @@ require ( github.com/stretchr/testify v1.10.0 github.com/yuin/gopher-lua v1.1.1 golang.org/x/oauth2 v0.26.0 + gonum.org/v1/gonum v0.14.0 gopkg.in/urfave/cli.v1 v1.20.0 gopkg.in/yaml.v2 v2.4.0 layeh.com/gopher-luar v1.0.11 diff --git a/go.sum b/go.sum index 7e54ca9..7540a6c 100644 --- a/go.sum +++ b/go.sum @@ -109,6 +109,8 @@ github.com/yuin/gopher-lua v1.1.1/go.mod h1:GBR0iDaNXjAgGg9zfCvksxSRnQx76gclCIb7 golang.org/x/crypto v0.0.0-20190308221718-c2843e01d9a2/go.mod h1:djNgcEr1/C05ACkg1iLfiJU5Ep61QUkGW8qpdssI0+w= golang.org/x/crypto v0.0.0-20210921155107-089bfa567519/go.mod h1:GvvjBRRGRdwPK5ydBHafDWAxML/pGHZbMvKqRZ5+Abc= golang.org/x/crypto v0.14.0/go.mod h1:MVFd36DqK4CsrnJYDkBA3VC4m2GkXAM0PvzMCn4JQf4= +golang.org/x/exp v0.0.0-20230321023759-10a507213a29 h1:ooxPy7fPvB4kwsA2h+iBNHkAbp/4JxTSwCmvdjEYmug= +golang.org/x/exp v0.0.0-20230321023759-10a507213a29/go.mod h1:CxIveKay+FTh1D0yPZemJVgC/95VzuuOLq5Qi4xnoYc= golang.org/x/mod v0.6.0-dev.0.20220419223038-86c51ed26bb4/go.mod h1:jJ57K6gSWd91VN4djpZkiMVwK6gcyfeH4XE8wZrZaV4= golang.org/x/mod v0.8.0/go.mod h1:iBbtSCu2XBx23ZKBPSOrRkjjQPZFPuis4dIYUhu/chs= golang.org/x/net v0.0.0-20190311183353-d8887717615a/go.mod h1:t9HGtf8HONx5eT2rtn7q6eTqICYqUVnKs3thJo3Qplg= @@ -161,6 +163,8 @@ golang.org/x/tools v0.1.12/go.mod h1:hNGJHUnrk76NpqgfD5Aqm5Crs+Hm0VOH/i9J2+nxYbc golang.org/x/tools v0.6.0/go.mod h1:Xwgl3UAJ/d3gWutnCtw505GrjyAbvKui8lOU390QaIU= golang.org/x/xerrors v0.0.0-20190717185122-a985d3407aa7/go.mod h1:I/5z698sn9Ka8TeJc9MKroUUfqBBauWjQqLJ2OPfmY0= golang.org/x/xerrors v0.0.0-20191204190536-9bdfabe68543/go.mod h1:I/5z698sn9Ka8TeJc9MKroUUfqBBauWjQqLJ2OPfmY0= +gonum.org/v1/gonum v0.14.0 h1:2NiG67LD1tEH0D7kM+ps2V+fXmsAnpUeec7n8tcr4S0= +gonum.org/v1/gonum v0.14.0/go.mod h1:AoWeoz0becf9QMWtE8iWXNXc27fK4fNeHNf/oMejGfU= gopkg.in/check.v1 v0.0.0-20161208181325-20d25e280405/go.mod h1:Co6ibVJAznAaIkqp8huTwlJQCZ016jof/cbN4VW5Yz0= gopkg.in/check.v1 v1.0.0-20201130134442-10cb98267c6c h1:Hei/4ADfdWqJk1ZMxUNpqntNwaWcugrBjAiHlqqRiVk= gopkg.in/check.v1 v1.0.0-20201130134442-10cb98267c6c/go.mod h1:JHkPIbrfpd72SG/EVd6muEfDQjcINNoR0C8j2r3qZ4Q= From 71c77187239d0903fec7e97f097a3e40fce50903 Mon Sep 17 00:00:00 2001 From: azamora Date: Wed, 25 Mar 2026 16:35:31 -0700 Subject: [PATCH 2/4] feat: adding jitter as gogen validator --- jitter/README.md | 131 ++++++++++++++++ jitter/jitter.go | 338 ++++++++++++++++++++++++++++++++++++++++++ jitter/jitter_test.go | 159 ++++++++++++++++++++ main.go | 53 +++++++ 4 files changed, 681 insertions(+) create mode 100644 jitter/README.md create mode 100644 jitter/jitter.go create mode 100644 jitter/jitter_test.go diff --git a/jitter/README.md b/jitter/README.md new file mode 100644 index 0000000..fe07d5a --- /dev/null +++ b/jitter/README.md @@ -0,0 +1,131 @@ +# Jitter Analysis Package + +This package provides statistical analysis of inter-arrival time jitter for event timing validation. It's designed to detect whether event timing has appropriate randomness characteristics or if it exhibits unwanted patterns like periodicity or insufficient variance. + +## Usage + +### CLI Command + +```bash +# Analyze intervals from stdin +cat intervals.txt | gogen jitter + +# With custom thresholds +cat intervals.txt | gogen jitter --cv-threshold 0.6 --burst-threshold -0.4 + +# JSON output +cat intervals.txt | gogen jitter --json + +# Quiet mode (only pass/fail) +cat intervals.txt | gogen jitter --quiet +``` + +### As a Library + +```go +import "github.com/coccyx/gogen/jitter" + +intervals := []float64{1.5, 2.1, 0.8, 1.2, 3.4, 1.9, 2.5, 1.1, 0.9, 2.8} + +opts := jitter.Options{ + CVThreshold: 0.5, + BurstThreshold: -0.5, + SCRThreshold: 0.3, +} + +result, err := jitter.AnalyzeIntervals(intervals, opts) +if err != nil { + // Handle error +} + +if result.OverallPass { + // Intervals have good variance characteristics +} +``` + +## Metrics + +The analysis computes several statistical measures: + +### Coefficient of Variation (CV) +- **Formula**: σ/μ (standard deviation / mean) +- **Interpretation**: Measures relative variability +- **Target**: > 0.5 (for exponential-like distribution, CV ≈ 1.0) + +### Burstiness +- **Formula**: (σ - μ) / (σ + μ) +- **Interpretation**: Characterizes temporal pattern + - B > 0: Bursty (events clustered) + - B ≈ 0: Random (exponential-like) + - B < 0: Regular (periodic-like) +- **Target**: > -0.5 (avoid overly regular patterns) + +### Kolmogorov-Smirnov (KS) Test +- **Purpose**: Tests if intervals follow exponential distribution +- **Interpretation**: Informational only (reported but not pass/fail) + +### Spectral Concentration Ratio (SCR) +- **Method**: FFT-based periodicity detection +- **Formula**: peak_power / total_power +- **Interpretation**: + - High SCR: Energy concentrated at one frequency (periodic) + - Low SCR: Energy spread across frequencies (random) +- **Target**: < 0.3 (avoid strong periodicity) + +## Input Format + +The tool reads inter-arrival times (in seconds) from stdin, one value per line: + +``` +1.5 +2.1 +0.8 +1.2 +... +``` + +Minimum 5 intervals required; 16+ recommended for FFT analysis. + +## Exit Codes + +- 0: All tests passed +- 1: One or more tests failed +- 2: Insufficient data or parsing error + +## Flags + +- `--cv-threshold`: Minimum CV threshold (default: 0.5) +- `--burst-threshold`: Minimum burstiness threshold (default: -0.5) +- `--scr-threshold`: Maximum SCR threshold (default: 0.3) +- `--json`: Output results as JSON +- `--quiet`, `-q`: Only output PASS/FAIL + +## Example Output + +``` +========================================== +Jitter Analysis +========================================== + +Samples analyzed: 100 +Mean interval: 0.9775 s +Std deviation: 0.9642 s +Lambda (MLE): 1.0230 events/s + +CV (target > 0.5): 0.9864 ✓ +Burstiness (target > -0.5): -0.0069 ✓ + +KS Statistic (D): 0.089234 (informational) +KS Critical (95%): 0.136000 ✓ + +Periodicity Analysis (FFT): + SCR (target < 0.3): 0.0123 ✓ + Peak frequency: 0.0204 Hz + Dominant period: 49.0123 s + +RESULT: PASS - Event timing has good variance characteristics +``` + +## Background + +This tool is useful for validating that generated events have realistic timing patterns, especially for security testing scenarios where detection systems might flag perfectly periodic or deterministic event timing as suspicious. diff --git a/jitter/jitter.go b/jitter/jitter.go new file mode 100644 index 0000000..f3a4eca --- /dev/null +++ b/jitter/jitter.go @@ -0,0 +1,338 @@ +package jitter + +import ( + "bufio" + "encoding/json" + "fmt" + "io" + "math" + "os" + "sort" + "strconv" + + "gonum.org/v1/gonum/dsp/fourier" +) + +// Stats contains basic statistical measures for interval analysis +type Stats struct { + N int `json:"n"` + Mean float64 `json:"mean"` + StdDev float64 `json:"stddev"` + Lambda float64 `json:"lambda"` + CV float64 `json:"cv"` + Burstiness float64 `json:"burstiness"` +} + +// KSResult contains Kolmogorov-Smirnov test results +type KSResult struct { + Statistic float64 `json:"ks_stat"` + Critical float64 `json:"ks_critical"` + Pass bool `json:"ks_pass"` +} + +// FFTResult contains Fast Fourier Transform periodicity analysis results +type FFTResult struct { + SCR float64 `json:"scr"` + PeakFreq float64 `json:"peak_freq"` + Period float64 `json:"period"` + Pass bool `json:"fft_pass"` +} + +// TestResult contains all analysis results +type TestResult struct { + Stats Stats `json:"stats"` + KS KSResult `json:"ks"` + FFT FFTResult `json:"fft"` + CVPass bool `json:"cv_pass"` + BurstPass bool `json:"burst_pass"` + OverallPass bool `json:"pass"` +} + +// Options configures the jitter analysis +type Options struct { + CVThreshold float64 + BurstThreshold float64 + SCRThreshold float64 + JSONOutput bool + Quiet bool + Reader io.Reader +} + +// Analyze performs jitter analysis on intervals read from the provided reader +// Returns an error if there's a problem reading/parsing data, or if tests fail. +func Analyze(opts Options) error { + // Read intervals from input + intervals, err := readIntervals(opts.Reader) + if err != nil { + return err + } + + result, err := AnalyzeIntervals(intervals, opts) + if err != nil { + return err + } + + // Output results + if opts.JSONOutput { + if err := printJSON(result); err != nil { + return err + } + } else if opts.Quiet { + if result.OverallPass { + fmt.Println("PASS") + } else { + fmt.Println("FAIL") + } + } else { + printResults(result, opts.CVThreshold, opts.BurstThreshold, opts.SCRThreshold) + } + + if !result.OverallPass { + return fmt.Errorf("jitter analysis failed") + } + return nil +} + +// AnalyzeIntervals performs jitter analysis on the provided intervals +func AnalyzeIntervals(intervals []float64, opts Options) (TestResult, error) { + n := len(intervals) + if n < 5 { + return TestResult{}, fmt.Errorf("need at least 5 intervals, got %d", n) + } + + // Calculate basic statistics + stats := calculateStats(intervals) + + // Calculate KS statistic + ks := calculateKS(intervals, stats.Lambda) + + // Calculate FFT periodicity + fft := calculateFFT(intervals, stats.Mean, opts.SCRThreshold) + + // Determine pass/fail + cvPass := stats.CV > opts.CVThreshold + burstPass := stats.Burstiness > opts.BurstThreshold + overallPass := cvPass && burstPass && fft.Pass + + result := TestResult{ + Stats: stats, + KS: ks, + FFT: fft, + CVPass: cvPass, + BurstPass: burstPass, + OverallPass: overallPass, + } + + return result, nil +} + +func readIntervals(r io.Reader) ([]float64, error) { + var intervals []float64 + scanner := bufio.NewScanner(r) + for scanner.Scan() { + val, err := strconv.ParseFloat(scanner.Text(), 64) + if err == nil && val > 0 { + intervals = append(intervals, val) + } + } + if err := scanner.Err(); err != nil { + return nil, err + } + return intervals, nil +} + +func calculateStats(intervals []float64) Stats { + n := len(intervals) + + var sum, sumSq float64 + for _, v := range intervals { + sum += v + sumSq += v * v + } + + mean := sum / float64(n) + variance := (sumSq / float64(n)) - (mean * mean) + if variance < 0 { + variance = 0 + } + stddev := math.Sqrt(variance) + + // Rate parameter lambda = 1/mean + lambda := 0.0 + if mean > 0 { + lambda = 1.0 / mean + } + + // CV (coefficient of variation) = stddev / mean + // For exponential distribution, CV should be ~1.0 + cv := 0.0 + if mean > 0 { + cv = stddev / mean + } + + // Burstiness B = (sigma - mu) / (sigma + mu) + // For exponential, B ≈ 0; for bursty B > 0; for regular B < 0 + burstiness := 0.0 + if stddev+mean > 0 { + burstiness = (stddev - mean) / (stddev + mean) + } + + return Stats{ + N: n, + Mean: mean, + StdDev: stddev, + Lambda: lambda, + CV: cv, + Burstiness: burstiness, + } +} + +func calculateKS(intervals []float64, lambda float64) KSResult { + n := len(intervals) + + // Sort intervals for empirical CDF + sorted := make([]float64, n) + copy(sorted, intervals) + sort.Float64s(sorted) + + // Compute KS statistic: D = max|F_n(x) - F(x)| + // where F(x) = 1 - e^(-lambda*x) for exponential + maxD := 0.0 + for i, x := range sorted { + // Empirical CDF at x + fn := float64(i+1) / float64(n) + fnPrev := float64(i) / float64(n) + + // Theoretical CDF for exponential + fx := 0.0 + if lambda > 0 && x > 0 { + fx = 1.0 - math.Exp(-lambda*x) + } + + // Check both D+ and D- + d1 := math.Abs(fn - fx) + d2 := math.Abs(fnPrev - fx) + if d1 > maxD { + maxD = d1 + } + if d2 > maxD { + maxD = d2 + } + } + + // Critical value at 95% confidence: D_crit = 1.36 / sqrt(n) + critical := 1.36 / math.Sqrt(float64(n)) + + return KSResult{ + Statistic: maxD, + Critical: critical, + Pass: maxD < critical, + } +} + +func calculateFFT(intervals []float64, meanInterval float64, threshold float64) FFTResult { + n := len(intervals) + + if n < 16 { + return FFTResult{SCR: 0, PeakFreq: 0, Period: 0, Pass: true} + } + + // Compute FFT using gonum + // For real input of length N, returns N/2+1 complex coefficients + fft := fourier.NewFFT(n) + coeffs := fft.Coefficients(nil, intervals) + + // Compute power spectrum (exclude DC component at index 0) + // Power at each frequency bin = |coefficient|² = real² + imag² + var totalPower, peakPower float64 + var peakIdx int + for i := 1; i < len(coeffs); i++ { + power := real(coeffs[i])*real(coeffs[i]) + imag(coeffs[i])*imag(coeffs[i]) + totalPower += power + if power > peakPower { + peakPower = power + peakIdx = i + } + } + + // SCR (Spectral Concentration Ratio) = peak_power / total_power + // High SCR = energy concentrated at one frequency (periodic) + // Low SCR = energy spread across frequencies (aperiodic/random) + scr := 0.0 + if totalPower > 0 { + scr = peakPower / totalPower + } + + // Convert FFT bin index to physical frequency and period + // Frequency: f = k / (N × mean_interval) Hz + peakFreq := float64(peakIdx) / (float64(n) * meanInterval) + + // Period = 1 / frequency (in seconds) + period := 0.0 + if peakFreq > 0 { + period = 1.0 / peakFreq + } + + return FFTResult{ + SCR: scr, + PeakFreq: peakFreq, + Period: period, + Pass: scr < threshold, + } +} + +func printResults(r TestResult, cvThresh, burstThresh, scrThresh float64) { + fmt.Println() + fmt.Println("==========================================") + fmt.Println("Jitter Analysis") + fmt.Println("==========================================") + fmt.Println() + fmt.Printf("Samples analyzed: %d\n", r.Stats.N) + fmt.Printf("Mean interval: %.4f s\n", r.Stats.Mean) + fmt.Printf("Std deviation: %.4f s\n", r.Stats.StdDev) + fmt.Printf("Lambda (MLE): %.4f events/s\n", r.Stats.Lambda) + fmt.Println() + + cvMark := mark(r.CVPass) + burstMark := mark(r.BurstPass) + ksMark := "" + if r.KS.Pass { + ksMark = "✓" + } + + fmt.Printf("CV (target > %.1f): %.4f %s\n", cvThresh, r.Stats.CV, cvMark) + fmt.Printf("Burstiness (target > %.1f): %.4f %s\n", burstThresh, r.Stats.Burstiness, burstMark) + fmt.Println() + fmt.Printf("KS Statistic (D): %.6f (informational)\n", r.KS.Statistic) + fmt.Printf("KS Critical (95%%): %.6f %s\n", r.KS.Critical, ksMark) + fmt.Println() + + fmt.Println("Periodicity Analysis (FFT):") + fftMark := mark(r.FFT.Pass) + fmt.Printf(" SCR (target < %.1f): %.4f %s\n", scrThresh, r.FFT.SCR, fftMark) + fmt.Printf(" Peak frequency: %.4f Hz\n", r.FFT.PeakFreq) + fmt.Printf(" Dominant period: %.4f s\n", r.FFT.Period) + fmt.Println() + + if r.OverallPass { + fmt.Println("RESULT: PASS - Event timing has good variance characteristics") + } else { + fmt.Println("RESULT: FAIL - Event timing lacks sufficient variance or has periodicity") + } + fmt.Println() +} + +func printJSON(r TestResult) error { + enc := json.NewEncoder(os.Stdout) + if err := enc.Encode(r); err != nil { + return err + } + return nil +} + +func mark(pass bool) string { + if pass { + return "✓" + } + return "✗" +} diff --git a/jitter/jitter_test.go b/jitter/jitter_test.go new file mode 100644 index 0000000..1121be4 --- /dev/null +++ b/jitter/jitter_test.go @@ -0,0 +1,159 @@ +package jitter + +import ( + "bytes" + "fmt" + "strings" + "testing" +) + +func TestAnalyze_HighVariance(t *testing.T) { + // Test with data that has high coefficient of variation + data := "0.05\n4.2\n0.1\n5.3\n0.2\n3.8\n0.15\n6.1\n0.3\n2.9\n" + + "0.08\n4.7\n0.25\n3.2\n0.12\n5.5\n0.18\n4.1\n0.22\n3.5\n" + + intervals := []float64{} + for _, line := range strings.Split(strings.TrimSpace(data), "\n") { + var val float64 + if _, err := fmt.Sscanf(line, "%f", &val); err == nil { + intervals = append(intervals, val) + } + } + + opts := Options{ + CVThreshold: 0.5, + BurstThreshold: -0.5, + SCRThreshold: 0.3, + } + + result, err := AnalyzeIntervals(intervals, opts) + if err != nil { + t.Fatalf("Expected analysis to complete, got error: %v", err) + } + + // Verify the result has all expected fields + if result.Stats.N != len(intervals) { + t.Errorf("Expected N=%d, got %d", len(intervals), result.Stats.N) + } + + if result.Stats.Mean <= 0 { + t.Error("Expected positive mean") + } + + if result.Stats.CV <= 0 { + t.Error("Expected positive CV") + } +} + +func TestAnalyze_JSONOutput(t *testing.T) { + data := "1.5\n2.1\n0.8\n1.2\n3.4\n1.9\n2.5\n1.1\n0.9\n2.8\n" + reader := strings.NewReader(data) + + opts := Options{ + CVThreshold: 0.5, + BurstThreshold: -0.5, + SCRThreshold: 0.3, + JSONOutput: true, + Quiet: false, + Reader: reader, + } + + // Capture stderr (we'll ignore the exit behavior in tests) + var buf bytes.Buffer + + // This will likely fail but should produce JSON output + // We're just testing that it doesn't panic + _ = Analyze(opts) + + // Test passed if we got here without panic + _ = buf.String() +} + +func TestAnalyze_InsufficientData(t *testing.T) { + data := "1.5\n2.1\n" + reader := strings.NewReader(data) + + opts := Options{ + CVThreshold: 0.5, + BurstThreshold: -0.5, + SCRThreshold: 0.3, + JSONOutput: false, + Quiet: true, + Reader: reader, + } + + err := Analyze(opts) + if err == nil { + t.Fatal("Expected error for insufficient data, got nil") + } + + if !strings.Contains(err.Error(), "at least 5 intervals") { + t.Errorf("Expected error about insufficient intervals, got: %v", err) + } +} + +func TestCalculateStats(t *testing.T) { + intervals := []float64{1.0, 2.0, 3.0, 4.0, 5.0} + + stats := calculateStats(intervals) + + if stats.N != 5 { + t.Errorf("Expected N=5, got %d", stats.N) + } + + if stats.Mean != 3.0 { + t.Errorf("Expected Mean=3.0, got %f", stats.Mean) + } + + // Check CV is calculated + if stats.CV <= 0 { + t.Errorf("Expected positive CV, got %f", stats.CV) + } +} + +func TestCalculateKS(t *testing.T) { + intervals := []float64{1.0, 2.0, 3.0, 4.0, 5.0} + lambda := 1.0 / 3.0 // Mean of intervals is 3.0 + + ks := calculateKS(intervals, lambda) + + if ks.Statistic < 0 || ks.Statistic > 1 { + t.Errorf("KS statistic should be between 0 and 1, got %f", ks.Statistic) + } + + if ks.Critical <= 0 { + t.Errorf("KS critical value should be positive, got %f", ks.Critical) + } +} + +func TestCalculateFFT(t *testing.T) { + // Create 20 intervals for FFT analysis + intervals := make([]float64, 20) + for i := range intervals { + intervals[i] = float64(i%3 + 1) // Simple varying pattern + } + + meanInterval := 2.0 + threshold := 0.3 + + fft := calculateFFT(intervals, meanInterval, threshold) + + // SCR should be between 0 and 1 + if fft.SCR < 0 || fft.SCR > 1 { + t.Errorf("SCR should be between 0 and 1, got %f", fft.SCR) + } +} + +func TestCalculateFFT_SmallDataset(t *testing.T) { + // FFT should handle small datasets gracefully + intervals := []float64{1.0, 2.0, 3.0} + meanInterval := 2.0 + threshold := 0.3 + + fft := calculateFFT(intervals, meanInterval, threshold) + + // Should pass with small dataset + if !fft.Pass { + t.Error("Small dataset should automatically pass FFT test") + } +} diff --git a/main.go b/main.go index 134e6c1..a9b2eac 100644 --- a/main.go +++ b/main.go @@ -13,6 +13,7 @@ import ( "time" config "github.com/coccyx/gogen/internal" + "github.com/coccyx/gogen/jitter" log "github.com/coccyx/gogen/logger" "github.com/coccyx/gogen/run" "github.com/olekukonko/tablewriter" @@ -533,6 +534,58 @@ func main() { return nil }, }, + { + Name: "jitter", + Usage: "Analyze inter-arrival time jitter from stdin", + ArgsUsage: " ", + SkipFlagParsing: false, + Flags: []cli.Flag{ + cli.Float64Flag{ + Name: "cv-threshold", + Value: 0.5, + Usage: "Minimum coefficient of variation threshold", + }, + cli.Float64Flag{ + Name: "burst-threshold", + Value: -0.5, + Usage: "Minimum burstiness threshold", + }, + cli.Float64Flag{ + Name: "scr-threshold", + Value: 0.3, + Usage: "Maximum spectral concentration ratio threshold for periodicity", + }, + cli.BoolFlag{ + Name: "json", + Usage: "Output results as JSON", + }, + cli.BoolFlag{ + Name: "quiet, q", + Usage: "Quiet mode - only output pass/fail", + }, + }, + Before: func(clic *cli.Context) error { + // Skip normal Setup for jitter command + return nil + }, + Action: func(clic *cli.Context) error { + opts := jitter.Options{ + CVThreshold: clic.Float64("cv-threshold"), + BurstThreshold: clic.Float64("burst-threshold"), + SCRThreshold: clic.Float64("scr-threshold"), + JSONOutput: clic.Bool("json"), + Quiet: clic.Bool("quiet"), + Reader: os.Stdin, + } + if err := jitter.Analyze(opts); err != nil { + if !clic.Bool("json") && !clic.Bool("quiet") { + fmt.Fprintf(os.Stderr, "Error: %v\n", err) + } + os.Exit(1) + } + return nil + }, + }, } app.Before = func(clic *cli.Context) error { Setup(clic) From f278e62a6827fa20e9b60ed91dee7dfedf286235 Mon Sep 17 00:00:00 2001 From: azamora Date: Wed, 25 Mar 2026 16:47:23 -0700 Subject: [PATCH 3/4] docs: adding dsp modulest to modules.txt --- vendor/modules.txt | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/vendor/modules.txt b/vendor/modules.txt index 01bae51..2f9b9e4 100644 --- a/vendor/modules.txt +++ b/vendor/modules.txt @@ -161,6 +161,10 @@ golang.org/x/oauth2/internal ## explicit; go 1.23.0 golang.org/x/sys/unix golang.org/x/sys/windows +# gonum.org/v1/gonum v0.14.0 +## explicit; go 1.20 +gonum.org/v1/gonum/dsp/fourier +gonum.org/v1/gonum/dsp/fourier/internal/fftpack # gopkg.in/check.v1 v1.0.0-20201130134442-10cb98267c6c ## explicit; go 1.11 # gopkg.in/urfave/cli.v1 v1.20.0 From b58744afb2ff6b925baeba8554e58204fae54fac Mon Sep 17 00:00:00 2001 From: azamora Date: Sat, 28 Mar 2026 10:19:40 -0400 Subject: [PATCH 4/4] Add gonum vendor dependency for jitter package --- vendor/gonum.org/v1/gonum/AUTHORS | 128 ++ vendor/gonum.org/v1/gonum/CONTRIBUTORS | 131 ++ vendor/gonum.org/v1/gonum/LICENSE | 23 + vendor/gonum.org/v1/gonum/dsp/fourier/doc.go | 6 + .../gonum.org/v1/gonum/dsp/fourier/fourier.go | 260 ++++ .../internal/fftpack/array_bounds_checks.go | 119 ++ .../fftpack/array_no_bounds_checks.go | 83 ++ .../dsp/fourier/internal/fftpack/cfft.go | 588 +++++++++ .../dsp/fourier/internal/fftpack/cosq.go | 219 ++++ .../dsp/fourier/internal/fftpack/cost.go | 143 +++ .../gonum/dsp/fourier/internal/fftpack/doc.go | 7 + .../dsp/fourier/internal/fftpack/helpers.go | 15 + .../dsp/fourier/internal/fftpack/rfft.go | 1079 +++++++++++++++++ .../dsp/fourier/internal/fftpack/sinq.go | 179 +++ .../dsp/fourier/internal/fftpack/sint.go | 146 +++ .../gonum.org/v1/gonum/dsp/fourier/quarter.go | 133 ++ .../gonum.org/v1/gonum/dsp/fourier/radix24.go | 343 ++++++ .../gonum.org/v1/gonum/dsp/fourier/sincos.go | 112 ++ 18 files changed, 3714 insertions(+) create mode 100644 vendor/gonum.org/v1/gonum/AUTHORS create mode 100644 vendor/gonum.org/v1/gonum/CONTRIBUTORS create mode 100644 vendor/gonum.org/v1/gonum/LICENSE create mode 100644 vendor/gonum.org/v1/gonum/dsp/fourier/doc.go create mode 100644 vendor/gonum.org/v1/gonum/dsp/fourier/fourier.go create mode 100644 vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/array_bounds_checks.go create mode 100644 vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/array_no_bounds_checks.go create mode 100644 vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/cfft.go create mode 100644 vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/cosq.go create mode 100644 vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/cost.go create mode 100644 vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/doc.go create mode 100644 vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/helpers.go create mode 100644 vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/rfft.go create mode 100644 vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/sinq.go create mode 100644 vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/sint.go create mode 100644 vendor/gonum.org/v1/gonum/dsp/fourier/quarter.go create mode 100644 vendor/gonum.org/v1/gonum/dsp/fourier/radix24.go create mode 100644 vendor/gonum.org/v1/gonum/dsp/fourier/sincos.go diff --git a/vendor/gonum.org/v1/gonum/AUTHORS b/vendor/gonum.org/v1/gonum/AUTHORS new file mode 100644 index 0000000..91467ea --- /dev/null +++ b/vendor/gonum.org/v1/gonum/AUTHORS @@ -0,0 +1,128 @@ +# This is the official list of Gonum authors for copyright purposes. +# This file is distinct from the CONTRIBUTORS files. +# See the latter for an explanation. + +# Names should be added to this file as +# Name or Organization +# The email address is not required for organizations. + +# Please keep the list sorted. + +Alexander Egurnov +Andrei Blinnikov +antichris +Bailey Lissington +Bill Gray +Bill Noon +Brendan Tracey +Brent Pedersen +Bulat Khasanov +Chad Kunde +Chan Kwan Yin +Chih-Wei Chang +Chong-Yeol Nah +Chris Tessum +Christophe Meessen +Christopher Waldon +Clayton Northey +Dan Kortschak +Daniel Fireman +Dario Heinisch +David Kleiven +David Samborski +Davor Kapsa +DeepMind Technologies +Delaney Gillilan +Dezmond Goff +Dong-hee Na +Dustin Spicuzza +Egon Elbre +Ekaterina Efimova +Ethan Burns +Evert Lammerts +Evgeny Savinov +Fabian Wickborn +Facundo Gaich +Fazlul Shahriar +Francesc Campoy +Google Inc +Gustaf Johansson +Hossein Zolfi +Iakov Davydov +Igor Mikushkin +Iskander Sharipov +Jalem Raj Rohit +James Bell +James Bowman +James Holmes <32bitkid@gmail.com> +Janne Snabb +Jeremy Atkinson +Jes Cok +Jinesi Yelizati +Jonas Kahler +Jonas Schulze +Jonathan J Lawlor +Jonathan Reiter +Jonathan Schroeder +Joost van Amersfoort +Joseph Watson +Josh Wilson +Julien Roland +Kai Trukenmüller +Kent English +Kevin C. Zimmerman +Kirill Motkov +Konstantin Shaposhnikov +Leonid Kneller +Lyron Winderbaum +Marco Leogrande +Mark Canning +Mark Skilbeck +Martin Diz +Matthew Connelly +Matthieu Di Mercurio +Max Halford +Maxim Sergeev +Microsoft Corporation +MinJae Kwon +Nathan Edwards +Nick Potts +Nils Wogatzky +Olivier Wulveryck +Or Rikon +Patricio Whittingslow +Patrick DeVivo +Pontus Melke +Renee French +Rishi Desai +Robin Eklind +Roger Welin +Rondall Jones +Sam Zaydel +Samuel Kelemen +Saran Ahluwalia +Scott Holden +Scott Kiesel +Sebastien Binet +Shawn Smith +Sintela Ltd +source{d} +Spencer Lyon +Steve McCoy +Taesu Pyo +Takeshi Yoneda +Tamir Hyman +The University of Adelaide +The University of Minnesota +The University of Washington +Thomas Berg +Tobin Harding +Valentin Deleplace +Vincent Thiery +Vladimír Chalupecký +Will Tekulve +Yasuhiro Matsumoto +Yevgeniy Vahlis +Yucheng Zhu +Yunomi +Zoe Juozapaitis diff --git a/vendor/gonum.org/v1/gonum/CONTRIBUTORS b/vendor/gonum.org/v1/gonum/CONTRIBUTORS new file mode 100644 index 0000000..93e01a2 --- /dev/null +++ b/vendor/gonum.org/v1/gonum/CONTRIBUTORS @@ -0,0 +1,131 @@ +# This is the official list of people who can contribute +# (and typically have contributed) code to the Gonum +# project. +# +# The AUTHORS file lists the copyright holders; this file +# lists people. For example, Google employees would be listed here +# but not in AUTHORS, because Google would hold the copyright. +# +# When adding J Random Contributor's name to this file, +# either J's name or J's organization's name should be +# added to the AUTHORS file. +# +# Names should be added to this file like so: +# Name +# +# Please keep the list sorted. + +Alexander Egurnov +Andrei Blinnikov +Andrew Brampton +antichris +Bailey Lissington +Bill Gray +Bill Noon +Brendan Tracey +Brent Pedersen +Bulat Khasanov +Chad Kunde +Chan Kwan Yin +Chih-Wei Chang +Chong-Yeol Nah +Chris Tessum +Christophe Meessen +Christopher Waldon +Clayton Northey +Dan Kortschak +Dan Lorenc +Daniel Fireman +Dario Heinisch +David Kleiven +David Samborski +Davor Kapsa +Delaney Gillilan +Dezmond Goff +Dong-hee Na +Dustin Spicuzza +Egon Elbre +Ekaterina Efimova +Ethan Burns +Evert Lammerts +Evgeny Savinov +Fabian Wickborn +Facundo Gaich +Fazlul Shahriar +Francesc Campoy +Gustaf Johansson +Hossein Zolfi +Iakov Davydov +Igor Mikushkin +Iskander Sharipov +Jalem Raj Rohit +James Bell +James Bowman +James Holmes <32bitkid@gmail.com> +Janne Snabb +Jeremy Atkinson +Jes Cok +Jinesi Yelizati +Jon Richards +Jonas Kahler +Jonas Schulze +Jonathan J Lawlor +Jonathan Reiter +Jonathan Schroeder +Joost van Amersfoort +Joseph Watson +Josh Wilson +Julien Roland +Kai Trukenmüller +Kent English +Kevin C. Zimmerman +Kirill Motkov +Konstantin Shaposhnikov +Leonid Kneller +Lyron Winderbaum +Marco Leogrande +Mark Canning +Mark Skilbeck +Martin Diz +Matthew Connelly +Matthieu Di Mercurio +Max Halford +Maxim Sergeev +MinJae Kwon +Nathan Edwards +Nick Potts +Nils Wogatzky +Olivier Wulveryck +Or Rikon +Patricio Whittingslow +Patrick DeVivo +Pontus Melke +Renee French +Rishi Desai +Robin Eklind +Roger Welin +Roman Werpachowski +Rondall Jones +Sam Zaydel +Samuel Kelemen +Saran Ahluwalia +Scott Holden +Scott Kiesel +Sebastien Binet +Shawn Smith +Spencer Lyon +Steve McCoy +Taesu Pyo +Takeshi Yoneda +Tamir Hyman +Thomas Berg +Tobin Harding +Valentin Deleplace +Vincent Thiery +Vladimír Chalupecký +Will Tekulve +Yasuhiro Matsumoto +Yevgeniy Vahlis +Yucheng Zhu +Yunomi +Zoe Juozapaitis diff --git a/vendor/gonum.org/v1/gonum/LICENSE b/vendor/gonum.org/v1/gonum/LICENSE new file mode 100644 index 0000000..ed477e5 --- /dev/null +++ b/vendor/gonum.org/v1/gonum/LICENSE @@ -0,0 +1,23 @@ +Copyright ©2013 The Gonum Authors. All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of the Gonum project nor the names of its authors and + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. \ No newline at end of file diff --git a/vendor/gonum.org/v1/gonum/dsp/fourier/doc.go b/vendor/gonum.org/v1/gonum/dsp/fourier/doc.go new file mode 100644 index 0000000..a33e7e7 --- /dev/null +++ b/vendor/gonum.org/v1/gonum/dsp/fourier/doc.go @@ -0,0 +1,6 @@ +// 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. + +// Package fourier provides functions to perform Discrete Fourier Transforms. +package fourier // import "gonum.org/v1/gonum/dsp/fourier" diff --git a/vendor/gonum.org/v1/gonum/dsp/fourier/fourier.go b/vendor/gonum.org/v1/gonum/dsp/fourier/fourier.go new file mode 100644 index 0000000..1ac7dcc --- /dev/null +++ b/vendor/gonum.org/v1/gonum/dsp/fourier/fourier.go @@ -0,0 +1,260 @@ +// 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. + +package fourier + +import "gonum.org/v1/gonum/dsp/fourier/internal/fftpack" + +// FFT implements Fast Fourier Transform and its inverse for real sequences. +type FFT struct { + work []float64 + ifac [15]int + + // real temporarily store complex data as + // pairs of real values to allow passing to + // the backing code. The length of real + // must always be half the length of work. + real []float64 +} + +// NewFFT returns an FFT initialized for work on sequences of length n. +func NewFFT(n int) *FFT { + var t FFT + t.Reset(n) + return &t +} + +// Len returns the length of the acceptable input. +func (t *FFT) Len() int { return len(t.real) } + +// Reset reinitializes the FFT for work on sequences of length n. +func (t *FFT) Reset(n int) { + if 2*n <= cap(t.work) { + t.work = t.work[:2*n] + t.real = t.real[:n] + } else { + t.work = make([]float64, 2*n) + t.real = make([]float64, n) + } + fftpack.Rffti(n, t.work, t.ifac[:]) +} + +// Coefficients computes the Fourier coefficients of the input sequence, +// converting the time series in seq into the frequency spectrum, placing +// the result in dst and returning it. This transform is unnormalized; a +// call to Coefficients followed by a call of Sequence will multiply the +// input sequence by the length of the sequence. +// +// If the length of seq is not t.Len(), Coefficients will panic. +// If dst is nil, a new slice is allocated and returned. If dst is not nil and +// the length of dst does not equal t.Len()/2+1, Coefficients will panic. +func (t *FFT) Coefficients(dst []complex128, seq []float64) []complex128 { + if len(seq) != t.Len() { + panic("fourier: sequence length mismatch") + } + if dst == nil { + dst = make([]complex128, t.Len()/2+1) + } else if len(dst) != t.Len()/2+1 { + panic("fourier: destination length mismatch") + } + copy(t.real, seq) + fftpack.Rfftf(len(t.real), t.real, t.work, t.ifac[:]) + dst[0] = complex(t.real[0], 0) + if len(seq) < 2 { + return dst + } + if len(seq)%2 == 1 { + dst[len(dst)-1] = complex(t.real[len(t.real)-2], t.real[len(t.real)-1]) + } else { + dst[len(dst)-1] = complex(t.real[len(t.real)-1], 0) + } + for i := 1; i < len(dst)-1; i++ { + dst[i] = complex(t.real[2*i-1], t.real[2*i]) + } + return dst +} + +// Sequence computes the real periodic sequence from the Fourier coefficients, +// converting the frequency spectrum in coeff into a time series, placing the +// result in dst and returning it. This transform is unnormalized; a call to +// Coefficients followed by a call of Sequence will multiply the input sequence +// by the length of the sequence. +// +// If the length of coeff is not t.Len()/2+1, Sequence will panic. +// If dst is nil, a new slice is allocated and returned. If dst is not nil and +// the length of dst does not equal the length of coeff, Sequence will panic. +func (t *FFT) Sequence(dst []float64, coeff []complex128) []float64 { + if len(coeff) != t.Len()/2+1 { + panic("fourier: coefficients length mismatch") + } + if dst == nil { + dst = make([]float64, t.Len()) + } else if len(dst) != t.Len() { + panic("fourier: destination length mismatch") + } + dst[0] = real(coeff[0]) + if len(dst) < 2 { + return dst + } + nf := coeff[len(coeff)-1] + if len(dst)%2 == 1 { + dst[len(dst)-2] = real(nf) + dst[len(dst)-1] = imag(nf) + } else { + dst[len(dst)-1] = real(nf) + } + + for i, cv := range coeff[1 : len(coeff)-1] { + dst[2*i+1] = real(cv) + dst[2*i+2] = imag(cv) + } + fftpack.Rfftb(len(dst), dst, t.work, t.ifac[:]) + return dst +} + +// Freq returns the relative frequency center for coefficient i. +// Freq will panic if i is negative or greater than or equal to t.Len(). +func (t *FFT) Freq(i int) float64 { + if i < 0 || t.Len() <= i { + panic("fourier: index out of range") + } + step := 1 / float64(t.Len()) + return step * float64(i) +} + +// CmplxFFT implements Fast Fourier Transform and its inverse for complex sequences. +type CmplxFFT struct { + work []float64 + ifac [15]int + + // real temporarily store complex data as + // pairs of real values to allow passing to + // the backing code. The length of real + // must always be half the length of work. + real []float64 +} + +// NewCmplxFFT returns an CmplxFFT initialized for work on sequences of length n. +func NewCmplxFFT(n int) *CmplxFFT { + var t CmplxFFT + t.Reset(n) + return &t +} + +// Len returns the length of the acceptable input. +func (t *CmplxFFT) Len() int { return len(t.work) / 4 } + +// Reset reinitializes the FFT for work on sequences of length n. +func (t *CmplxFFT) Reset(n int) { + if 4*n <= cap(t.work) { + t.work = t.work[:4*n] + t.real = t.real[:2*n] + } else { + t.work = make([]float64, 4*n) + t.real = make([]float64, 2*n) + } + fftpack.Cffti(n, t.work, t.ifac[:]) +} + +// Coefficients computes the Fourier coefficients of a complex input sequence, +// converting the time series in seq into the frequency spectrum, placing +// the result in dst and returning it. This transform is unnormalized; a call +// to Coefficients followed by a call of Sequence will multiply the input +// sequence by the length of the sequence. +// +// If the length of seq is not t.Len(), Coefficients will panic. +// If dst is nil, a new slice is allocated and returned. If dst is not nil and +// the length of dst does not equal the length of seq, Coefficients will panic. +// It is safe to use the same slice for dst and seq. +func (t *CmplxFFT) Coefficients(dst, seq []complex128) []complex128 { + if len(seq) != t.Len() { + panic("fourier: sequence length mismatch") + } + if dst == nil { + dst = make([]complex128, len(seq)) + } else if len(dst) != len(seq) { + panic("fourier: destination length mismatch") + } + for i, cv := range seq { + t.real[2*i] = real(cv) + t.real[2*i+1] = imag(cv) + } + fftpack.Cfftf(len(dst), t.real, t.work, t.ifac[:]) + for i := range dst { + dst[i] = complex(t.real[2*i], t.real[2*i+1]) + } + return dst +} + +// Sequence computes the complex periodic sequence from the Fourier coefficients, +// converting the frequency spectrum in coeff into a time series, placing the +// result in dst and returning it. This transform is unnormalized; a call to +// Coefficients followed by a call of Sequence will multiply the input sequence +// by the length of the sequence. +// +// If the length of coeff is not t.Len(), Sequence will panic. +// If dst is nil, a new slice is allocated and returned. If dst is not nil and +// the length of dst does not equal the length of coeff, Sequence will panic. +// It is safe to use the same slice for dst and coeff. +func (t *CmplxFFT) Sequence(dst, coeff []complex128) []complex128 { + if len(coeff) != t.Len() { + panic("fourier: coefficients length mismatch") + } + if dst == nil { + dst = make([]complex128, len(coeff)) + } else if len(dst) != len(coeff) { + panic("fourier: destination length mismatch") + } + for i, cv := range coeff { + t.real[2*i] = real(cv) + t.real[2*i+1] = imag(cv) + } + fftpack.Cfftb(len(dst), t.real, t.work, t.ifac[:]) + for i := range dst { + dst[i] = complex(t.real[2*i], t.real[2*i+1]) + } + return dst +} + +// Freq returns the relative frequency center for coefficient i. +// Freq will panic if i is negative or greater than or equal to t.Len(). +func (t *CmplxFFT) Freq(i int) float64 { + if i < 0 || t.Len() <= i { + panic("fourier: index out of range") + } + step := 1 / float64(t.Len()) + if i < (t.Len()-1)/2+1 { + return step * float64(i) + } + return step * float64(i-t.Len()) +} + +// ShiftIdx returns a shifted index into a slice of coefficients +// returned by the CmplxFFT so that indexing into the coefficients +// places the zero frequency component at the center of the spectrum. +// ShiftIdx will panic if i is negative or greater than or equal to +// t.Len(). +func (t *CmplxFFT) ShiftIdx(i int) int { + if i < 0 || t.Len() <= i { + panic("fourier: index out of range") + } + h := t.Len() / 2 + if i < h { + return i + (t.Len()+1)/2 + } + return i - h +} + +// UnshiftIdx returns inverse of ShiftIdx. UnshiftIdx will panic if i is +// negative or greater than or equal to t.Len(). +func (t *CmplxFFT) UnshiftIdx(i int) int { + if i < 0 || t.Len() <= i { + panic("fourier: index out of range") + } + h := (t.Len() + 1) / 2 + if i < h { + return i + t.Len()/2 + } + return i - h +} diff --git a/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/array_bounds_checks.go b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/array_bounds_checks.go new file mode 100644 index 0000000..b53b817 --- /dev/null +++ b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/array_bounds_checks.go @@ -0,0 +1,119 @@ +// 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 file must be kept in sync with array_no_bound_checks.go. + +//go:build bounds +// +build bounds + +package fftpack + +import "fmt" + +// The types in array.go implement Fortran-like arrays for bootstrapping +// the implementation of the FFT functions translated from FFTPACK; they +// are column-major. + +type twoArray struct { + i, j int + jStride int + data []float64 +} + +func newTwoArray(i, j int, data []float64) twoArray { + if len(data) < i*j { + panic(fmt.Sprintf("short data: len(data)=%d, i=%d, j=%d", len(data), i, j)) + } + return twoArray{ + i: i, + j: j, + jStride: i, + data: data[:i*j], + } +} + +func (a twoArray) at(i, j int) float64 { + if i < 0 || a.i <= i || j < 0 || a.j <= j { + panic(fmt.Sprintf("out of bounds at(%d, %d): bounds i=%d, j=%d", i, j, a.i, a.j)) + } + return a.data[i+a.jStride*j] +} + +func (a twoArray) atCmplx(i, j int) complex128 { + if i < 0 || a.i <= i || j < 0 || a.j <= j { + panic(fmt.Sprintf("out of bounds at(%d, %d): bounds i=%d, j=%d", i, j, a.i, a.j)) + } + return complex(a.data[i+a.jStride*j], a.data[i+a.jStride*j+1]) +} + +func (a twoArray) set(i, j int, v float64) { + if i < 0 || a.i <= i || j < 0 || a.j <= j { + panic(fmt.Sprintf("out of bounds set(%d, %d): bounds i=%d, j=%d", i, j, a.i, a.j)) + } + a.data[i+a.jStride*j] = v +} + +func (a twoArray) setCmplx(i, j int, v complex128) { + if i < 0 || a.i <= i || j < 0 || a.j <= j { + panic(fmt.Sprintf("out of bounds set(%d, %d): bounds i=%d, j=%d", i, j, a.i, a.j)) + } + a.data[i+a.jStride*j] = real(v) + a.data[i+a.jStride*j+1] = imag(v) +} + +func (a twoArray) add(i, j int, v float64) { + if i < 0 || a.i <= i || j < 0 || a.j <= j { + panic(fmt.Sprintf("out of bounds set(%d, %d): bounds i=%d, j=%d", i, j, a.i, a.j)) + } + a.data[i+a.jStride*j] += v +} + +type threeArray struct { + i, j, k int + jStride, kStride int + data []float64 +} + +func newThreeArray(i, j, k int, data []float64) threeArray { + if len(data) < i*j*k { + panic(fmt.Sprintf("short data: len(data)=%d, i=%d, j=%d, k=%d", len(data), i, j, k)) + } + return threeArray{ + i: i, + j: j, + k: k, + jStride: i, + kStride: i * j, + data: data[:i*j*k], + } +} + +func (a threeArray) at(i, j, k int) float64 { + if i < 0 || a.i <= i || j < 0 || a.j <= j || k < 0 || a.k <= k { + panic(fmt.Sprintf("out of bounds at(%d, %d, %d): bounds i=%d, j=%d, k=%d", i, j, k, a.i, a.j, a.k)) + } + return a.data[i+a.jStride*j+a.kStride*k] +} + +func (a threeArray) atCmplx(i, j, k int) complex128 { + if i < 0 || a.i <= i || j < 0 || a.j <= j || k < 0 || a.k <= k { + panic(fmt.Sprintf("out of bounds at(%d, %d, %d): bounds i=%d, j=%d, k=%d", i, j, k, a.i, a.j, a.k)) + } + return complex(a.data[i+a.jStride*j+a.kStride*k], a.data[i+a.jStride*j+a.kStride*k+1]) +} + +func (a threeArray) set(i, j, k int, v float64) { + if i < 0 || a.i <= i || j < 0 || a.j <= j || k < 0 || a.k <= k { + panic(fmt.Sprintf("out of bounds set(%d, %d, %d): bounds i=%d, j=%d, k=%d", i, j, k, a.i, a.j, a.k)) + } + a.data[i+a.jStride*j+a.kStride*k] = v +} + +func (a threeArray) setCmplx(i, j, k int, v complex128) { + if i < 0 || a.i <= i || j < 0 || a.j <= j || k < 0 || a.k <= k { + panic(fmt.Sprintf("out of bounds set(%d, %d, %d): bounds i=%d, j=%d, k=%d", i, j, k, a.i, a.j, a.k)) + } + a.data[i+a.jStride*j+a.kStride*k] = real(v) + a.data[i+a.jStride*j+a.kStride*k+1] = imag(v) +} diff --git a/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/array_no_bounds_checks.go b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/array_no_bounds_checks.go new file mode 100644 index 0000000..52f81b0 --- /dev/null +++ b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/array_no_bounds_checks.go @@ -0,0 +1,83 @@ +// 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 file must be kept in sync with array_bound_checks.go. + +//go:build !bounds +// +build !bounds + +package fftpack + +// The types in array.go implement Fortran-like arrays for bootstrapping +// the implementation of the FFT functions translated from FFTPACK; they +// are column-major. + +type twoArray struct { + jStride int + data []float64 +} + +func newTwoArray(i, j int, data []float64) twoArray { + if len(data) < i*j { + panic("fourier: short data") + } + return twoArray{ + jStride: i, + data: data[:i*j], + } +} + +func (a twoArray) at(i, j int) float64 { + return a.data[i+a.jStride*j] +} + +func (a twoArray) atCmplx(i, j int) complex128 { + return complex(a.data[i+a.jStride*j], a.data[i+a.jStride*j+1]) +} + +func (a twoArray) set(i, j int, v float64) { + a.data[i+a.jStride*j] = v +} + +func (a twoArray) setCmplx(i, j int, v complex128) { + a.data[i+a.jStride*j] = real(v) + a.data[i+a.jStride*j+1] = imag(v) +} + +func (a twoArray) add(i, j int, v float64) { + a.data[i+a.jStride*j] += v +} + +type threeArray struct { + jStride, kStride int + data []float64 +} + +func newThreeArray(i, j, k int, data []float64) threeArray { + if len(data) < i*j*k { + panic("fourier: short data") + } + return threeArray{ + jStride: i, + kStride: i * j, + data: data[:i*j*k], + } +} + +func (a threeArray) at(i, j, k int) float64 { + return a.data[i+a.jStride*j+a.kStride*k] +} + +func (a threeArray) atCmplx(i, j, k int) complex128 { + return complex(a.data[i+a.jStride*j+a.kStride*k], a.data[i+a.jStride*j+a.kStride*k+1]) +} + +func (a threeArray) set(i, j, k int, v float64) { + a.data[i+a.jStride*j+a.kStride*k] = v +} + +func (a threeArray) setCmplx(i, j, k int, v complex128) { + a.data[i+a.jStride*j+a.kStride*k] = real(v) + a.data[i+a.jStride*j+a.kStride*k+1] = imag(v) +} diff --git a/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/cfft.go b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/cfft.go new file mode 100644 index 0000000..099b280 --- /dev/null +++ b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/cfft.go @@ -0,0 +1,588 @@ +// 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 cfft functions by +// Paul N Swarztrauber, placed in the public domain at +// http://www.netlib.org/fftpack/. + +package fftpack + +import ( + "math" + "math/cmplx" +) + +// Cffti initializes the array work which is used in both Cfftf +// and Cfftb. 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. +// +// Output parameters: +// +// work A work array which must be dimensioned at least 4*n. +// the same work array can be used for both Cfftf and Cfftb +// 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 Cfftf or Cfftb. +// +// ifac A work array containing the factors of n. ifac must have +// length 15. +func Cffti(n int, work []float64, ifac []int) { + if len(work) < 4*n { + panic("fourier: short work") + } + if len(ifac) < 15 { + panic("fourier: short ifac") + } + if n == 1 { + return + } + cffti1(n, work[2*n:4*n], ifac[:15]) +} + +func cffti1(n int, wa []float64, ifac []int) { + ntryh := [4]int{3, 4, 2, 5} + + nl := n + nf := 0 + +outer: + for j, ntry := 0, 0; ; j++ { + if j < 4 { + ntry = ntryh[j] + } else { + ntry += 2 + } + for { + if nl%ntry != 0 { + continue outer + } + + ifac[nf+2] = ntry + nl /= ntry + nf++ + + if ntry == 2 && nf != 1 { + for i := 1; i < nf; i++ { + ib := nf - i + 1 + ifac[ib+1] = ifac[ib] + } + ifac[2] = 2 + } + + if nl == 1 { + break outer + } + } + } + + ifac[0] = n + ifac[1] = nf + + argh := 2 * math.Pi / float64(n) + i := 1 + l1 := 1 + for k1 := 0; k1 < nf; k1++ { + ip := ifac[k1+2] + ld := 0 + l2 := l1 * ip + ido := n / l2 + idot := 2*ido + 2 + for j := 0; j < ip-1; j++ { + i1 := i + wa[i-1] = 1 + wa[i] = 0 + ld += l1 + var fi float64 + argld := float64(ld) * argh + for ii := 3; ii < idot; ii += 2 { + i += 2 + fi++ + arg := fi * argld + wa[i-1] = math.Cos(arg) + wa[i] = math.Sin(arg) + } + if ip > 5 { + wa[i1-1] = wa[i-1] + wa[i1] = wa[i] + } + } + l1 = l2 + } +} + +// Cfftf computes the forward complex Discrete Fourier transform +// (the Fourier analysis). Equivalently, Cfftf computes the +// Fourier coefficients of a complex periodic sequence. The +// transform is defined below at output parameter c. +// +// Input parameters: +// +// n The length of the array c to be transformed. The method +// is most efficient when n is a product of small primes. +// n may change so long as different work arrays are provided. +// +// c A complex array of length n which contains the sequence +// to be transformed. +// +// work A real work array which must be dimensioned at least 4*n. +// in the program that calls Cfftf. The work array must be +// initialized by calling subroutine Cffti(n,work,ifac) 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. +// the same work array can be used by Cfftf and Cfftb. +// +// ifac A work array containing the factors of n. ifac must have +// length of at least 15. +// +// Output parameters: +// +// c for j=0, ..., n-1 +// c[j]=the sum from k=0, ..., n-1 of +// c[k]*exp(-i*j*k*2*pi/n) +// +// where i=sqrt(-1) +// +// This transform is unnormalized since a call of Cfftf +// followed by a call of Cfftb will multiply the input +// sequence by n. +// +// The n elements of c are represented in n pairs of real +// values in r where c[j] = r[j*2]+r[j*2+1]i. +// +// work Contains results which must not be destroyed between +// calls of Cfftf or Cfftb. +// ifac Contains results which must not be destroyed between +// calls of Cfftf or Cfftb. +func Cfftf(n int, r, work []float64, ifac []int) { + if len(r) < 2*n { + panic("fourier: short sequence") + } + if len(work) < 4*n { + panic("fourier: short work") + } + if len(ifac) < 15 { + panic("fourier: short ifac") + } + if n == 1 { + return + } + cfft1(n, r[:2*n], work[:2*n], work[2*n:4*n], ifac[:15], -1) +} + +// Cfftb computes the backward complex Discrete Fourier Transform +// (the Fourier synthesis). Equivalently, Cfftf computes the computes +// a complex periodic sequence from its Fourier coefficients. The +// transform is defined below at output parameter c. +// +// Input parameters: +// +// n The length of the array c to be transformed. The method +// is most efficient when n is a product of small primes. +// n may change so long as different work arrays are provided. +// +// c A complex array of length n which contains the sequence +// to be transformed. +// +// work A real work array which must be dimensioned at least 4*n. +// in the program that calls Cfftb. The work array must be +// initialized by calling subroutine Cffti(n,work,ifac) 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. +// The same work array can be used by Cfftf and Cfftb. +// +// ifac A work array containing the factors of n. ifac must have +// length of at least 15. +// +// Output parameters: +// +// c for j=0, ..., n-1 +// c[j]=the sum from k=0, ..., n-1 of +// c[k]*exp(i*j*k*2*pi/n) +// +// where i=sqrt(-1) +// +// This transform is unnormalized since a call of Cfftf +// followed by a call of Cfftb will multiply the input +// sequence by n. +// +// The n elements of c are represented in n pairs of real +// values in r where c[j] = r[j*2]+r[j*2+1]i. +// +// work Contains results which must not be destroyed between +// calls of Cfftf or Cfftb. +// ifac Contains results which must not be destroyed between +// calls of Cfftf or Cfftb. +func Cfftb(n int, c, work []float64, ifac []int) { + if len(c) < 2*n { + panic("fourier: short sequence") + } + if len(work) < 4*n { + panic("fourier: short work") + } + if len(ifac) < 15 { + panic("fourier: short ifac") + } + if n == 1 { + return + } + cfft1(n, c[:2*n], work[:2*n], work[2*n:4*n], ifac[:15], 1) +} + +// cfft1 implements cfftf1 and cfftb1 depending on sign. +func cfft1(n int, c, ch, wa []float64, ifac []int, sign float64) { + nf := ifac[1] + na := false + l1 := 1 + iw := 0 + + for k1 := 1; k1 <= nf; k1++ { + ip := ifac[k1+1] + l2 := ip * l1 + ido := n / l2 + idot := 2 * ido + idl1 := idot * l1 + + switch ip { + case 4: + ix2 := iw + idot + ix3 := ix2 + idot + if na { + pass4(idot, l1, ch, c, wa[iw:], wa[ix2:], wa[ix3:], sign) + } else { + pass4(idot, l1, c, ch, wa[iw:], wa[ix2:], wa[ix3:], sign) + } + na = !na + case 2: + if na { + pass2(idot, l1, ch, c, wa[iw:], sign) + } else { + pass2(idot, l1, c, ch, wa[iw:], sign) + } + na = !na + case 3: + ix2 := iw + idot + if na { + pass3(idot, l1, ch, c, wa[iw:], wa[ix2:], sign) + } else { + pass3(idot, l1, c, ch, wa[iw:], wa[ix2:], sign) + } + na = !na + case 5: + ix2 := iw + idot + ix3 := ix2 + idot + ix4 := ix3 + idot + if na { + pass5(idot, l1, ch, c, wa[iw:], wa[ix2:], wa[ix3:], wa[ix4:], sign) + } else { + pass5(idot, l1, c, ch, wa[iw:], wa[ix2:], wa[ix3:], wa[ix4:], sign) + } + na = !na + default: + var nac bool + if na { + nac = pass(idot, ip, l1, idl1, ch, ch, ch, c, c, wa[iw:], sign) + } else { + nac = pass(idot, ip, l1, idl1, c, c, c, ch, ch, wa[iw:], sign) + } + if nac { + na = !na + } + } + + l1 = l2 + iw += (ip - 1) * idot + } + + if na { + for i := 0; i < 2*n; i++ { + c[i] = ch[i] + } + } +} + +// pass2 implements passf2 and passb2 depending on sign. +func pass2(ido, l1 int, cc, ch, wa1 []float64, sign float64) { + cc3 := newThreeArray(ido, 2, l1, cc) + ch3 := newThreeArray(ido, l1, 2, ch) + + if ido <= 2 { + for k := 0; k < l1; k++ { + ch3.setCmplx(0, k, 0, cc3.atCmplx(0, 0, k)+cc3.atCmplx(0, 1, k)) + ch3.setCmplx(0, k, 1, cc3.atCmplx(0, 0, k)-cc3.atCmplx(0, 1, k)) + } + return + } + for k := 0; k < l1; k++ { + for i := 1; i < ido; i += 2 { + ch3.setCmplx(i-1, k, 0, cc3.atCmplx(i-1, 0, k)+cc3.atCmplx(i-1, 1, k)) + t2 := cc3.atCmplx(i-1, 0, k) - cc3.atCmplx(i-1, 1, k) + ch3.setCmplx(i-1, k, 1, complex(wa1[i-1], sign*wa1[i])*t2) + } + } +} + +// pass3 implements passf3 and passb3 depending on sign. +func pass3(ido, l1 int, cc, ch, wa1, wa2 []float64, sign float64) { + const ( + taur = -0.5 + taui = 0.866025403784439 // sqrt(3)/2 + ) + + cc3 := newThreeArray(ido, 3, l1, cc) + ch3 := newThreeArray(ido, l1, 3, ch) + + if ido == 2 { + for k := 0; k < l1; k++ { + t2 := cc3.atCmplx(0, 1, k) + cc3.atCmplx(0, 2, k) + ch3.setCmplx(0, k, 0, cc3.atCmplx(0, 0, k)+t2) + + c2 := cc3.atCmplx(0, 0, k) + scale(taur, t2) + c3 := cmplx.Conj(swap(scale(sign*taui, cc3.atCmplx(0, 1, k)-cc3.atCmplx(0, 2, k)))) + ch3.setCmplx(0, k, 1, c2-c3) + ch3.setCmplx(0, k, 2, c2+c3) + } + return + } + for k := 0; k < l1; k++ { + for i := 1; i < ido; i += 2 { + t2 := cc3.atCmplx(i-1, 1, k) + cc3.atCmplx(i-1, 2, k) + ch3.setCmplx(i-1, k, 0, cc3.atCmplx(i-1, 0, k)+t2) + + c2 := cc3.atCmplx(i-1, 0, k) + scale(taur, t2) + c3 := cmplx.Conj(swap(scale(sign*taui, cc3.atCmplx(i-1, 1, k)-cc3.atCmplx(i-1, 2, k)))) + d2 := c2 - c3 + d3 := c2 + c3 + ch3.setCmplx(i-1, k, 1, complex(wa1[i-1], sign*wa1[i])*d2) + ch3.setCmplx(i-1, k, 2, complex(wa2[i-1], sign*wa2[i])*d3) + } + } +} + +// pass4 implements passf4 and passb4 depending on sign. +func pass4(ido, l1 int, cc, ch, wa1, wa2, wa3 []float64, sign float64) { + cc3 := newThreeArray(ido, 4, l1, cc) + ch3 := newThreeArray(ido, l1, 4, ch) + + if ido == 2 { + for k := 0; k < l1; k++ { + t1 := cc3.atCmplx(0, 0, k) - cc3.atCmplx(0, 2, k) + t2 := cc3.atCmplx(0, 0, k) + cc3.atCmplx(0, 2, k) + t3 := cc3.atCmplx(0, 1, k) + cc3.atCmplx(0, 3, k) + t4 := cmplx.Conj(swap(scale(sign, cc3.atCmplx(0, 3, k)-cc3.atCmplx(0, 1, k)))) + + ch3.setCmplx(0, k, 0, t2+t3) + ch3.setCmplx(0, k, 1, t1+t4) + ch3.setCmplx(0, k, 2, t2-t3) + ch3.setCmplx(0, k, 3, t1-t4) + } + return + } + for k := 0; k < l1; k++ { + for i := 1; i < ido; i += 2 { + t1 := cc3.atCmplx(i-1, 0, k) - cc3.atCmplx(i-1, 2, k) + t2 := cc3.atCmplx(i-1, 0, k) + cc3.atCmplx(i-1, 2, k) + t3 := cc3.atCmplx(i-1, 1, k) + cc3.atCmplx(i-1, 3, k) + t4 := cmplx.Conj(swap(scale(sign, cc3.atCmplx(i-1, 3, k)-cc3.atCmplx(i-1, 1, k)))) + ch3.setCmplx(i-1, k, 0, t2+t3) + + c2 := t1 + t4 + c3 := t2 - t3 + c4 := t1 - t4 + ch3.setCmplx(i-1, k, 1, complex(wa1[i-1], sign*wa1[i])*c2) + ch3.setCmplx(i-1, k, 2, complex(wa2[i-1], sign*wa2[i])*c3) + ch3.setCmplx(i-1, k, 3, complex(wa3[i-1], sign*wa3[i])*c4) + } + } +} + +// pass5 implements passf5 and passb5 depending on sign. +func pass5(ido, l1 int, cc, ch, wa1, wa2, wa3, wa4 []float64, sign float64) { + const ( + tr11 = 0.309016994374947 + ti11 = 0.951056516295154 + tr12 = -0.809016994374947 + ti12 = 0.587785252292473 + ) + + cc3 := newThreeArray(ido, 5, l1, cc) + ch3 := newThreeArray(ido, l1, 5, ch) + + if ido == 2 { + for k := 0; k < l1; k++ { + t2 := cc3.atCmplx(0, 1, k) + cc3.atCmplx(0, 4, k) + t3 := cc3.atCmplx(0, 2, k) + cc3.atCmplx(0, 3, k) + t4 := cc3.atCmplx(0, 2, k) - cc3.atCmplx(0, 3, k) + t5 := cc3.atCmplx(0, 1, k) - cc3.atCmplx(0, 4, k) + ch3.setCmplx(0, k, 0, cc3.atCmplx(0, 0, k)+t2+t3) + + c2 := cc3.atCmplx(0, 0, k) + scale(tr11, t2) + scale(tr12, t3) + c3 := cc3.atCmplx(0, 0, k) + scale(tr12, t2) + scale(tr11, t3) + c4 := cmplx.Conj(swap(scale(sign, scale(ti12, t5)-scale(ti11, t4)))) + c5 := cmplx.Conj(swap(scale(sign, scale(ti11, t5)+scale(ti12, t4)))) + ch3.setCmplx(0, k, 1, c2-c5) + ch3.setCmplx(0, k, 2, c3-c4) + ch3.setCmplx(0, k, 3, c3+c4) + ch3.setCmplx(0, k, 4, c2+c5) + } + return + } + for k := 0; k < l1; k++ { + for i := 1; i < ido; i += 2 { + t2 := cc3.atCmplx(i-1, 1, k) + cc3.atCmplx(i-1, 4, k) + t3 := cc3.atCmplx(i-1, 2, k) + cc3.atCmplx(i-1, 3, k) + t4 := cc3.atCmplx(i-1, 2, k) - cc3.atCmplx(i-1, 3, k) + t5 := cc3.atCmplx(i-1, 1, k) - cc3.atCmplx(i-1, 4, k) + ch3.setCmplx(i-1, k, 0, cc3.atCmplx(i-1, 0, k)+t2+t3) + + c2 := cc3.atCmplx(i-1, 0, k) + scale(tr11, t2) + scale(tr12, t3) + c3 := cc3.atCmplx(i-1, 0, k) + scale(tr12, t2) + scale(tr11, t3) + c4 := cmplx.Conj(swap(scale(sign, scale(ti12, t5)-scale(ti11, t4)))) + c5 := cmplx.Conj(swap(scale(sign, scale(ti11, t5)+scale(ti12, t4)))) + d2 := c2 - c5 + d3 := c3 - c4 + d4 := c3 + c4 + d5 := c2 + c5 + ch3.setCmplx(i-1, k, 1, complex(wa1[i-1], sign*wa1[i])*d2) + ch3.setCmplx(i-1, k, 2, complex(wa2[i-1], sign*wa2[i])*d3) + ch3.setCmplx(i-1, k, 3, complex(wa3[i-1], sign*wa3[i])*d4) + ch3.setCmplx(i-1, k, 4, complex(wa4[i-1], sign*wa4[i])*d5) + } + } +} + +// pass implements passf and passb depending on sign. +func pass(ido, ip, l1, idl1 int, cc, c1, c2, ch, ch2, wa []float64, sign float64) (nac bool) { + cc3 := newThreeArray(ido, ip, l1, cc) + c13 := newThreeArray(ido, l1, ip, c1) + ch3 := newThreeArray(ido, l1, ip, ch) + c2m := newTwoArray(idl1, ip, c2) + ch2m := newTwoArray(idl1, ip, ch2) + + idot := ido / 2 + ipph := (ip + 1) / 2 + idp := ip * ido + + if ido < l1 { + for j := 1; j < ipph; j++ { + jc := ip - j + for i := 0; i < ido; i++ { + for k := 0; k < l1; k++ { + ch3.set(i, k, j, cc3.at(i, j, k)+cc3.at(i, jc, k)) + ch3.set(i, k, jc, cc3.at(i, j, k)-cc3.at(i, jc, k)) + } + } + } + for i := 0; i < ido; i++ { + for k := 0; k < l1; k++ { + ch3.set(i, k, 0, cc3.at(i, 0, k)) + } + } + } else { + for j := 1; j < ipph; j++ { + jc := ip - j + for k := 0; k < l1; k++ { + for i := 0; i < ido; i++ { + ch3.set(i, k, j, cc3.at(i, j, k)+cc3.at(i, jc, k)) + ch3.set(i, k, jc, cc3.at(i, j, k)-cc3.at(i, jc, k)) + } + } + } + for k := 0; k < l1; k++ { + for i := 0; i < ido; i++ { + ch3.set(i, k, 0, cc3.at(i, 0, k)) + } + } + } + + idl := 1 - ido + inc := 0 + for l := 1; l < ipph; l++ { + lc := ip - l + idl += ido + for ik := 0; ik < idl1; ik++ { + c2m.set(ik, l, ch2m.at(ik, 0)+wa[idl-1]*ch2m.at(ik, 1)) + c2m.set(ik, lc, sign*wa[idl]*ch2m.at(ik, ip-1)) + } + idlj := idl + inc += ido + for j := 2; j < ipph; j++ { + jc := ip - j + idlj += inc + if idlj > idp { + idlj -= idp + } + war := wa[idlj-1] + wai := wa[idlj] + for ik := 0; ik < idl1; ik++ { + c2m.add(ik, l, war*ch2m.at(ik, j)) + c2m.add(ik, lc, sign*wai*ch2m.at(ik, jc)) + } + } + } + + for j := 1; j < ipph; j++ { + for ik := 0; ik < idl1; ik++ { + ch2m.add(ik, 0, ch2m.at(ik, j)) + } + } + + for j := 1; j < ipph; j++ { + jc := ip - j + for ik := 1; ik < idl1; ik += 2 { + ch2m.setCmplx(ik-1, j, c2m.atCmplx(ik-1, j)-cmplx.Conj(swap(c2m.atCmplx(ik-1, jc)))) + ch2m.setCmplx(ik-1, jc, c2m.atCmplx(ik-1, j)+cmplx.Conj(swap(c2m.atCmplx(ik-1, jc)))) + } + } + + if ido == 2 { + return true + } + + for ik := 0; ik < idl1; ik++ { + c2m.set(ik, 0, ch2m.at(ik, 0)) + } + + for j := 1; j < ip; j++ { + for k := 0; k < l1; k++ { + c13.setCmplx(0, k, j, ch3.atCmplx(0, k, j)) + } + } + + if idot > l1 { + idj := 1 - ido + for j := 1; j < ip; j++ { + idj += ido + for k := 0; k < l1; k++ { + idij := idj + for i := 3; i < ido; i += 2 { + idij += 2 + c13.setCmplx(i-1, k, j, complex(wa[idij-1], sign*wa[idij])*ch3.atCmplx(i-1, k, j)) + } + } + } + + return false + } + + idij := -1 + for j := 1; j < ip; j++ { + idij += 2 + for i := 3; i < ido; i += 2 { + idij += 2 + for k := 0; k < l1; k++ { + c13.setCmplx(i-1, k, j, complex(wa[idij-1], sign*wa[idij])*ch3.atCmplx(i-1, k, j)) + } + } + } + return false +} diff --git a/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/cosq.go b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/cosq.go new file mode 100644 index 0000000..58ef516 --- /dev/null +++ b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/cosq.go @@ -0,0 +1,219 @@ +// 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 +} diff --git a/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/cost.go b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/cost.go new file mode 100644 index 0000000..46e808a --- /dev/null +++ b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/cost.go @@ -0,0 +1,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 + } + } +} diff --git a/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/doc.go b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/doc.go new file mode 100644 index 0000000..74c2965 --- /dev/null +++ b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/doc.go @@ -0,0 +1,7 @@ +// 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. + +// Package fftpack implements Discrete Fourier Transform functions +// ported from the Fortran implementation of FFTPACK. +package fftpack // import "gonum.org/v1/gonum/dsp/fourier/internal/fftpack" diff --git a/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/helpers.go b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/helpers.go new file mode 100644 index 0000000..0278414 --- /dev/null +++ b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/helpers.go @@ -0,0 +1,15 @@ +// Copyright ©2020 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. + +package fftpack + +// swap returns c with the real and imaginary parts swapped. +func swap(c complex128) complex128 { + return complex(imag(c), real(c)) +} + +// scale scales the complex number c by f. +func scale(f float64, c complex128) complex128 { + return complex(f*real(c), f*imag(c)) +} diff --git a/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/rfft.go b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/rfft.go new file mode 100644 index 0000000..1472cca --- /dev/null +++ b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/rfft.go @@ -0,0 +1,1079 @@ +// 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 rfft functions by +// Paul N Swarztrauber, placed in the public domain at +// http://www.netlib.org/fftpack/. + +package fftpack + +import ( + "math" + "math/cmplx" +) + +// Rffti initializes the array work which is used in both Rfftf +// and Rfftb. 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. +// +// Output parameters: +// +// work A work array which must be dimensioned at least 2*n. +// The same work array can be used for both Rfftf and Rfftb +// 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 Rfftf or Rfftb. +// +// ifac A work array containing the factors of n. ifac must have +// length of at least 15. +func Rffti(n int, work []float64, ifac []int) { + if len(work) < 2*n { + panic("fourier: short work") + } + if len(ifac) < 15 { + panic("fourier: short ifac") + } + if n == 1 { + return + } + rffti1(n, work[n:2*n], ifac[:15]) +} + +func rffti1(n int, wa []float64, ifac []int) { + ntryh := [4]int{4, 2, 3, 5} + + nl := n + nf := 0 + +outer: + for j, ntry := 0, 0; ; j++ { + if j < 4 { + ntry = ntryh[j] + } else { + ntry += 2 + } + for { + if nl%ntry != 0 { + continue outer + } + + ifac[nf+2] = ntry + nl /= ntry + nf++ + + if ntry == 2 && nf != 1 { + for i := 1; i < nf; i++ { + ib := nf - i + 1 + ifac[ib+1] = ifac[ib] + } + ifac[2] = 2 + } + + if nl == 1 { + break outer + } + } + } + + ifac[0] = n + ifac[1] = nf + if nf == 1 { + return + } + argh := 2 * math.Pi / float64(n) + + is := 0 + l1 := 1 + for k1 := 0; k1 < nf-1; k1++ { + ip := ifac[k1+2] + ld := 0 + l2 := l1 * ip + ido := n / l2 + for j := 0; j < ip-1; j++ { + ld += l1 + i := is + fi := 0.0 + argld := float64(ld) * argh + for ii := 2; ii < ido; ii += 2 { + fi++ + arg := fi * argld + wa[i] = math.Cos(arg) + wa[i+1] = math.Sin(arg) + i += 2 + } + is += ido + } + l1 = l2 + } +} + +// Rfftf computes the Fourier coefficients of a real perodic sequence +// (Fourier analysis). The transform is defined below at output +// parameter r. +// +// Input parameters: +// +// n The length of the array r to be transformed. The method +// is most efficient when n is a product of small primes. +// n may change so long as different work arrays are provided. +// +// r A real array of length n which contains the sequence +// to be transformed. +// +// work a work array which must be dimensioned at least 2*n. +// in the program that calls Rfftf. the work array must be +// initialized by calling subroutine rffti(n,work,ifac) 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. +// The same work array can be used by Rfftf and Rfftb. +// +// ifac A work array containing the factors of n. ifac must have +// length of at least 15. +// +// Output parameters: +// +// r r[0] = the sum from i=0 to i=n-1 of r[i] +// +// if n is even set l=n/2, if n is odd set l = (n+1)/2 +// then for k = 1, ..., l-1 +// r[2*k-1] = the sum from i = 0 to i = n-1 of +// r[i]*cos(k*i*2*pi/n) +// r[2*k] = the sum from i = 0 to i = n-1 of +// -r[i]*sin(k*i*2*pi/n) +// +// if n is even +// r[n-1] = the sum from i = 0 to i = n-1 of +// (-1)^i*r[i] +// +// This transform is unnormalized since a call of Rfftf +// followed by a call of Rfftb will multiply the input +// sequence by n. +// +// work contains results which must not be destroyed between +// calls of Rfftf or Rfftb. +// ifac contains results which must not be destroyed between +// calls of Rfftf or Rfftb. +func Rfftf(n int, r, work []float64, ifac []int) { + if len(r) < n { + panic("fourier: short sequence") + } + if len(work) < 2*n { + panic("fourier: short work") + } + if len(ifac) < 15 { + panic("fourier: short ifac") + } + if n == 1 { + return + } + rfftf1(n, r[:n], work[:n], work[n:2*n], ifac[:15]) +} + +func rfftf1(n int, c, ch, wa []float64, ifac []int) { + nf := ifac[1] + na := true + l2 := n + iw := n - 1 + + for k1 := 1; k1 <= nf; k1++ { + kh := nf - k1 + ip := ifac[kh+2] + l1 := l2 / ip + ido := n / l2 + idl1 := ido * l1 + iw -= (ip - 1) * ido + na = !na + + switch ip { + case 4: + ix2 := iw + ido + ix3 := ix2 + ido + if na { + radf4(ido, l1, ch, c, wa[iw:], wa[ix2:], wa[ix3:]) + } else { + radf4(ido, l1, c, ch, wa[iw:], wa[ix2:], wa[ix3:]) + } + case 2: + if na { + radf2(ido, l1, ch, c, wa[iw:]) + } else { + radf2(ido, l1, c, ch, wa[iw:]) + } + case 3: + ix2 := iw + ido + if na { + radf3(ido, l1, ch, c, wa[iw:], wa[ix2:]) + } else { + radf3(ido, l1, c, ch, wa[iw:], wa[ix2:]) + } + case 5: + ix2 := iw + ido + ix3 := ix2 + ido + ix4 := ix3 + ido + if na { + radf5(ido, l1, ch, c, wa[iw:], wa[ix2:], wa[ix3:], wa[ix4:]) + } else { + radf5(ido, l1, c, ch, wa[iw:], wa[ix2:], wa[ix3:], wa[ix4:]) + } + default: + if ido == 1 { + na = !na + } + if na { + radfg(ido, ip, l1, idl1, ch, ch, ch, c, c, wa[iw:]) + na = false + } else { + radfg(ido, ip, l1, idl1, c, c, c, ch, ch, wa[iw:]) + na = true + } + } + + l2 = l1 + } + + if na { + return + } + for i := 0; i < n; i++ { + c[i] = ch[i] + } +} + +func radf2(ido, l1 int, cc, ch, wa1 []float64) { + cc3 := newThreeArray(ido, l1, 2, cc) + ch3 := newThreeArray(ido, 2, l1, ch) + + for k := 0; k < l1; k++ { + ch3.set(0, 0, k, cc3.at(0, k, 0)+cc3.at(0, k, 1)) + ch3.set(ido-1, 1, k, cc3.at(0, k, 0)-cc3.at(0, k, 1)) + } + if ido < 2 { + return + } + if ido > 2 { + idp2 := ido + 1 + for k := 0; k < l1; k++ { + for i := 2; i < ido; i += 2 { + ic := idp2 - (i + 1) + + t2 := complex(wa1[i-2], -wa1[i-1]) * cc3.atCmplx(i-1, k, 1) + ch3.setCmplx(i-1, 0, k, cc3.atCmplx(i-1, k, 0)+t2) + + // This is left as conj(z1)-conj(z2) rather than conj(z1-z2) + // to retain current signed zero behaviour. + ch3.setCmplx(ic-1, 1, k, cmplx.Conj(cc3.atCmplx(i-1, k, 0))-cmplx.Conj(t2)) + } + } + if ido%2 == 1 { + return + } + } + for k := 0; k < l1; k++ { + ch3.set(0, 1, k, -cc3.at(ido-1, k, 1)) + ch3.set(ido-1, 0, k, cc3.at(ido-1, k, 0)) + } +} + +func radf3(ido, l1 int, cc, ch, wa1, wa2 []float64) { + const ( + taur = -0.5 + taui = 0.866025403784439 // sqrt(3)/2 + ) + + cc3 := newThreeArray(ido, l1, 3, cc) + ch3 := newThreeArray(ido, 3, l1, ch) + + for k := 0; k < l1; k++ { + cr2 := cc3.at(0, k, 1) + cc3.at(0, k, 2) + ch3.set(0, 0, k, cc3.at(0, k, 0)+cr2) + ch3.set(0, 2, k, taui*(cc3.at(0, k, 2)-cc3.at(0, k, 1))) + ch3.set(ido-1, 1, k, cc3.at(0, k, 0)+taur*cr2) + } + if ido < 2 { + return + } + idp2 := ido + 1 + for k := 0; k < l1; k++ { + for i := 2; i < ido; i += 2 { + ic := idp2 - (i + 1) + + d2 := complex(wa1[i-2], -wa1[i-1]) * cc3.atCmplx(i-1, k, 1) + d3 := complex(wa2[i-2], -wa2[i-1]) * cc3.atCmplx(i-1, k, 2) + c2 := d2 + d3 + ch3.setCmplx(i-1, 0, k, cc3.atCmplx(i-1, k, 0)+c2) + + t2 := cc3.atCmplx(i-1, k, 0) + scale(taur, c2) + t3 := scale(taui, cmplx.Conj(swap(d2-d3))) + ch3.setCmplx(i-1, 2, k, t2+t3) + ch3.setCmplx(ic-1, 1, k, cmplx.Conj(t2-t3)) + } + } +} + +func radf4(ido, l1 int, cc, ch, wa1, wa2, wa3 []float64) { + const hsqt2 = math.Sqrt2 / 2 + + cc3 := newThreeArray(ido, l1, 4, cc) + ch3 := newThreeArray(ido, 4, l1, ch) + + for k := 0; k < l1; k++ { + tr1 := cc3.at(0, k, 1) + cc3.at(0, k, 3) + tr2 := cc3.at(0, k, 0) + cc3.at(0, k, 2) + ch3.set(0, 0, k, tr1+tr2) + ch3.set(ido-1, 3, k, tr2-tr1) + ch3.set(ido-1, 1, k, cc3.at(0, k, 0)-cc3.at(0, k, 2)) + ch3.set(0, 2, k, cc3.at(0, k, 3)-cc3.at(0, k, 1)) + } + if ido < 2 { + return + } + if ido > 2 { + idp2 := ido + 1 + for k := 0; k < l1; k++ { + for i := 2; i < ido; i += 2 { + ic := idp2 - (i + 1) + + c2 := complex(wa1[i-2], -wa1[i-1]) * cc3.atCmplx(i-1, k, 1) + c3 := complex(wa2[i-2], -wa2[i-1]) * cc3.atCmplx(i-1, k, 2) + c4 := complex(wa3[i-2], -wa3[i-1]) * cc3.atCmplx(i-1, k, 3) + t1 := c2 + c4 + t2 := cc3.atCmplx(i-1, k, 0) + c3 + t3 := cc3.atCmplx(i-1, k, 0) - c3 + t4 := cmplx.Conj(c4 - c2) + ch3.setCmplx(i-1, 0, k, t1+t2) + ch3.setCmplx(ic-1, 3, k, cmplx.Conj(t2-t1)) + ch3.setCmplx(i-1, 2, k, swap(t4)+t3) + ch3.setCmplx(ic-1, 1, k, cmplx.Conj(t3-swap(t4))) + } + } + + if ido%2 == 1 { + return + } + } + for k := 0; k < l1; k++ { + ti1 := -hsqt2 * (cc3.at(ido-1, k, 1) + cc3.at(ido-1, k, 3)) + tr1 := hsqt2 * (cc3.at(ido-1, k, 1) - cc3.at(ido-1, k, 3)) + ch3.set(ido-1, 0, k, tr1+cc3.at(ido-1, k, 0)) + ch3.set(ido-1, 2, k, cc3.at(ido-1, k, 0)-tr1) + ch3.set(0, 1, k, ti1-cc3.at(ido-1, k, 2)) + ch3.set(0, 3, k, ti1+cc3.at(ido-1, k, 2)) + } +} + +func radf5(ido, l1 int, cc, ch, wa1, wa2, wa3, wa4 []float64) { + const ( + tr11 = 0.309016994374947 + ti11 = 0.951056516295154 + tr12 = -0.809016994374947 + ti12 = 0.587785252292473 + ) + + cc3 := newThreeArray(ido, l1, 5, cc) + ch3 := newThreeArray(ido, 5, l1, ch) + + for k := 0; k < l1; k++ { + cr2 := cc3.at(0, k, 4) + cc3.at(0, k, 1) + cr3 := cc3.at(0, k, 3) + cc3.at(0, k, 2) + ci4 := cc3.at(0, k, 3) - cc3.at(0, k, 2) + ci5 := cc3.at(0, k, 4) - cc3.at(0, k, 1) + ch3.set(0, 0, k, cc3.at(0, k, 0)+cr2+cr3) + ch3.set(ido-1, 1, k, cc3.at(0, k, 0)+tr11*cr2+tr12*cr3) + ch3.set(0, 2, k, ti11*ci5+ti12*ci4) + ch3.set(ido-1, 3, k, cc3.at(0, k, 0)+tr12*cr2+tr11*cr3) + ch3.set(0, 4, k, ti12*ci5-ti11*ci4) + } + + if ido < 2 { + return + } + idp2 := ido + 1 + for k := 0; k < l1; k++ { + for i := 2; i < ido; i += 2 { + ic := idp2 - (i + 1) + + d2 := complex(wa1[i-2], -wa1[i-1]) * cc3.atCmplx(i-1, k, 1) + d3 := complex(wa2[i-2], -wa2[i-1]) * cc3.atCmplx(i-1, k, 2) + d4 := complex(wa3[i-2], -wa3[i-1]) * cc3.atCmplx(i-1, k, 3) + d5 := complex(wa4[i-2], -wa4[i-1]) * cc3.atCmplx(i-1, k, 4) + c2 := d2 + d5 + c3 := d3 + d4 + c4 := cmplx.Conj(swap(d3 - d4)) + c5 := cmplx.Conj(swap(d2 - d5)) + ch3.setCmplx(i-1, 0, k, cc3.atCmplx(i-1, k, 0)+c2+c3) + + t2 := cc3.atCmplx(i-1, k, 0) + scale(tr11, c2) + scale(tr12, c3) + t3 := cc3.atCmplx(i-1, k, 0) + scale(tr12, c2) + scale(tr11, c3) + t4 := scale(ti12, c5) - scale(ti11, c4) + t5 := scale(ti11, c5) + scale(ti12, c4) + ch3.setCmplx(ic-1, 1, k, cmplx.Conj(t2-t5)) + ch3.setCmplx(i-1, 2, k, t2+t5) + ch3.setCmplx(ic-1, 3, k, cmplx.Conj(t3-t4)) + ch3.setCmplx(i-1, 4, k, t3+t4) + } + } +} + +func radfg(ido, ip, l1, idl1 int, cc, c1, c2, ch, ch2, wa []float64) { + cc3 := newThreeArray(ido, ip, l1, cc) + c13 := newThreeArray(ido, l1, ip, c1) + ch3 := newThreeArray(ido, l1, ip, ch) + c2m := newTwoArray(idl1, ip, c2) + ch2m := newTwoArray(idl1, ip, ch2) + + arg := 2 * math.Pi / float64(ip) + dcp := math.Cos(arg) + dsp := math.Sin(arg) + ipph := (ip + 1) / 2 + nbd := (ido - 1) / 2 + + if ido == 1 { + for ik := 0; ik < idl1; ik++ { + c2m.set(ik, 0, ch2m.at(ik, 0)) + } + } else { + for ik := 0; ik < idl1; ik++ { + ch2m.set(ik, 0, c2m.at(ik, 0)) + } + for j := 1; j < ip; j++ { + for k := 0; k < l1; k++ { + ch3.set(0, k, j, c13.at(0, k, j)) + } + } + + is := -ido - 1 + if nbd > l1 { + for j := 1; j < ip; j++ { + is += ido + for k := 0; k < l1; k++ { + idij := is + for i := 2; i < ido; i += 2 { + idij += 2 + ch3.setCmplx(i-1, k, j, complex(wa[idij-1], -wa[idij])*c13.atCmplx(i-1, k, j)) + } + } + } + } else { + for j := 1; j < ip; j++ { + is += ido + idij := is + for i := 2; i < ido; i += 2 { + idij += 2 + for k := 0; k < l1; k++ { + ch3.setCmplx(i-1, k, j, complex(wa[idij-1], -wa[idij])*c13.atCmplx(i-1, k, j)) + } + } + } + } + if nbd < l1 { + for j := 1; j < ipph; j++ { + jc := ip - j + for i := 2; i < ido; i += 2 { + for k := 0; k < l1; k++ { + c13.setCmplx(i-1, k, j, ch3.atCmplx(i-1, k, j)+ch3.atCmplx(i-1, k, jc)) + c13.setCmplx(i-1, k, jc, cmplx.Conj(swap(ch3.atCmplx(i-1, k, j)-ch3.atCmplx(i-1, k, jc)))) + } + } + } + } else { + for j := 1; j < ipph; j++ { + jc := ip - j + for k := 0; k < l1; k++ { + for i := 2; i < ido; i += 2 { + c13.setCmplx(i-1, k, j, ch3.atCmplx(i-1, k, j)+ch3.atCmplx(i-1, k, jc)) + c13.setCmplx(i-1, k, jc, cmplx.Conj(swap(ch3.atCmplx(i-1, k, j)-ch3.atCmplx(i-1, k, jc)))) + } + } + } + } + } + + for j := 1; j < ipph; j++ { + jc := ip - j + for k := 0; k < l1; k++ { + c13.set(0, k, j, ch3.at(0, k, j)+ch3.at(0, k, jc)) + c13.set(0, k, jc, ch3.at(0, k, jc)-ch3.at(0, k, j)) + } + } + ar1 := 1.0 + ai1 := 0.0 + for l := 1; l < ipph; l++ { + lc := ip - l + ar1h := dcp*ar1 - dsp*ai1 + ai1 = dcp*ai1 + dsp*ar1 + ar1 = ar1h + for ik := 0; ik < idl1; ik++ { + ch2m.set(ik, l, c2m.at(ik, 0)+ar1*c2m.at(ik, 1)) + ch2m.set(ik, lc, ai1*c2m.at(ik, ip-1)) + } + dc2 := ar1 + ds2 := ai1 + ar2 := ar1 + ai2 := ai1 + for j := 2; j < ipph; j++ { + jc := ip - j + ar2h := dc2*ar2 - ds2*ai2 + ai2 = dc2*ai2 + ds2*ar2 + ar2 = ar2h + for ik := 0; ik < idl1; ik++ { + ch2m.add(ik, l, ar2*c2m.at(ik, j)) + ch2m.add(ik, lc, ai2*c2m.at(ik, jc)) + } + } + } + for j := 1; j < ipph; j++ { + for ik := 0; ik < idl1; ik++ { + ch2m.add(ik, 0, c2m.at(ik, j)) + } + } + + if ido < l1 { + for i := 0; i < ido; i++ { + for k := 0; k < l1; k++ { + cc3.set(i, 0, k, ch3.at(i, k, 0)) + } + } + } else { + for k := 0; k < l1; k++ { + for i := 0; i < ido; i++ { + cc3.set(i, 0, k, ch3.at(i, k, 0)) + } + } + } + for j := 1; j < ipph; j++ { + jc := ip - j + j2 := 2 * j + for k := 0; k < l1; k++ { + cc3.set(ido-1, j2-1, k, ch3.at(0, k, j)) + cc3.set(0, j2, k, ch3.at(0, k, jc)) + } + } + + if ido == 1 { + return + } + if nbd < l1 { + for j := 1; j < ipph; j++ { + jc := ip - j + j2 := 2 * j + for i := 2; i < ido; i += 2 { + ic := ido - i + for k := 0; k < l1; k++ { + cc3.setCmplx(i-1, j2, k, ch3.atCmplx(i-1, k, j)+ch3.atCmplx(i-1, k, jc)) + cc3.setCmplx(ic-1, j2-1, k, cmplx.Conj(ch3.atCmplx(i-1, k, j)-ch3.atCmplx(i-1, k, jc))) + } + } + } + return + } + for j := 1; j < ipph; j++ { + jc := ip - j + j2 := 2 * j + for k := 0; k < l1; k++ { + for i := 2; i < ido; i += 2 { + ic := ido - i + + cc3.setCmplx(i-1, j2, k, ch3.atCmplx(i-1, k, j)+ch3.atCmplx(i-1, k, jc)) + cc3.setCmplx(ic-1, j2-1, k, cmplx.Conj(ch3.atCmplx(i-1, k, j)-ch3.atCmplx(i-1, k, jc))) + } + } + } +} + +// Rfftb computes the real perodic sequence from its Fourier +// coefficients (Fourier synthesis). The transform is defined +// below at output parameter r. +// +// Input parameters +// +// n The length of the array r to be transformed. The method +// is most efficient when n is a product of small primes. +// n may change so long as different work arrays are provided. +// +// r A real array of length n which contains the sequence +// to be transformed. +// +// work A work array which must be dimensioned at least 2*n. +// in the program that calls Rfftb. The work array must be +// initialized by calling subroutine rffti(n,work,ifac) 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. +// The same work array can be used by Rfftf and Rfftb. +// +// ifac A work array containing the factors of n. ifac must have +// length of at least 15. +// +// output parameters +// +// r for n even and for i = 0, ..., n +// r[i] = r[0]+(-1)^i*r[n-1] +// plus the sum from k=1 to k=n/2-1 of +// 2*r(2*k-1)*cos(k*i*2*pi/n) +// -2*r(2*k)*sin(k*i*2*pi/n) +// +// for n odd and for i = 0, ..., n-1 +// r[i] = r[0] plus the sum from k=1 to k=(n-1)/2 of +// 2*r(2*k-1)*cos(k*i*2*pi/n) +// -2*r(2*k)*sin(k*i*2*pi/n) +// +// This transform is unnormalized since a call of Rfftf +// followed by a call of Rfftb will multiply the input +// sequence by n. +// +// work Contains results which must not be destroyed between +// calls of Rfftf or Rfftb. +// ifac Contains results which must not be destroyed between +// calls of Rfftf or Rfftb. +func Rfftb(n int, r, work []float64, ifac []int) { + if len(r) < n { + panic("fourier: short sequence") + } + if len(work) < 2*n { + panic("fourier: short work") + } + if len(ifac) < 15 { + panic("fourier: short ifac") + } + if n == 1 { + return + } + rfftb1(n, r[:n], work[:n], work[n:2*n], ifac[:15]) +} + +func rfftb1(n int, c, ch, wa []float64, ifac []int) { + nf := ifac[1] + na := false + l1 := 1 + iw := 0 + + for k1 := 1; k1 <= nf; k1++ { + ip := ifac[k1+1] + l2 := ip * l1 + ido := n / l2 + idl1 := ido * l1 + + switch ip { + case 4: + ix2 := iw + ido + ix3 := ix2 + ido + if na { + radb4(ido, l1, ch, c, wa[iw:], wa[ix2:], wa[ix3:]) + } else { + radb4(ido, l1, c, ch, wa[iw:], wa[ix2:], wa[ix3:]) + } + na = !na + case 2: + if na { + radb2(ido, l1, ch, c, wa[iw:]) + } else { + radb2(ido, l1, c, ch, wa[iw:]) + } + na = !na + case 3: + ix2 := iw + ido + if na { + radb3(ido, l1, ch, c, wa[iw:], wa[ix2:]) + } else { + radb3(ido, l1, c, ch, wa[iw:], wa[ix2:]) + } + na = !na + case 5: + ix2 := iw + ido + ix3 := ix2 + ido + ix4 := ix3 + ido + if na { + radb5(ido, l1, ch, c, wa[iw:], wa[ix2:], wa[ix3:], wa[ix4:]) + } else { + radb5(ido, l1, c, ch, wa[iw:], wa[ix2:], wa[ix3:], wa[ix4:]) + } + na = !na + default: + if na { + radbg(ido, ip, l1, idl1, ch, ch, ch, c, c, wa[iw:]) + } else { + radbg(ido, ip, l1, idl1, c, c, c, ch, ch, wa[iw:]) + } + if ido == 1 { + na = !na + } + } + + l1 = l2 + iw += (ip - 1) * ido + } + + if na { + for i := 0; i < n; i++ { + c[i] = ch[i] + } + } +} + +func radb2(ido, l1 int, cc, ch, wa1 []float64) { + cc3 := newThreeArray(ido, 2, l1, cc) + ch3 := newThreeArray(ido, l1, 2, ch) + + for k := 0; k < l1; k++ { + ch3.set(0, k, 0, cc3.at(0, 0, k)+cc3.at(ido-1, 1, k)) + ch3.set(0, k, 1, cc3.at(0, 0, k)-cc3.at(ido-1, 1, k)) + } + + if ido < 2 { + return + } + if ido > 2 { + idp2 := ido + 1 + for k := 0; k < l1; k++ { + for i := 2; i < ido; i += 2 { + ic := idp2 - (i + 1) + + ch3.setCmplx(i-1, k, 0, cc3.atCmplx(i-1, 0, k)+cmplx.Conj(cc3.atCmplx(ic-1, 1, k))) + + t2 := cc3.atCmplx(i-1, 0, k) - cmplx.Conj(cc3.atCmplx(ic-1, 1, k)) + ch3.setCmplx(i-1, k, 1, complex(wa1[i-2], wa1[i-1])*t2) + } + } + + if ido%2 == 1 { + return + } + } + for k := 0; k < l1; k++ { + ch3.set(ido-1, k, 0, 2*cc3.at(ido-1, 0, k)) + ch3.set(ido-1, k, 1, -2*cc3.at(0, 1, k)) + } +} + +func radb3(ido, l1 int, cc, ch, wa1, wa2 []float64) { + const ( + taur = -0.5 + taui = 0.866025403784439 // sqrt(3)/2 + ) + + cc3 := newThreeArray(ido, 3, l1, cc) + ch3 := newThreeArray(ido, l1, 3, ch) + + for k := 0; k < l1; k++ { + tr2 := cc3.at(ido-1, 1, k) + cc3.at(ido-1, 1, k) + cr2 := cc3.at(0, 0, k) + taur*tr2 + ch3.set(0, k, 0, cc3.at(0, 0, k)+tr2) + ci3 := taui * (cc3.at(0, 2, k) + cc3.at(0, 2, k)) + ch3.set(0, k, 1, cr2-ci3) + ch3.set(0, k, 2, cr2+ci3) + } + + if ido == 1 { + return + } + + idp2 := ido + 1 + for k := 0; k < l1; k++ { + for i := 2; i < ido; i += 2 { + ic := idp2 - (i + 1) + + t2 := cc3.atCmplx(i-1, 2, k) + cmplx.Conj(cc3.atCmplx(ic-1, 1, k)) + c2 := cc3.atCmplx(i-1, 0, k) + scale(taur, t2) + ch3.setCmplx(i-1, k, 0, cc3.atCmplx(i-1, 0, k)+t2) + + c3 := scale(taui, cc3.atCmplx(i-1, 2, k)-cmplx.Conj(cc3.atCmplx(ic-1, 1, k))) + d2 := c2 - cmplx.Conj(swap(c3)) + d3 := c2 + cmplx.Conj(swap(c3)) + ch3.setCmplx(i-1, k, 1, complex(wa1[i-2], wa1[i-1])*d2) + ch3.setCmplx(i-1, k, 2, complex(wa2[i-2], wa2[i-1])*d3) + } + } +} + +func radb4(ido, l1 int, cc, ch, wa1, wa2, wa3 []float64) { + cc3 := newThreeArray(ido, 4, l1, cc) + ch3 := newThreeArray(ido, l1, 4, ch) + + for k := 0; k < l1; k++ { + tr1 := cc3.at(0, 0, k) - cc3.at(ido-1, 3, k) + tr2 := cc3.at(0, 0, k) + cc3.at(ido-1, 3, k) + tr3 := cc3.at(ido-1, 1, k) + cc3.at(ido-1, 1, k) + tr4 := cc3.at(0, 2, k) + cc3.at(0, 2, k) + ch3.set(0, k, 0, tr2+tr3) + ch3.set(0, k, 1, tr1-tr4) + ch3.set(0, k, 2, tr2-tr3) + ch3.set(0, k, 3, tr1+tr4) + } + + if ido < 2 { + return + } + if ido > 2 { + idp2 := ido + 1 + for k := 0; k < l1; k++ { + for i := 2; i < ido; i += 2 { + ic := idp2 - (i + 1) + + t1 := cc3.atCmplx(i-1, 0, k) - cmplx.Conj(cc3.atCmplx(ic-1, 3, k)) + t2 := cc3.atCmplx(i-1, 0, k) + cmplx.Conj(cc3.atCmplx(ic-1, 3, k)) + t3 := cc3.atCmplx(i-1, 2, k) + cmplx.Conj(cc3.atCmplx(ic-1, 1, k)) + t4 := swap(cc3.atCmplx(i-1, 2, k) - cmplx.Conj(cc3.atCmplx(ic-1, 1, k))) + ch3.setCmplx(i-1, k, 0, t2+t3) + + c2 := t1 - cmplx.Conj(t4) + c3 := t2 - t3 + c4 := t1 + cmplx.Conj(t4) + ch3.setCmplx(i-1, k, 1, complex(wa1[i-2], wa1[i-1])*c2) + ch3.setCmplx(i-1, k, 2, complex(wa2[i-2], wa2[i-1])*c3) + ch3.setCmplx(i-1, k, 3, complex(wa3[i-2], wa3[i-1])*c4) + } + } + + if ido%2 == 1 { + return + } + } + for k := 0; k < l1; k++ { + tr1 := cc3.at(ido-1, 0, k) - cc3.at(ido-1, 2, k) + ti1 := cc3.at(0, 1, k) + cc3.at(0, 3, k) + tr2 := cc3.at(ido-1, 0, k) + cc3.at(ido-1, 2, k) + ti2 := cc3.at(0, 3, k) - cc3.at(0, 1, k) + ch3.set(ido-1, k, 0, tr2+tr2) + ch3.set(ido-1, k, 1, math.Sqrt2*(tr1-ti1)) + ch3.set(ido-1, k, 2, ti2+ti2) + ch3.set(ido-1, k, 3, -math.Sqrt2*(tr1+ti1)) + } +} + +func radb5(ido, l1 int, cc, ch, wa1, wa2, wa3, wa4 []float64) { + const ( + tr11 = 0.309016994374947 + ti11 = 0.951056516295154 + tr12 = -0.809016994374947 + ti12 = 0.587785252292473 + ) + + cc3 := newThreeArray(ido, 5, l1, cc) + ch3 := newThreeArray(ido, l1, 5, ch) + + for k := 0; k < l1; k++ { + tr2 := cc3.at(ido-1, 1, k) + cc3.at(ido-1, 1, k) + tr3 := cc3.at(ido-1, 3, k) + cc3.at(ido-1, 3, k) + ti4 := cc3.at(0, 4, k) + cc3.at(0, 4, k) + ti5 := cc3.at(0, 2, k) + cc3.at(0, 2, k) + ch3.set(0, k, 0, cc3.at(0, 0, k)+tr2+tr3) + + cr2 := cc3.at(0, 0, k) + tr11*tr2 + tr12*tr3 + + cr3 := cc3.at(0, 0, k) + tr12*tr2 + tr11*tr3 + + ci4 := ti12*ti5 - ti11*ti4 + + ci5 := ti11*ti5 + ti12*ti4 + + ch3.set(0, k, 1, cr2-ci5) + ch3.set(0, k, 2, cr3-ci4) + ch3.set(0, k, 3, cr3+ci4) + ch3.set(0, k, 4, cr2+ci5) + } + + if ido == 1 { + return + } + + idp2 := ido + 1 + for k := 0; k < l1; k++ { + for i := 2; i < ido; i += 2 { + ic := idp2 - (i + 1) + + t2 := cc3.atCmplx(i-1, 2, k) + cmplx.Conj(cc3.atCmplx(ic-1, 1, k)) + t3 := cc3.atCmplx(i-1, 4, k) + cmplx.Conj(cc3.atCmplx(ic-1, 3, k)) + t4 := cc3.atCmplx(i-1, 4, k) - cmplx.Conj(cc3.atCmplx(ic-1, 3, k)) + t5 := cc3.atCmplx(i-1, 2, k) - cmplx.Conj(cc3.atCmplx(ic-1, 1, k)) + ch3.setCmplx(i-1, k, 0, cc3.atCmplx(i-1, 0, k)+t2+t3) + + c2 := cc3.atCmplx(i-1, 0, k) + scale(tr11, t2) + scale(tr12, t3) + c3 := cc3.atCmplx(i-1, 0, k) + scale(tr12, t2) + scale(tr11, t3) + c4 := scale(ti12, t5) - scale(ti11, t4) + c5 := scale(ti11, t5) + scale(ti12, t4) + d2 := c2 - cmplx.Conj(swap(c5)) + d3 := c3 - cmplx.Conj(swap(c4)) + d4 := c3 + cmplx.Conj(swap(c4)) + d5 := c2 + cmplx.Conj(swap(c5)) + ch3.setCmplx(i-1, k, 1, complex(wa1[i-2], wa1[i-1])*d2) + ch3.setCmplx(i-1, k, 2, complex(wa2[i-2], wa2[i-1])*d3) + ch3.setCmplx(i-1, k, 3, complex(wa3[i-2], wa3[i-1])*d4) + ch3.setCmplx(i-1, k, 4, complex(wa4[i-2], wa4[i-1])*d5) + } + } +} + +func radbg(ido, ip, l1, idl1 int, cc, c1, c2, ch, ch2, wa []float64) { + cc3 := newThreeArray(ido, ip, l1, cc) + c13 := newThreeArray(ido, l1, ip, c1) + ch3 := newThreeArray(ido, l1, ip, ch) + c2m := newTwoArray(idl1, ip, c2) + ch2m := newTwoArray(idl1, ip, ch2) + + arg := 2 * math.Pi / float64(ip) + dcp := math.Cos(arg) + dsp := math.Sin(arg) + ipph := (ip + 1) / 2 + nbd := (ido - 1) / 2 + + if ido < l1 { + for i := 0; i < ido; i++ { + for k := 0; k < l1; k++ { + ch3.set(i, k, 0, cc3.at(i, 0, k)) + } + } + } else { + for k := 0; k < l1; k++ { + for i := 0; i < ido; i++ { + ch3.set(i, k, 0, cc3.at(i, 0, k)) + } + } + } + + for j := 1; j < ipph; j++ { + jc := ip - j + j2 := 2 * j + for k := 0; k < l1; k++ { + ch3.set(0, k, j, cc3.at(ido-1, j2-1, k)+cc3.at(ido-1, j2-1, k)) + ch3.set(0, k, jc, cc3.at(0, j2, k)+cc3.at(0, j2, k)) + } + } + + if ido != 1 { + if nbd < l1 { + for j := 1; j < ipph; j++ { + jc := ip - j + j2 := 2 * j + for i := 2; i < ido; i += 2 { + ic := ido - i + for k := 0; k < l1; k++ { + ch3.setCmplx(i-1, k, j, cc3.atCmplx(i-1, j2, k)+cmplx.Conj(cc3.atCmplx(ic-1, j2-1, k))) + ch3.setCmplx(i-1, k, jc, cc3.atCmplx(i-1, j2, k)-cmplx.Conj(cc3.atCmplx(ic-1, j2-1, k))) + } + } + } + } else { + for j := 1; j < ipph; j++ { + jc := ip - j + j2 := 2 * j + for k := 0; k < l1; k++ { + for i := 2; i < ido; i += 2 { + ic := ido - i + ch3.setCmplx(i-1, k, j, cc3.atCmplx(i-1, j2, k)+cmplx.Conj(cc3.atCmplx(ic-1, j2-1, k))) + ch3.setCmplx(i-1, k, jc, cc3.atCmplx(i-1, j2, k)-cmplx.Conj(cc3.atCmplx(ic-1, j2-1, k))) + } + } + } + } + } + + ar1 := 1.0 + ai1 := 0.0 + for l := 1; l < ipph; l++ { + lc := ip - l + ar1h := dcp*ar1 - dsp*ai1 + ai1 = dcp*ai1 + dsp*ar1 + ar1 = ar1h + for ik := 0; ik < idl1; ik++ { + c2m.set(ik, l, ch2m.at(ik, 0)+ar1*ch2m.at(ik, 1)) + c2m.set(ik, lc, ai1*ch2m.at(ik, ip-1)) + } + dc2 := ar1 + ds2 := ai1 + ar2 := ar1 + ai2 := ai1 + for j := 2; j < ipph; j++ { + jc := ip - j + ar2h := dc2*ar2 - ds2*ai2 + ai2 = dc2*ai2 + ds2*ar2 + ar2 = ar2h + for ik := 0; ik < idl1; ik++ { + c2m.add(ik, l, ar2*ch2m.at(ik, j)) + c2m.add(ik, lc, ai2*ch2m.at(ik, jc)) + } + } + } + + for j := 1; j < ipph; j++ { + for ik := 0; ik < idl1; ik++ { + ch2m.add(ik, 0, ch2m.at(ik, j)) + } + } + for j := 1; j < ipph; j++ { + jc := ip - j + for k := 0; k < l1; k++ { + ch3.set(0, k, j, c13.at(0, k, j)-c13.at(0, k, jc)) + ch3.set(0, k, jc, c13.at(0, k, j)+c13.at(0, k, jc)) + } + } + + if ido != 1 { + if nbd < l1 { + for j := 1; j < ipph; j++ { + jc := ip - j + for i := 2; i < ido; i += 2 { + for k := 0; k < l1; k++ { + ch3.setCmplx(i-1, k, j, c13.atCmplx(i-1, k, j)-cmplx.Conj(swap(c13.atCmplx(i-1, k, jc)))) + ch3.setCmplx(i-1, k, jc, c13.atCmplx(i-1, k, j)+cmplx.Conj(swap(c13.atCmplx(i-1, k, jc)))) + } + } + } + } else { + for j := 1; j < ipph; j++ { + jc := ip - j + for k := 0; k < l1; k++ { + for i := 2; i < ido; i += 2 { + ch3.setCmplx(i-1, k, j, c13.atCmplx(i-1, k, j)-cmplx.Conj(swap(c13.atCmplx(i-1, k, jc)))) + ch3.setCmplx(i-1, k, jc, c13.atCmplx(i-1, k, j)+cmplx.Conj(swap(c13.atCmplx(i-1, k, jc)))) + } + } + } + } + } + + if ido == 1 { + return + } + for ik := 0; ik < idl1; ik++ { + c2m.set(ik, 0, ch2m.at(ik, 0)) + } + for j := 1; j < ip; j++ { + for k := 0; k < l1; k++ { + c13.set(0, k, j, ch3.at(0, k, j)) + } + } + + is := -ido - 1 + if nbd > l1 { + for j := 1; j < ip; j++ { + is += ido + for k := 0; k < l1; k++ { + idij := is + for i := 2; i < ido; i += 2 { + idij += 2 + c13.setCmplx(i-1, k, j, complex(wa[idij-1], wa[idij])*ch3.atCmplx(i-1, k, j)) + } + } + } + return + } + for j := 1; j < ip; j++ { + is += ido + idij := is + for i := 2; i < ido; i += 2 { + idij += 2 + for k := 0; k < l1; k++ { + c13.setCmplx(i-1, k, j, complex(wa[idij-1], wa[idij])*ch3.atCmplx(i-1, k, j)) + } + } + } +} diff --git a/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/sinq.go b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/sinq.go new file mode 100644 index 0000000..a2240c0 --- /dev/null +++ b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/sinq.go @@ -0,0 +1,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] + } + } +} diff --git a/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/sint.go b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/sint.go new file mode 100644 index 0000000..a3bf424 --- /dev/null +++ b/vendor/gonum.org/v1/gonum/dsp/fourier/internal/fftpack/sint.go @@ -0,0 +1,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] + } +} diff --git a/vendor/gonum.org/v1/gonum/dsp/fourier/quarter.go b/vendor/gonum.org/v1/gonum/dsp/fourier/quarter.go new file mode 100644 index 0000000..3f640de --- /dev/null +++ b/vendor/gonum.org/v1/gonum/dsp/fourier/quarter.go @@ -0,0 +1,133 @@ +// 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. + +package fourier + +import "gonum.org/v1/gonum/dsp/fourier/internal/fftpack" + +// QuarterWaveFFT implements Fast Fourier Transform for quarter wave data. +type QuarterWaveFFT struct { + work []float64 + ifac [15]int +} + +// NewQuarterWaveFFT returns a QuarterWaveFFT initialized for work on sequences of length n. +func NewQuarterWaveFFT(n int) *QuarterWaveFFT { + var t QuarterWaveFFT + t.Reset(n) + return &t +} + +// Len returns the length of the acceptable input. +func (t *QuarterWaveFFT) Len() int { return len(t.work) / 3 } + +// Reset reinitializes the QuarterWaveFFT for work on sequences of length n. +func (t *QuarterWaveFFT) Reset(n int) { + if 3*n <= cap(t.work) { + t.work = t.work[:3*n] + } else { + t.work = make([]float64, 3*n) + } + fftpack.Cosqi(n, t.work, t.ifac[:]) +} + +// CosCoefficients computes the Fast Fourier Transform of quarter wave data for +// the input sequence, seq, placing the cosine series coefficients in dst and +// returning it. +// This transform is unnormalized; a call to CosCoefficients followed by a call +// to CosSequence will multiply the input sequence by 4*n, where n is the length +// of the sequence. +// +// If the length of seq is not t.Len(), CosCoefficients will panic. +// If dst is nil, a new slice is allocated and returned. If dst is not nil and +// the length of dst does not equal t.Len(), CosCoefficients will panic. +// It is safe to use the same slice for dst and seq. +func (t *QuarterWaveFFT) CosCoefficients(dst, seq []float64) []float64 { + if len(seq) != t.Len() { + panic("fourier: sequence length mismatch") + } + if dst == nil { + dst = make([]float64, t.Len()) + } else if len(dst) != t.Len() { + panic("fourier: destination length mismatch") + } + copy(dst, seq) + fftpack.Cosqf(len(dst), dst, t.work, t.ifac[:]) + return dst +} + +// CosSequence computes the Inverse Fast Fourier Transform of quarter wave data for +// the input cosine series coefficients, coeff, placing the sequence data in dst +// and returning it. +// This transform is unnormalized; a call to CosSequence followed by a call +// to CosCoefficients will multiply the input sequence by 4*n, where n is the length +// of the sequence. +// +// If the length of seq is not t.Len(), CosSequence will panic. +// If dst is nil, a new slice is allocated and returned. If dst is not nil and +// the length of dst does not equal t.Len(), CosSequence will panic. +// It is safe to use the same slice for dst and seq. +func (t *QuarterWaveFFT) CosSequence(dst, coeff []float64) []float64 { + if len(coeff) != t.Len() { + panic("fourier: coefficients length mismatch") + } + if dst == nil { + dst = make([]float64, t.Len()) + } else if len(dst) != t.Len() { + panic("fourier: destination length mismatch") + } + copy(dst, coeff) + fftpack.Cosqb(len(dst), dst, t.work, t.ifac[:]) + return dst +} + +// SinCoefficients computes the Fast Fourier Transform of quarter wave data for +// the input sequence, seq, placing the sine series coefficients in dst and +// returning it. +// This transform is unnormalized; a call to SinCoefficients followed by a call +// to SinSequence will multiply the input sequence by 4*n, where n is the length +// of the sequence. +// +// If the length of seq is not t.Len(), SinCoefficients will panic. +// If dst is nil, a new slice is allocated and returned. If dst is not nil and +// the length of dst does not equal t.Len(), SinCoefficients will panic. +// It is safe to use the same slice for dst and seq. +func (t *QuarterWaveFFT) SinCoefficients(dst, seq []float64) []float64 { + if len(seq) != t.Len() { + panic("fourier: sequence length mismatch") + } + if dst == nil { + dst = make([]float64, t.Len()) + } else if len(dst) != t.Len() { + panic("fourier: destination length mismatch") + } + copy(dst, seq) + fftpack.Sinqf(len(dst), dst, t.work, t.ifac[:]) + return dst +} + +// SinSequence computes the Inverse Fast Fourier Transform of quarter wave data for +// the input sine series coefficients, coeff, placing the sequence data in dst +// and returning it. +// This transform is unnormalized; a call to SinSequence followed by a call +// to SinCoefficients will multiply the input sequence by 4*n, where n is the length +// of the sequence. +// +// If the length of seq is not t.Len(), SinSequence will panic. +// If dst is nil, a new slice is allocated and returned. If dst is not nil and +// the length of dst does not equal t.Len(), SinSequence will panic. +// It is safe to use the same slice for dst and seq. +func (t *QuarterWaveFFT) SinSequence(dst, coeff []float64) []float64 { + if len(coeff) != t.Len() { + panic("fourier: coefficients length mismatch") + } + if dst == nil { + dst = make([]float64, t.Len()) + } else if len(dst) != t.Len() { + panic("fourier: destination length mismatch") + } + copy(dst, coeff) + fftpack.Sinqb(len(dst), dst, t.work, t.ifac[:]) + return dst +} diff --git a/vendor/gonum.org/v1/gonum/dsp/fourier/radix24.go b/vendor/gonum.org/v1/gonum/dsp/fourier/radix24.go new file mode 100644 index 0000000..5641332 --- /dev/null +++ b/vendor/gonum.org/v1/gonum/dsp/fourier/radix24.go @@ -0,0 +1,343 @@ +// Copyright ©2020 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. + +package fourier + +import ( + "math" + "math/bits" + "math/cmplx" +) + +// CoefficientsRadix2 computes the Fourier coefficients of the input +// sequence, converting the time series in seq into the frequency spectrum, +// in place and returning it. This transform is unnormalized; a call to +// CoefficientsRadix2 followed by a call of SequenceRadix2 will multiply the +// input sequence by the length of the sequence. +// +// CoefficientsRadix2 does not allocate, requiring that FFT twiddle factors +// be calculated lazily. For performance reasons, this is done by successive +// multiplication, so numerical accuracies can accumulate for large inputs. +// If accuracy is needed, the FFT or CmplxFFT types should be used. +// +// If the length of seq is not an integer power of 2, CoefficientsRadix2 will +// panic. +func CoefficientsRadix2(seq []complex128) []complex128 { + x := seq + switch len(x) { + default: + if bits.OnesCount(uint(len(x))) != 1 { + panic("fourier: radix-2 fft called with non-power 2 length") + } + + case 0, 1: + return x + + case 2: + x[0], x[1] = + x[0]+x[1], + x[0]-x[1] + return x + + case 4: + t := x[1] + x[3] + u := x[2] + v := negi(x[1] - x[3]) + x[0], x[1], x[2], x[3] = + x[0]+u+t, + x[0]-u+v, + x[0]+u-t, + x[0]-u-v + return x + } + + bitReversePermute(x) + + for k := 0; k < len(x); k += 4 { + t := x[k+2] + x[k+3] + u := x[k+1] + v := negi(x[k+2] - x[k+3]) + x[k], x[k+1], x[k+2], x[k+3] = + x[k]+u+t, + x[k]-u+v, + x[k]+u-t, + x[k]-u-v + } + + for m := 4; m < len(x); m *= 2 { + f := swap(complex(math.Sincos(-math.Pi / float64(m)))) + for k := 0; k < len(x); k += 2 * m { + w := 1 + 0i + for j := 0; j < m; j++ { + i := j + k + + u := w * x[i+m] + x[i], x[i+m] = + x[i]+u, + x[i]-u + + w *= f + } + } + } + + return x +} + +// bitReversePermute performs a bit-reversal permutation on x. +func bitReversePermute(x []complex128) { + if len(x) < 2 || bits.OnesCount(uint(len(x))) != 1 { + panic("fourier: invalid bitReversePermute call") + } + lz := bits.LeadingZeros(uint(len(x) - 1)) + i := 0 + for ; i < len(x)/2; i++ { + j := int(bits.Reverse(uint(i)) >> lz) + if i < j { + x[i], x[j] = x[j], x[i] + } + } + for i++; i < len(x); i += 2 { + j := int(bits.Reverse(uint(i)) >> lz) + if i < j { + x[i], x[j] = x[j], x[i] + } + } +} + +// SequenceRadix2 computes the real periodic sequence from the Fourier +// coefficients, converting the frequency spectrum in coeff into a time +// series, in place and returning it. This transform is unnormalized; a +// call to CoefficientsRadix2 followed by a call of SequenceRadix2 will +// multiply the input sequence by the length of the sequence. +// +// SequenceRadix2 does not allocate, requiring that FFT twiddle factors +// be calculated lazily. For performance reasons, this is done by successive +// multiplication, so numerical accuracies can accumulate for large inputs. +// If accuracy is needed, the FFT or CmplxFFT types should be used. +// +// If the length of coeff is not an integer power of 2, SequenceRadix2 +// will panic. +func SequenceRadix2(coeff []complex128) []complex128 { + x := coeff + for i, j := 1, len(x)-1; i < j; i, j = i+1, j-1 { + x[i], x[j] = x[j], x[i] + } + + CoefficientsRadix2(x) + return x +} + +// PadRadix2 returns the values in x in a slice that is an integer +// power of 2 long. If x already has an integer power of 2 length +// it is returned unaltered. +func PadRadix2(x []complex128) []complex128 { + if len(x) == 0 { + return x + } + b := bits.Len(uint(len(x))) + if len(x) == 1<<(b-1) { + return x + } + p := make([]complex128, 1<> lz) + if i < j { + x[i], x[j] = x[j], x[i] + } + } + for i++; i < len(x); i += 2 { + j := int(reversePairs(uint(i)) >> lz) + if i < j { + x[i], x[j] = x[j], x[i] + } + } +} + +// SequenceRadix4 computes the real periodic sequence from the Fourier +// coefficients, converting the frequency spectrum in coeff into a time +// series, in place and returning it. This transform is unnormalized; a +// call to CoefficientsRadix4 followed by a call of SequenceRadix4 will +// multiply the input sequence by the length of the sequence. +// +// SequenceRadix4 does not allocate, requiring that FFT twiddle factors +// be calculated lazily. For performance reasons, this is done by successive +// multiplication, so numerical accuracies can accumulate for large inputs. +// If accuracy is needed, the FFT or CmplxFFT types should be used. +// +// If the length of coeff is not an integer power of 4, SequenceRadix4 +// will panic. +func SequenceRadix4(coeff []complex128) []complex128 { + x := coeff + for i, j := 1, len(x)-1; i < j; i, j = i+1, j-1 { + x[i], x[j] = x[j], x[i] + } + + CoefficientsRadix4(x) + return x +} + +// PadRadix4 returns the values in x in a slice that is an integer +// power of 4 long. If x already has an integer power of 4 length +// it is returned unaltered. +func PadRadix4(x []complex128) []complex128 { + if len(x) == 0 { + return x + } + b := bits.Len(uint(len(x))) + if len(x) == 1<<(b-1) && b&0x1 == 1 { + return x + } + p := make([]complex128, 1<<((b+1)&^0x1)) + copy(p, x) + return p +} + +// TrimRadix4 returns the largest slice of x that is has an integer +// power of 4 length, and a slice holding the remaining elements. +func TrimRadix4(x []complex128) (even, remains []complex128) { + if len(x) == 0 { + return x, nil + } + n := 1 << ((bits.Len(uint(len(x))) - 1) &^ 0x1) + return x[:n], x[n:] +} + +// reversePairs returns the value of x with its bit pairs in reversed order. +func reversePairs(x uint) uint { + if bits.UintSize == 32 { + return uint(reversePairs32(uint32(x))) + } + return uint(reversePairs64(uint64(x))) +} + +const ( + m1 = 0x3333333333333333 + m2 = 0x0f0f0f0f0f0f0f0f +) + +// reversePairs32 returns the value of x with its bit pairs in reversed order. +func reversePairs32(x uint32) uint32 { + const m = 1<<32 - 1 + x = x>>2&(m1&m) | x&(m1&m)<<2 + x = x>>4&(m2&m) | x&(m2&m)<<4 + return bits.ReverseBytes32(x) +} + +// reversePairs64 returns the value of x with its bit pairs in reversed order. +func reversePairs64(x uint64) uint64 { + const m = 1<<64 - 1 + x = x>>2&(m1&m) | x&(m1&m)<<2 + x = x>>4&(m2&m) | x&(m2&m)<<4 + return bits.ReverseBytes64(x) +} + +func negi(c complex128) complex128 { + return cmplx.Conj(swap(c)) +} + +func swap(c complex128) complex128 { + return complex(imag(c), real(c)) +} diff --git a/vendor/gonum.org/v1/gonum/dsp/fourier/sincos.go b/vendor/gonum.org/v1/gonum/dsp/fourier/sincos.go new file mode 100644 index 0000000..fb765e9 --- /dev/null +++ b/vendor/gonum.org/v1/gonum/dsp/fourier/sincos.go @@ -0,0 +1,112 @@ +// 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. + +package fourier + +import "gonum.org/v1/gonum/dsp/fourier/internal/fftpack" + +// DCT implements Discrete Cosine Transform for real sequences. +type DCT struct { + work []float64 + ifac [15]int +} + +// NewDCT returns a DCT initialized for work on sequences of length n. +// NewDCT will panic is n is not greater than 1. +func NewDCT(n int) *DCT { + var t DCT + t.Reset(n) + return &t +} + +// Len returns the length of the acceptable input. +func (t *DCT) Len() int { return len(t.work) / 3 } + +// Reset reinitializes the DCT for work on sequences of length n. +// Reset will panic is n is not greater than 1. +func (t *DCT) Reset(n int) { + if n <= 1 { + panic("fourier: n less than 2") + } + if 3*n <= cap(t.work) { + t.work = t.work[:3*n] + } else { + t.work = make([]float64, 3*n) + } + fftpack.Costi(n, t.work, t.ifac[:]) +} + +// Transform computes the Discrete Fourier Cosine Transform of +// the input data, src, placing the result in dst and returning it. +// This transform is unnormalized; a call to Transform followed by +// another call to Transform will multiply the input sequence by 2*(n-1), +// where n is the length of the sequence. +// +// If the length of src is not t.Len(), Transform will panic. +// If dst is nil, a new slice is allocated and returned. If dst is not nil and +// the length of dst does not equal t.Len(), FFT will panic. +// It is safe to use the same slice for dst and src. +func (t *DCT) Transform(dst, src []float64) []float64 { + if len(src) != t.Len() { + panic("fourier: sequence length mismatch") + } + if dst == nil { + dst = make([]float64, t.Len()) + } else if len(dst) != t.Len() { + panic("fourier: destination length mismatch") + } + copy(dst, src) + fftpack.Cost(len(dst), dst, t.work, t.ifac[:]) + return dst +} + +// DST implements Discrete Sine Transform for real sequences. +type DST struct { + work []float64 + ifac [15]int +} + +// NewDST returns a DST initialized for work on sequences of length n. +func NewDST(n int) *DST { + var t DST + t.Reset(n) + return &t +} + +// Len returns the length of the acceptable input. +func (t *DST) Len() int { return (2*len(t.work)+1)/5 - 1 } + +// Reset reinitializes the DCT for work on sequences of length n. +func (t *DST) Reset(n int) { + if 5*(n+1)/2 <= cap(t.work) { + t.work = t.work[:5*(n+1)/2] + } else { + t.work = make([]float64, 5*(n+1)/2) + } + fftpack.Sinti(n, t.work, t.ifac[:]) +} + +// Transform computes the Discrete Fourier Sine Transform of the input +// data, src, placing the result in dst and returning it. +// This transform is unnormalized; a call to Transform followed by +// another call to Transform will multiply the input sequence by 2*(n-1), +// where n is the length of the sequence. +// +// If the length of src is not t.Len(), Transform will panic. +// If dst is nil, a new slice is allocated and returned. If dst is not nil and +// the length of dst does not equal t.Len(), FFT will panic. +// It is safe to use the same slice for dst and src. +func (t *DST) Transform(dst, src []float64) []float64 { + if len(src) != t.Len() { + panic("fourier: sequence length mismatch") + } + if dst == nil { + dst = make([]float64, t.Len()) + } else if len(dst) != t.Len() { + panic("fourier: destination length mismatch") + } + copy(dst, src) + fftpack.Sint(len(dst), dst, t.work, t.ifac[:]) + return dst +}