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);
85 phase = atanf (Im / Re);
86 if (Re < 0.0 && Im > 0.0) {
88 } else if (Re < 0.0 && Im < 0.0) {
93 _power_at_bin[i] += power;
94 _phase_at_bin[i] += phase;
103 if (_iterations > 1) {
104 for (uint32_t i=0; i < _data_size - 1; i++) {
105 _power_at_bin[i] /= (float)_iterations;
106 _phase_at_bin[i] /= (float)_iterations;
113 FFT::get_hann_window ()
119 _hann_window = (float*) malloc (sizeof (float) * _window_size);
123 for (uint32_t i = 0; i < _window_size; ++i) {
124 _hann_window[i] = 0.81f * (0.5f - (0.5f * (float) cos (2.0f * M_PI * (float)i / (float)(_window_size))));
125 sum += _hann_window[i];
128 double isum = 1.0 / sum;
130 for (uint32_t i = 0; i < _window_size; ++i) {
131 _hann_window[i] *= isum;
142 fftwf_destroy_plan(_plan);
145 fftwf_free(_fftOutput);
146 fftwf_free(_fftInput);