2 Copyright (C) 2008 Paul Davis
3 Author: Sampo Savolainen
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2 of the License, or
8 (at your option) any later version.
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
26 using namespace GTKArdour
;
28 FFT::FFT(uint32_t windowSize
)
29 : _window_size(windowSize
),
30 _data_size(_window_size
/2),
34 _fftInput
= (float *) fftwf_malloc(sizeof(float) * _window_size
);
36 _fftOutput
= (float *) fftwf_malloc(sizeof(float) * _window_size
);
38 _power_at_bin
= (float *) malloc(sizeof(float) * _data_size
);
39 _phase_at_bin
= (float *) malloc(sizeof(float) * _data_size
);
41 _plan
= fftwf_plan_r2r_1d(_window_size
, _fftInput
, _fftOutput
, FFTW_R2HC
, FFTW_ESTIMATE
);
49 memset(_power_at_bin
, 0, sizeof(float) * _data_size
);
50 memset(_phase_at_bin
, 0, sizeof(float) * _data_size
);
56 FFT::analyze(ARDOUR::Sample
*input
, WindowingType windowing_type
)
60 memcpy(_fftInput
, input
, sizeof(float) * _window_size
);
62 if (windowing_type
== HANN
) {
63 float *window
= get_hann_window();
64 for (uint32_t i
= 0; i
< _window_size
; i
++) {
65 _fftInput
[i
] *= window
[i
];
71 _power_at_bin
[0] += _fftOutput
[0] * _fftOutput
[0];
72 _phase_at_bin
[0] += 0.0;
77 #define Re (_fftOutput[i])
78 #define Im (_fftOutput[_window_size-i])
79 for (uint32_t i
=1; i
< _data_size
- 1; i
++) {
81 power
= (Re
* Re
) + (Im
* Im
);
82 phase
= atanf(Im
/ Re
);
84 if (Re
< 0.0 && Im
> 0.0) {
86 } else if (Re
< 0.0 && Im
< 0.0) {
90 _power_at_bin
[i
] += power
;
91 _phase_at_bin
[i
] += phase
;
100 if (_iterations
> 1) {
101 for (uint32_t i
=0; i
< _data_size
- 1; i
++) {
102 _power_at_bin
[i
] /= (float)_iterations
;
103 _phase_at_bin
[i
] /= (float)_iterations
;
110 FFT::get_hann_window()
116 _hann_window
= (float *) malloc(sizeof(float) * _window_size
);
120 for (uint32_t i
=0; i
< _window_size
; i
++) {
121 _hann_window
[i
]=0.81f
* ( 0.5f
- (0.5f
* (float) cos(2.0f
* M_PI
* (float)i
/ (float)(_window_size
))));
122 sum
+= _hann_window
[i
];
125 double isum
= 1.0 / sum
;
127 for (uint32_t i
=0; i
< _window_size
; i
++) {
128 _hann_window
[i
] *= isum
;
140 fftwf_destroy_plan(_plan
);