Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactor FFT helpers #5309

Merged
merged 9 commits into from
Mar 11, 2020
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 0 additions & 12 deletions include/fft_helpers.h
Original file line number Diff line number Diff line change
Expand Up @@ -103,16 +103,4 @@ int LMMS_EXPORT absspec(const fftwf_complex *complex_buffer, float *absspec_buff
int LMMS_EXPORT compressbands(const float * _absspec_buffer, float * _compressedband,
int _num_old, int _num_new, int _bottom, int _top);


int LMMS_EXPORT calc13octaveband31(float * _absspec_buffer, float * _subbands,
int _num_spec, float _max_frequency);


/** Compute power of finite time sequence.
* Take care num_values is length of timesignal[].
*
* @return power on success, else -1
*/
float LMMS_EXPORT signalpower(const float *timesignal, int num_values);

#endif
127 changes: 20 additions & 107 deletions src/core/fft_helpers.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -35,13 +35,12 @@
*/
float maximum(const float *abs_spectrum, unsigned int spec_size)
{
float maxi = 0;
unsigned int i;

if (abs_spectrum == NULL) {return -1;}
if (spec_size <= 0) {return -1;}
if (spec_size == 0) {return -1;}

float maxi = 0;

for (i = 0; i < spec_size; i++)
for (unsigned int i = 0; i < spec_size; i++)
{
if (abs_spectrum[i] > maxi) {maxi = abs_spectrum[i];}
}
Expand All @@ -61,13 +60,11 @@ float maximum(const std::vector<float> &abs_spectrum)
*/
int normalize(const float *abs_spectrum, float *norm_spectrum, unsigned int bin_count, unsigned int block_size)
{
int i;

if (abs_spectrum == NULL || norm_spectrum == NULL) {return -1;}
if (bin_count == 0 || block_size == 0) {return -1;}

block_size /= 2;
for (i = 0; i < bin_count; i++)
for (unsigned int i = 0; i < bin_count; i++)
{
norm_spectrum[i] = abs_spectrum[i] / block_size;
}
Expand All @@ -89,9 +86,9 @@ int normalize(const std::vector<float> &abs_spectrum, std::vector<float> &norm_s
*/
int notEmpty(const std::vector<float> &spectrum)
{
for (int i = 0; i < spectrum.size(); i++)
for (float s : spectrum)
{
if (spectrum[i] != 0) {return 1;}
if (s != 0) {return 1;}
}
return 0;
}
Expand All @@ -103,22 +100,21 @@ int notEmpty(const std::vector<float> &spectrum)
*/
int precomputeWindow(float *window, unsigned int length, FFT_WINDOWS type, bool normalized)
{
unsigned int i;
if (window == NULL) {return -1;}

float gain = 0;
float a0;
float a1;
float a2;
float a3;

if (window == NULL) {return -1;}

// constants taken from
// https://en.wikipedia.org/wiki/Window_function#AList_of_window_functions
switch (type)
{
default:
case RECTANGULAR:
for (i = 0; i < length; i++) {window[i] = 1.0;}
for (unsigned int i = 0; i < length; i++) {window[i] = 1.0;}
gain = 1;
return 0;
case BLACKMAN_HARRIS:
Expand All @@ -142,7 +138,7 @@ int precomputeWindow(float *window, unsigned int length, FFT_WINDOWS type, bool
}

// common computation for cosine-sum based windows
for (i = 0; i < length; i++)
for (unsigned int i = 0; i < length; i++)
{
window[i] = (a0 - a1 * cos(2 * F_PI * i / ((float)length - 1.0))
+ a2 * cos(4 * F_PI * i / ((float)length - 1.0))
Expand All @@ -152,7 +148,7 @@ int precomputeWindow(float *window, unsigned int length, FFT_WINDOWS type, bool

// apply amplitude correction
gain /= (float) length;
for (i = 0; i < length; i++) {window[i] /= gain;}
for (unsigned int i = 0; i < length; i++) {window[i] /= gain;}

return 0;
}
Expand All @@ -166,12 +162,10 @@ int precomputeWindow(float *window, unsigned int length, FFT_WINDOWS type, bool
*/
int absspec(const fftwf_complex *complex_buffer, float *absspec_buffer, unsigned int compl_length)
{
int i;

if (complex_buffer == NULL || absspec_buffer == NULL) {return -1;}
if (compl_length <= 0) {return -1;}
if (compl_length == 0) {return -1;}

for (i = 0; i < compl_length; i++)
for (unsigned int i = 0; i < compl_length; i++)
{
absspec_buffer[i] = (float)sqrt(complex_buffer[i][0] * complex_buffer[i][0]
+ complex_buffer[i][1] * complex_buffer[i][1]);
Expand All @@ -189,113 +183,32 @@ int absspec(const fftwf_complex *complex_buffer, float *absspec_buffer, unsigned
*/
int compressbands(const float *absspec_buffer, float *compressedband, int num_old, int num_new, int bottom, int top)
{
float ratio;
int i, usefromold;
float j;
float j_min, j_max;

if (absspec_buffer == NULL || compressedband == NULL) {return -1;}
if (num_old < num_new) {return -1;}
if (num_old <= 0 || num_new <= 0) {return -1;}
if (bottom < 0) {bottom = 0;}
if (top >= num_old) {top = num_old - 1;}

usefromold = num_old - (num_old - top) - bottom;
int usefromold = num_old - (num_old - top) - bottom;

ratio = (float)usefromold / (float)num_new;
float ratio = (float)usefromold / (float)num_new;

// for each new subband
for (i = 0; i < num_new; i++)
for (int i = 0; i < num_new; i++)
{
compressedband[i] = 0;

j_min = (i * ratio) + bottom;
float j_min = (i * ratio) + bottom;

if (j_min < 0) {j_min = bottom;}

j_max = j_min + ratio;
float j_max = j_min + ratio;

for (j = (int)j_min; j <= j_max; j++)
for (float j = (int)j_min; j <= j_max; j++)
{
compressedband[i] += absspec_buffer[(int)j];
}
}

return 0;
}


int calc13octaveband31(float *absspec_buffer, float *subbands, int num_spec, float max_frequency)
{
static const int onethirdoctavecenterfr[] = {20, 25, 31, 40, 50, 63, 80, 100, 125, 160, 200, 250, 315, 400, 500, 630, 800, 1000, 1250, 1600, 2000, 2500, 3150, 4000, 5000, 6300, 8000, 10000, 12500, 16000, 20000};
int i, j;
float f_min, f_max, frequency, bandwidth;
int j_min, j_max = 0;
float fpower;

if (absspec_buffer == NULL || subbands == NULL) {return -1;}
if (num_spec < 31) {return -1;}
if (max_frequency <= 0) {return -1;}

/*** energy ***/
fpower = 0;
for (i = 0; i < num_spec; i++)
{
absspec_buffer[i] = (absspec_buffer[i] * absspec_buffer[i]) / FFT_BUFFER_SIZE;
fpower = fpower + (2 * absspec_buffer[i]);
}
fpower = fpower - (absspec_buffer[0]); //dc not mirrored

/*** for each subband: sum up power ***/
for (i = 0; i < 31; i++)
{
subbands[i] = 0;

// calculate bandwidth for subband
frequency = onethirdoctavecenterfr[i];

bandwidth = (pow(2, 1.0/3.0)-1) * frequency;

f_min = frequency - bandwidth / 2.0;
f_max = frequency + bandwidth / 2.0;

j_min = (int)(f_min / max_frequency * (float)num_spec);
j_max = (int)(f_max / max_frequency * (float)num_spec);


if (j_min < 0 || j_max < 0)
{
fprintf(stderr, "Error: calc13octaveband31() in fft_helpers.cpp line %d failed.\n", __LINE__);
return -1;
}

for (j = j_min; j <= j_max; j++)
{
if (j_max < num_spec) {subbands[i] += absspec_buffer[j];}
}

} //for

return 0;
}

/* Compute power of finite time sequence.
* Take care num_values is length of timesignal[]
*
* return power on success, else -1
*/
float signalpower(const float *timesignal, int num_values)
{
if (num_values <= 0) {return -1;}

if (timesignal == NULL) {return -1;}

float power = 0;
for (int i = 0; i < num_values; i++)
{
power += timesignal[i] * timesignal[i];
}

return power;
}