File: fft.cc

package info (click to toggle)
din 5.2.1-3
  • links: PTS, VCS
  • area: main
  • in suites: jessie, jessie-kfreebsd
  • size: 2,152 kB
  • ctags: 2,490
  • sloc: cpp: 9,369; sh: 6,563; ansic: 2,977; tcl: 1,770; makefile: 285
file content (101 lines) | stat: -rw-r--r-- 2,606 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
/*
 * This file is part of din.
 *
 * din is copyright (c) 2006 - 2012 S Jagannathan <jag@dinisnoise.org>
 * For more information, please visit http://dinisnoise.org
 *
 * din is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 2 of the License, or
 * (at your option) any later version.
 *
 * din is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
 * for more details.
 *
 * You should have received a copy of the GNU General Public License along
 * with din.  If not, see <http://www.gnu.org/licenses/>.
 *
*/
#include "fft.h"
#include "multi_curve.h"
#include "solver.h"
#include "math.h"
#include "ansi_color_codes.h"

#include <algorithm>
#include <iostream>
using namespace std;

fft::fft () : lev ("fft_levels") {

  n = 20; // 2 x number of harmonics
  halfn = n / 2;
  dx = 1. / n;
  x = 0;

  // use fftw library for fft
  in = (float *) fftwf_malloc (sizeof(float) * n);
  out = (fftwf_complex *) fftwf_malloc (sizeof (fftwf_complex) * (halfn + 1));

  extern string dotdin;
  string fname = dotdin + "fftwf_wisdom";
  FILE* wis = fopen (fname.c_str(), "r");
  if (wis) {
    cout << DOING << "+++ using saved FFT wisdom +++" << ENDL;
    fftwf_import_wisdom_from_file (wis);
    fclose (wis);
  }

  plan = fftwf_plan_dft_r2c_1d (n, in, out, FFTW_MEASURE);

  cout << PASS << "+++ FFT setup complete +++ " << ENDL;
}

fft::~fft () {
    extern string dotdin;
    string fname = dotdin + "fftwf_wisdom";
    FILE* wis = fopen (fname.c_str(), "w");
    if (wis) {
      fftwf_export_wisdom_to_file (wis);
      fclose (wis);
    }
  fftwf_destroy_plan (plan);
  fftwf_free (in);
  fftwf_free (out);
}

void fft::go (multi_curve* crv) {

  // do fft on a bezier multi curve

  solver s (crv);

  x = 0;
  s (x, dx, n, in);

  fftwf_execute (plan);

  harms.clear ();
  for (int i = 1, j = halfn +1; i < j; ++i) {
    float x = out[i][0];
    float y = out[i][1];
    float x2 = x*x;
    float y2 = y*y;
    harms.push_back (sqrt (x2 + y2));
    // ignore phase (atan2 (y2/x2)) for now, but have to add in later revision.
  }

  float m0 = *min_element (harms.begin(), harms.end());
  float m1 = *max_element (harms.begin(), harms.end());
  float mm = m1 - m0;
  if (m1 == m0) return;

  for (int i = 0, j = harms.size(); i < j; ++i) {
    float dm = harms[i] - m0;
    float y = dm / mm;
    lev.set (i, y);
  }

}