1 /*
2 * fft_helpers.cpp - some functions around FFT analysis
3 *
4 * Copyright (c) 2008-2012 Tobias Doerffel <tobydox/at/users.sourceforge.net>
5 *
6 * This file is part of LMMS - https://lmms.io
7 *
8 * This program is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public
10 * License as published by the Free Software Foundation; either
11 * version 2 of the License, or (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 * General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public
19 * License along with this program (see COPYING); if not, write to the
20 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
21 * Boston, MA 02110-1301 USA.
22 *
23 */
24
25
26 #include "fft_helpers.h"
27
28 #include <cmath>
29 #include "lmms_constants.h"
30
31 /* returns biggest value from abs_spectrum[spec_size] array
32
33 returns -1 on error
34 */
maximum(float * abs_spectrum,unsigned int spec_size)35 float maximum(float *abs_spectrum, unsigned int spec_size)
36 {
37 float maxi=0;
38 unsigned int i;
39
40 if ( abs_spectrum==NULL )
41 return -1;
42
43 if (spec_size<=0)
44 return -1;
45
46 for ( i=0; i<spec_size; i++ )
47 {
48 if ( abs_spectrum[i]>maxi )
49 maxi=abs_spectrum[i];
50 }
51
52 return maxi;
53 }
54
55
56 /* apply hanning or hamming window to channel
57
58 returns -1 on error */
hanming(float * timebuffer,int length,WINDOWS type)59 int hanming(float *timebuffer, int length, WINDOWS type)
60 {
61 int i;
62 float alpha;
63
64 if ( (timebuffer==NULL)||(length<=0) )
65 return -1;
66
67 switch (type)
68 {
69 case HAMMING: alpha=0.54; break;
70 case HANNING:
71 default: alpha=0.5; break;
72 }
73
74 for ( i=0; i<length; i++ )
75 {
76 timebuffer[i]=timebuffer[i]*(alpha+(1-alpha)*cos(2*F_PI*i/((float)length-1.0)));
77 }
78
79 return 0;
80 }
81
82
83 /* compute absolute values of complex_buffer, save to absspec_buffer
84 take care that - compl_len is not bigger than complex_buffer!
85 - absspec buffer is big enough!
86
87 returns 0 on success, else -1 */
absspec(fftwf_complex * complex_buffer,float * absspec_buffer,int compl_length)88 int absspec(fftwf_complex *complex_buffer, float *absspec_buffer, int compl_length)
89 {
90 int i;
91
92 if ( (complex_buffer==NULL)||(absspec_buffer==NULL) )
93 return -1;
94 if ( compl_length<=0 )
95 return -1;
96
97 for (i=0; i<compl_length; i++)
98 {
99 absspec_buffer[i]=(float )sqrt(complex_buffer[i][0]*complex_buffer[i][0] + complex_buffer[i][1]*complex_buffer[i][1]);
100 }
101
102 return 0;
103 }
104
105
106 /* build fewer subbands from many absolute spectrum values
107 take care that - compressedbands[] array num_new elements long
108 - num_old > num_new
109
110 returns 0 on success, else -1 */
compressbands(float * absspec_buffer,float * compressedband,int num_old,int num_new,int bottom,int top)111 int compressbands(float *absspec_buffer, float *compressedband, int num_old, int num_new, int bottom, int top)
112 {
113 float ratio;
114 int i, usefromold;
115 float j;
116 float j_min, j_max;
117
118 if ( (absspec_buffer==NULL)||(compressedband==NULL) )
119 return -1;
120
121 if ( num_old<num_new )
122 return -1;
123
124 if ( (num_old<=0)||(num_new<=0) )
125 return -1;
126
127 if ( bottom<0 )
128 bottom=0;
129
130 if ( top>=num_old )
131 top=num_old-1;
132
133 usefromold=num_old-(num_old-top)-bottom;
134
135 ratio=(float)usefromold/(float)num_new;
136
137 // for each new subband
138 for ( i=0; i<num_new; i++ )
139 {
140 compressedband[i]=0;
141
142 j_min=(i*ratio)+bottom;
143
144 if ( j_min<0 )
145 j_min=bottom;
146
147 j_max=j_min+ratio;
148
149 for ( j=(int)j_min; j<=j_max; j++ )
150 {
151 compressedband[i]+=absspec_buffer[(int)j];
152 }
153 }
154
155 return 0;
156 }
157
158
calc13octaveband31(float * absspec_buffer,float * subbands,int num_spec,float max_frequency)159 int calc13octaveband31(float *absspec_buffer, float *subbands, int num_spec, float max_frequency)
160 {
161 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};
162 int i, j;
163 float f_min, f_max, frequency, bandwidth;
164 int j_min, j_max=0;
165 float fpower;
166
167
168 if ( (absspec_buffer==NULL)||(subbands==NULL) )
169 return -1;
170
171 if ( num_spec<31 )
172 return -1;
173
174 if ( max_frequency<=0 )
175 return -1;
176
177 /*** energy ***/
178 fpower=0;
179 for ( i=0; i<num_spec; i++ )
180 {
181 absspec_buffer[i]=(absspec_buffer[i]*absspec_buffer[i])/FFT_BUFFER_SIZE;
182 fpower=fpower+(2*absspec_buffer[i]);
183 }
184 fpower=fpower-(absspec_buffer[0]); //dc not mirrored
185
186
187 /*** for each subband: sum up power ***/
188 for ( i=0; i<31; i++ )
189 {
190 subbands[i]=0;
191
192 // calculate bandwidth for subband
193 frequency=onethirdoctavecenterfr[i];
194
195 bandwidth=(pow(2, 1.0/3.0)-1)*frequency;
196
197 f_min=frequency-bandwidth/2.0;
198 f_max=frequency+bandwidth/2.0;
199
200 j_min=(int)(f_min/max_frequency*(float)num_spec);
201
202 j_max=(int)(f_max/max_frequency*(float)num_spec);
203
204
205 if ( (j_min<0)||(j_max<0) )
206 {
207 fprintf(stderr, "Error: calc13octaveband31() in fft_helpers.cpp line %d failed.\n", __LINE__);
208 return -1;
209 }
210
211 for ( j=j_min; j<=j_max; j++ )
212 {
213 if( j_max<num_spec )
214 subbands[i]+=absspec_buffer[j];
215 }
216
217 } //for
218
219
220 return 0;
221 }
222
223 /* compute power of finite time sequence
224 take care num_values is length of timesignal[]
225
226 returns power on success, else -1 */
signalpower(float * timesignal,int num_values)227 float signalpower(float *timesignal, int num_values)
228 {
229 if ( num_values<=0 )
230 return -1;
231
232 if( timesignal==NULL )
233 return -1;
234
235 float power=0;
236 for ( int i=0; i<num_values; i++ )
237 {
238 power+=timesignal[i]*timesignal[i];
239 }
240
241 return power;
242 }
243
244