Marsyas
0.6.0-alpha
|
00001 /* 00002 ** Copyright (C) 1998-2005 George Tzanetakis <gtzan@cs.uvic.ca> 00003 ** 00004 ** This program is free software; you can redistribute it and/or modify 00005 ** it under the terms of the GNU General Public License as published by 00006 ** the Free Software Foundation; either version 2 of the License, or 00007 ** (at your option) any later version. 00008 ** 00009 ** This program is distributed in the hope that it will be useful, 00010 ** but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00012 ** GNU General Public License for more details. 00013 ** 00014 ** You should have received a copy of the GNU General Public License 00015 ** along with this program; if not, write to the Free Software 00016 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 00017 */ 00018 00019 #include <marsyas/Conversions.h> 00020 #include <marsyas/realvec.h> 00021 00022 using std::ostringstream; 00023 using namespace Marsyas; 00024 00025 mrs_real 00026 Marsyas::pitch2hertz(mrs_real pitch) { 00027 return mrs_real(440.0 * pow(2.0, ((pitch - 69.0) / 12.0))); 00028 } 00029 00030 mrs_real 00031 Marsyas::hertz2pitch(mrs_real hz) { 00032 return (hz == 0.0) ? (mrs_real)0.0 : (mrs_real)(69.0 + 12.0 * (log(hz/440.0)/log(2.0))); 00033 } 00034 00035 mrs_real 00036 Marsyas::samples2hertz(mrs_natural samples, mrs_real srate) { 00037 return (samples <= 0) ? (mrs_real)0.0 : (mrs_real) (srate * 1.0) / (samples); 00038 } 00039 00040 mrs_real 00041 Marsyas::samples2hertz(mrs_real samples, mrs_real srate) { 00042 return (samples <= 0.001) ? (mrs_real)0.0 : (mrs_real) (srate * 1.0) / (samples); 00043 } 00044 00045 00046 mrs_natural 00047 Marsyas::hertz2samples(mrs_real hz, mrs_real srate) { 00048 return (hz == 0.0) ? (mrs_natural)0 : (mrs_natural) (srate / hz); 00049 } 00050 00051 /* convert a string representing time to number of samples base on the 00052 given sample rate. Format "123.456#" where # is the time division. 00053 Valid time divisions: { h, m, s, ms, us }. 00054 On a format error, 00055 Errors: -1 is returned. ie more than 1 decimal point, invalid time 00056 division. 00057 */ 00058 mrs_natural 00059 Marsyas::time2samples(mrs_string time, mrs_real srate) { 00060 //example times: { "10us", "10ms", "10s", "10m", "10h" } 00061 if (time=="") { return 0; } 00062 // calculate time value 00063 mrs_real samples=0; 00064 int i=0; 00065 int len=(int)time.length(); 00066 bool decimal_point=false; 00067 mrs_real divisor = 10.0; 00068 for (i=0; i<len && (time[i]=='.' || (time[i]>='0' && time[i]<='9')); ++i) { 00069 if (decimal_point) { 00070 if (time[i]=='.') { return -1; } 00071 samples = samples + ((mrs_real)(time[i]-'0'))/divisor; 00072 divisor = divisor * 10.0; 00073 } else if (time[i]=='.') { 00074 decimal_point=true; 00075 } else { 00076 samples = samples * 10.0 + (time[i]-'0'); 00077 } 00078 } 00079 // 00080 if (i<len) { 00081 char a=time[++i]; 00082 if (i>=len) { 00083 if (a=='h') { // hours 00084 samples= 120.0*samples*srate; 00085 } else if (a=='m') { // minutes 00086 samples= 60.0*samples*srate; 00087 } else if (a=='s') { // seconds 00088 samples= samples*srate; 00089 } else { 00090 return -1; 00091 } 00092 } else { 00093 char b=time[i]; 00094 if ((i+1)>=len) { 00095 if (a=='u' && b=='s') { // micro-seconds 00096 samples= samples/1000000.0*srate; 00097 } else if (a=='m' && b=='s') { // milli-seconds 00098 samples= samples/1000.0*srate; 00099 } else { 00100 return -1; 00101 } 00102 } 00103 } 00104 } 00105 return (mrs_natural)samples; 00106 } 00107 mrs_natural 00108 Marsyas::time2usecs(mrs_string time) { 00109 //example times: { "10us", "10ms", "10s", "10m", "10h" } 00110 if (time=="") { return 0; } 00111 // calculate time value 00112 mrs_real usecs=0; 00113 int i=0; 00114 int len=(int)time.length(); 00115 bool decimal_point=false; 00116 mrs_real divisor = 10.0; 00117 for (i=0; i<len && (time[i]=='.' || (time[i]>='0' && time[i]<='9')); ++i) { 00118 if (decimal_point) { 00119 if (time[i]=='.') { return -1; } 00120 usecs = usecs + ((mrs_real)(time[i]-'0'))/divisor; 00121 divisor = divisor * 10.0; 00122 } else if (time[i]=='.') { 00123 decimal_point=true; 00124 } else { 00125 usecs = usecs * 10.0 + (time[i]-'0'); 00126 } 00127 } 00128 // 00129 if (i<len) { 00130 char a=time[++i]; 00131 if (i>=len) { 00132 if (a=='h') { // hours 00133 usecs *= 1000.0 * 1000.0 * 60.0 * 60.0; 00134 } else if (a=='m') { // minutes 00135 usecs *= 1000.0 * 1000.0 * 60.0; 00136 } else if (a=='s') { // seconds 00137 usecs *= 1000.0 * 1000.0; 00138 } else { 00139 return -1; 00140 } 00141 } else { 00142 char b=time[i]; 00143 if ((i+1)>=len) { 00144 if (a=='u' && b=='s') { // micro-seconds 00145 ; 00146 } else if (a=='m' && b=='s') { // milli-seconds 00147 usecs *= 1000.0; 00148 } else { 00149 return -1; 00150 } 00151 } 00152 } 00153 } 00154 return (mrs_natural)usecs; 00155 } 00156 00157 mrs_real Marsyas::amplitude2dB(mrs_real a) 00158 { 00159 MRSASSERT (a > 0); 00160 return 20*log10(a); 00161 } 00162 00163 mrs_real Marsyas::dB2amplitude(mrs_real a) 00164 { 00165 return pow(10.0, a/20); 00166 } 00167 00168 mrs_real 00169 Marsyas::hertz2octs(mrs_real f, mrs_real middleAfreq) 00170 { 00171 //adapted from Dan Ellis fft2chromamx.m MATLAB routine 00172 // 00173 // octs = hz2octs(freq, A440) 00174 // Convert a frequency in Hz into a real number counting 00175 // the octaves above A0. So hz2octs(440) = 4.0 00176 // Optional A440 specifies the Hz to be treated as middle A (default 440). 00177 // 2006-06-29 dpwe@ee.columbia.edu for fft2chromamx 00178 // 00179 // A4 = 440 Hz, so A0 = 440/16 Hz 00180 // octs = log(freq./(A440/16))./log(2); 00181 // 00182 return log(f/(middleAfreq/16.0))/log(2.0); 00183 } 00184 00185 mrs_real Marsyas::hertz2bark(mrs_real f, mrs_natural mode) 00186 { 00187 switch (mode) 00188 { 00189 case 0: 00190 default: 00191 return 6 * log(f/600 + sqrt(1+ (pow(f/600,2)))); // 6*asinh(f/600); 00192 00193 case 1: // zwicker 00194 return 13 * atan (0.00076*f) + 3.5 * atan ((f*0.000133333333333333)*(f*0.000133333333333333)); 00195 00196 case 2: // terhardt 00197 return 13.3 * atan (0.00075*f); 00198 00199 case 3: // schroeder 00200 return 7.0 * log (f*0.00153846153846154 + sqrt ((f*0.00153846153846154)*(f*0.00153846153846154) + 1)); // 7*asinh (fFrequency/650); 00201 } 00202 } 00203 00204 mrs_real Marsyas::bark2hertz(mrs_real f, mrs_natural mode) 00205 { 00206 switch (mode) 00207 { 00208 case 0: 00209 case 1: 00210 default: 00211 return 600*sinh(f/6.0); 00212 00213 case 2: // terhardt 00214 return 4*1000.0/3.0 * tan (f*0.075187969924812); 00215 case 3: // schroeder 00216 return 650.0 * sinh (f*0.142857142857143); 00217 } 00218 } 00219 00220 mrs_real Marsyas::hertz2erb(mrs_real hz) 00221 { 00222 return 21.4 * log10(1+ (hz / 229.0)); 00223 } 00224 00225 mrs_real Marsyas::erb2hertz(mrs_real erb) 00226 { 00227 return (pow(10, (erb/21.4)) - 1.0) * 229.0; 00228 } 00229 00230 00231 mrs_real 00232 Marsyas::hertz2mel(mrs_real f, bool htk) 00233 { 00234 // z = hertz2mel(f,htk) 00235 // Convert frequencies f (in Hz) to mel 'scale'. 00236 // Optional htk = 1 uses the mel axis defined in the HTKBook 00237 // otherwise use Slaney's formula. 00238 // 00239 // adapted from Dan Ellis fft2melmx.m MATLAB code 00240 00241 if(htk) 00242 { 00243 return 2595.0 * log10(1.0 + f / 700.0); 00244 } 00245 else 00246 { 00247 // Mel fn to match Slaney's Auditory Toolbox mfcc.m 00248 mrs_real f_0 = 0.0; //133.33333; 00249 mrs_real f_sp = 200.0/3.0; //66.66667; 00250 mrs_real brkfrq = 1000.0; 00251 mrs_real brkpt = (brkfrq - f_0)/f_sp; //starting mel value for log region 00252 00253 //the magic 1.0711703 which is the ratio needed to get from 00254 //1000 Hz to 6400 Hz in 27 steps, and is *almost* the ratio between 00255 //1000 Hz and the preceding linear filter center at 933.33333 Hz 00256 //(actually 1000/933.33333 = 1.07142857142857 and 00257 //exp(log(6.4)/27) = 1.07117028749447) 00258 mrs_real logstep = exp(log(6.4)/27.0); 00259 00260 if(f < brkfrq) 00261 return (f - f_0) / f_sp; //linear 00262 else 00263 return brkpt + log(f / brkfrq) / log(logstep); //non-linear 00264 } 00265 } 00266 00267 mrs_real 00268 Marsyas::mel2hertz(mrs_real z, bool htk) 00269 { 00270 // f = mel2hz(z, htk) 00271 // Convert 'mel scale' frequencies into Hz 00272 // Optional htk = 1 means use the HTK formula 00273 // else use the formula from Slaney's mfcc.m 00274 // 00275 // Adapted from Dan Ellis fft2melmx.m MATLAB code 00276 00277 if(htk) 00278 { 00279 return 700.0 * (pow(10.0, z/2595.0) - 1.0); 00280 } 00281 else 00282 { 00283 mrs_real f_0 = 0.0; //133.33333; 00284 mrs_real f_sp = 200.0 / 3.0; //66.66667; 00285 mrs_real brkfrq = 1000.0; 00286 mrs_real brkpt = (brkfrq - f_0)/f_sp; //starting mel value for log region 00287 00288 //the magic 1.0711703 which is the ratio needed to get from 1000 Hz 00289 //to 6400 Hz in 27 steps, and is *almost* the ratio between 1000 Hz 00290 //and the preceding linear filter center at 933.33333 Hz 00291 //(actually 1000/933.33333 = 1.07142857142857 and 00292 //exp(log(6.4)/27) = 1.07117028749447) 00293 mrs_real logstep = exp(log(6.4)/27.0); 00294 00295 if(z < brkpt) 00296 return f_0 + f_sp * z; 00297 else 00298 return brkfrq * exp(log(logstep)*(z-brkpt)); 00299 } 00300 } 00301 00302 // Returns the frequency in Hz, represented by a PowerSpectrum bin number 00303 // (as returned from the PowerSpectrum Marsystem). 00304 // * bin is the bin number 00305 // * nyquistFreq is the Nyquist frequency (the sampling rate / 2) 00306 // * framerate is the samplerate / total # bins (usually we can pass in israte_) 00307 // 00308 // For example, if the nyquistFreq = 4000, the original sample rate was 8000Hz, 00309 // and we are using a framesize of 8 samples (8 FFT bins), then we should pass 00310 // in framerate=1000(Hz). Then there are 5 PowerSpect bins ((8/2)+1) 00311 // corresponding to frequencies as follows: 00312 // bin=0: returns 0 00313 // bin=1: returns 4000 (+/- 4000 Hz) 00314 // bin=2: returns 1000 (+/- 1000 Hz) 00315 // bin=3: returns 2000 (+/- 2000 Hz) 00316 // bin=4: returns 3000 (+/- 3000 Hz) 00317 mrs_real 00318 Marsyas::bin2hertz(mrs_natural bin, mrs_real nyquistFreq, mrs_real framerate) 00319 { 00320 if (bin == 0) 00321 return 0; 00322 if (bin == 1) 00323 return nyquistFreq; 00324 else 00325 return (bin-1) * framerate; 00326 } 00327 00328 mrs_natural Marsyas::powerOfTwo(mrs_real v) 00329 { 00330 mrs_natural n=1, res=0; 00331 while(res < v) 00332 { 00333 res = (mrs_natural) pow(2.0, n+.0); 00334 n++; 00335 } 00336 return res; 00337 } 00338 00339 void 00340 Marsyas::string2parameters(mrs_string s, realvec &v, char d) 00341 { 00342 mrs_natural i =0, pos=0, newPos=0; 00343 mrs_string tmp; 00344 while(newPos != -1 ) 00345 { 00346 newPos = (mrs_natural) s.find_first_of(&d, pos, 1); 00347 tmp = s.substr(pos, newPos); 00348 v(i++) = atof(tmp.c_str()); 00349 pos = newPos+1; 00350 } 00351 }