recursion;
+public:
+ FFT () {
+#ifdef USE_VDSP
+ m_fftSize = N;
+ m_fftSizeOver2 = m_fftSize / 2;
+ m_log2n = P; //log2f (m_fftSize); // bins
+ m_log2nOver2 = m_log2n / 2;
+
+ m_splitData.realp = new T[m_fftSize];
+ m_splitData.imagp = new T[m_fftSize];
+
+ m_scale = 1.0f / (T) (4.0f * m_fftSize);
+
+ // allocate the fft object once
+ m_fftSetup = vDSP_create_fftsetup (m_log2n, FFT_RADIX2);
+ if (m_fftSetup == NULL || /*m_inReal == NULL || m_outReal == NULL || */
+ m_splitData.realp == NULL || m_splitData.imagp == NULL /*|| m_window == NULL*/) {
+ throw std::runtime_error ("FFT_Setup failed to allocate enough memory");
+ }
+#endif
+ }
+
+ virtual ~FFT () {
+#ifdef USE_VDSP
+ delete [] (m_splitData.realp);
+ delete [] (m_splitData.imagp);
+#endif
+ }
+ void conjugate (T* data) {
+ for (int i = 0; i < N; ++i) {
+ data[2 * i + 1] *= -1;
+ }
+ }
+ void scramble (T* data) {
+ int i, m, j = 1;
+ for (i = 1; i < 2 * N; i += 2) {
+ if (j > i) {
+ std::swap (data[j - 1], data[i - 1]);
+ std::swap (data[j], data[i]);
+ }
+ m = N;
+ while (m >= 2 && j > m) {
+ j -= m;
+ m >>= 1;
+ }
+ j += m;
+ }
+ }
+ public:
+ enum { id = P };
+ void forward (T* data) {
+#ifdef USE_VDSP
+ //convert to split complex format with evens in real and odds in imag
+ vDSP_ctoz ((COMPLEX *) data, 2, &m_splitData, 1, m_fftSize);
+
+ //calc fft
+ vDSP_fft_zip (m_fftSetup, &m_splitData, 1, m_log2n, FFT_FORWARD);
+
+ //m_splitData.imagp[0] = 0.0;
+
+ vDSP_ztoc (&m_splitData, 1, (COMPLEX*) data, 2, m_fftSize);
+#else
+ scramble (data);
+ recursion.apply (data);
+#endif
+ }
+ void inverse (T* data) {
+#ifdef USE_VDSP
+ //convert to split complex format with evens in real and odds in imag
+ vDSP_ctoz ((COMPLEX *) data, 2, &m_splitData, 1, m_fftSize);
+
+ vDSP_fft_zip (m_fftSetup, &m_splitData, 1, m_log2n, FFT_INVERSE);
+ vDSP_ztoc (&m_splitData, 1, (COMPLEX*) data, 2, m_fftSize);
+
+ //vDSP_vsmul (buffer, 1, &m_scale, buffer, 1, m_fftSize);
+#else
+ conjugate (data);
+ scramble (data);
+ recursion.apply (data);
+
+#endif
+ }
+private:
+#ifdef USE_VDSP
+ int m_fftSize;
+ int m_fftSizeOver2;
+ int m_log2n;
+ int m_log2nOver2;
+ int m_windowSize;
+
+ T* m_inReal;
+ T* m_outReal;
+ T* m_window;
+
+ T m_scale;
+
+ FFTSetup m_fftSetup;
+ COMPLEX_SPLIT m_splitData;
+#endif
+};
+
+template
+ AbstractFFT* createFFT (int N) {
+ switch (N) {
+ case 4:
+ return new FFT<2, T> ();
+ break;
+ case 8:
+ return new FFT<3, T> ();
+ break;
+ case 16:
+ return new FFT<4, T> ();
+ break;
+ case 32:
+ return new FFT<5, T> ();
+ break;
+ case 64:
+ return new FFT<6, T> ();
+ break;
+ case 128:
+ return new FFT<7, T> ();
+ break;
+ case 256:
+ return new FFT<8, T> ();
+ break;
+ case 512:
+ return new FFT<9, T> ();
+ break;
+ case 1024:
+ return new FFT<10, T> ();
+ break;
+ case 2048:
+ return new FFT<11, T> ();
+ break;
+ case 4096:
+ return new FFT<12, T> ();
+ break;
+ case 8192:
+ return new FFT<13, T> ();
+ break;
+ case 16384:
+ return new FFT<14, T> ();
+ break;
+ case 32768:
+ return new FFT<15, T> ();
+ break;
+ case 65536:
+ return new FFT<16, T> ();
+ break;
+ case 131072:
+ return new FFT<17, T> ();
+ break;
+ case 262144:
+ return new FFT<18, T> ();
+ break;
+ case 524288:
+ return new FFT<19, T> ();
+ break;
+ case 1048576:
+ return new FFT<20, T> ();
+ break;
+ case 2097152:
+ return new FFT<21, T> ();
+ break;
+ case 4194304:
+ return new FFT<22, T> ();
+ break;
+ case 8388608:
+ return new FFT<23, T> ();
+ break;
+ case 16777216:
+ return new FFT<24, T> ();
+ break;
+ default:
+ throw std::runtime_error ("invalid size requested for fft");
+ }
+}
+
+template
+void fft (T *fftBuffer, long fftFrameSize, long sign) {
+ T wr, wi, arg, *p1, *p2, temp;
+ T tr, ti, ur, ui, *p1r, *p1i, *p2r, *p2i;
+ long i, bitm, j, le, le2, k;
+
+ int fftFrameSize2 = fftFrameSize << 1;
+
+ for (i = 2; i < fftFrameSize2 - 2; i += 2) {
+ for (bitm = 2, j = 0; bitm < fftFrameSize2; bitm <<= 1) {
+ if (i & bitm) j++;
+ j <<= 1;
+ }
+ if (i < j) {
+ p1 = fftBuffer + i;
+ p2 = fftBuffer + j;
+ temp = *p1;
+ *(p1++) = *p2;
+ *(p2++) = temp;
+ temp = *p1;
+ *p1 = *p2;
+ *p2 = temp;
+ }
+ }
+
+ for (k = 0, le = 2; k < ceil (log2f (fftFrameSize)); k++) {
+ le <<= 1;
+ le2 = le >> 1;
+ ur = 1.;
+ ui = 0.;
+ arg = PI / (le2 >> 1);
+ wr = cosf (arg);
+ wi = sign* sinf (arg);
+ for (j = 0; j < le2; j += 2) {
+ p1r = fftBuffer + j;
+ p1i = p1r + 1;
+ p2r = p1r + le2;
+ p2i = p2r + 1;
+ for (i = j; i < fftFrameSize2; i += le) {
+ tr = *p2r * ur - *p2i * ui;
+ ti = *p2r * ui + *p2i * ur;
+ *p2r = *p1r - tr;
+ *p2i = *p1i - ti;
+ *p1r += tr;
+ *p1i += ti;
+ p1r += le;
+ p1i += le;
+ p2r += le;
+ p2i += le;
+ }
+ tr = ur*wr - ui*wi;
+ ui = ur*wi + ui*wr;
+ ur = tr;
+ }
+ }
+}
+
+// ---------------------------------------------------------------------- //
+template
+void hanningz (T* out, int N) {
+ for (int i = 0; i < N; ++i) {
+ out[i] = .5 * (1. - cos (TWOPI * (T) i / N));
+ }
+}
+
+template
+void makeWindow (T* out, int N, T a0, T a1, T a2) {
+ // .5, .5, 0 --> hanning
+ // .54, .46, 0 --> hamming
+ // .42, .5, 0.08 --> blackman
+ for (int i = 0; i < N; ++i) {
+ out[i] = a0 - a1 * cos ((TWOPI * (T) i) / (N - 1)) + a2 *
+ cos ((2 * TWOPI * (T) i) / (N - 1)); // hamming, hann or blackman
+ }
+}
+
+template
+void fftshift (T* cbuf, int N) {
+ int N2 = N >> 1;
+ for (int i = 0; i < (N2); ++i) {
+ //cout << "\t" << i << endl;
+ T tmp = cbuf[i];
+ cbuf[i] = cbuf[i + N2];
+ cbuf[i + N2] = tmp;
+ }
+}
+
+template
+T princarg (T phi) {
+ T a = phi / TWOPI;
+ int k = (int) round (a);
+ return phi - (T) k * TWOPI;
+}
+
+template
+void pol2rect (T* cbuf, int N) {
+#ifdef USE_VDSP
+ vDSP_rect (cbuf, 2, cbuf, 2, N);
+#else
+ for (int i = 0; i < N; ++i) {
+ T amp = cbuf[2 * i];
+ T phi = cbuf[2 * i + 1];
+ T real = amp * cosf (phi);
+ T imag = amp * sinf (phi);
+ cbuf[2 * i] = real;
+ cbuf[2 * i + 1] = imag;
+ }
+#endif
+}
+
+template
+void rect2pol (T* cbuf, int N) {
+#ifdef USE_VDSP
+ vDSP_polar (cbuf, 2, cbuf, 2, N);
+#else
+ for (int i = 0; i < N; ++i) {
+ T real = cbuf[2 * i];
+ T imag = cbuf[2 * i + 1];
+ T amp = sqrtf (real * real + imag * imag);
+ //T amp = hypot (real, imag);
+ T phi = atan2f (imag, real);
+ cbuf[2 * i] = amp;
+ cbuf[2 * i + 1] = phi;
+ }
+#endif
+}
+
+template
+int locmax (const T* amp, int N, int* max) {
+ T maxPeak = amp[1];
+ int count = 0;
+ for (int i = 1; i < N - 1; ++i) {
+ T magCurr = amp[i];
+ T magPrev = amp[(i - 1)];
+ T magNext = amp[(i + 1)];
+
+ if (magCurr > magPrev && magCurr > magNext) {
+ max[count] = i;
+ ++count;
+ if (magCurr > maxPeak) maxPeak = magCurr;
+ }
+ }
+ return count;
+}
+
+template
+ int locmax2 (const T* cplx, int N, int* max) {
+ T maxPeak = cplx[2];
+ int count = 0;
+ for (int i = 2; i < N - 2; ++i) {
+ T magCurr = cplx[2 * i];
+ T magPrev = cplx[2 * (i - 1)];
+ T magPrevPrev = cplx[2 * (i - 2)];
+ T magNext = cplx[2 * (i + 1)];
+ T magNextNext = cplx[2 * (i + 2)];
+
+ if (magCurr > magPrev && magCurr > magNext &&
+ magCurr > magPrevPrev && magCurr > magNextNext) {
+ max[count] = i;
+ ++count;
+ if (magCurr > maxPeak) maxPeak = magCurr;
+ }
+ }
+ return count;
+}
+template
+T locmax2AmpFreq (const T* amp, const T* freq, int size, std::vector& max, T FUNDAMENTAL) {
+ T maxPeak = amp[2];
+ for (int i = 2; i < size - 2; ++i) {
+ T magCurr = amp[i]; //20 * log10 (amp[i] + .00000001);
+ T magPrev = amp[i - 1]; //20 * log10 (amp[i - 1] + .00000001);
+ T magPrevPrev = amp[i - 2]; //20 * log10 (amp[i - 1] + .00000001);
+ T magNext = amp[i + 1]; // 20 * log10 (amp[i + 1] + .00000001);
+ T magNextNext = amp[i + 2]; // 20 * log10 (amp[i + 1] + .00000001);
+
+ if (magCurr > magPrev && magCurr > magNext &&
+ magCurr > magPrevPrev && magCurr > magNextNext) {
+ max.push_back (i);
+ if (magCurr > maxPeak) maxPeak = magCurr;
+ }
+ }
+
+ return maxPeak;
+}
+
+// rect to amp-freq
+template
+inline void ampFreqPhaseDiff (const T* cbuffer, T* amp, T* freq, T* oldPhi,
+ int N, int hop, T R) {
+ T osamp = N / hop;
+ T freqPerBin = R / (T) N;
+ T expct = TWOPI * (T) hop / (T) N;
+ rect2pol (cbuffer, N);
+ for (int i = 0; i < N; ++i) {
+ amp[i] = cbuffer[2 * i];
+ T phase = cbuffer[2 * i + 1];
+ T tmp = phase - oldPhi[i];
+ oldPhi[i] = phase;
+ tmp -= (T) i * expct;
+ int qpd = (int) (tmp / PI);
+ if (qpd >= 0) qpd += qpd & 1;
+ else qpd -= qpd & 1;
+ tmp -= PI * (T) qpd;
+ tmp = osamp * tmp / (2. * PI);
+ tmp = (T) i * freqPerBin + tmp * freqPerBin;
+ freq[i] = tmp;
+ }
+}
+
+template
+void ampFreqParabolic (const T* cbuffer, T* amp, T* freq, int N, double R) {
+ T freqPerBin = (R ) / (T) N;
+ T min = 0;
+ for (int i = 0; i < N; ++i) {
+ amp[i] = sqrtf (cbuffer[2 * i] * cbuffer[2 * i] + cbuffer[2 * i + 1] * cbuffer[2 * i + 1]);
+ }
+
+ freq[0] = freqPerBin;
+ for (int i = 1; i < N - 1; ++i) {
+ freq[i] = parabolicInterpolate (freqPerBin * (i - 1), freqPerBin * (i), freqPerBin * (i + 1),
+ amp[i - 1], amp[i], amp[i + 1], &min);
+ }
+}
+
+template
+void ampFreqBins (const T* cbuffer, T* amp, T* freq, int N, double R) {
+ T freqPerBin = (R) / (T) N;
+ for (int i = 0; i < N; ++i) {
+ amp[i] = sqrtf (cbuffer[2 * i] * cbuffer[2 * i] + cbuffer[2 * i + 1] * cbuffer[2 * i + 1]);
+ }
+
+ for (int i = 1; i < N - 1; ++i) {
+ freq[i] = (T) i * freqPerBin;
+ }
+}
+template
+void sortSpectrum (Peak* peaks, int k) {
+ Peak p;
+ for (int i = 0; i < k; ++i) {
+ for (int j = i + 1; j < k; ++j) {
+ //if (peaks[i].freq > peaks[j].freq) {
+ if (peaks[i].amp < peaks[j].amp) {
+ p = peaks[i];
+ peaks[i] = peaks[j];
+ peaks[j] = p;
+ }
+ }
+ }
+}
+
+#endif // FFT_H
+
+// EOF
diff --git a/source/cella/WavFile.h b/source/cella/WavFile.h
new file mode 100644
index 0000000..2f6f86e
--- /dev/null
+++ b/source/cella/WavFile.h
@@ -0,0 +1,667 @@
+// WavFile.h
+//
+
+
+#ifndef WAVFILE_H
+#define WAVFILE_H
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#ifndef uint
+ typedef unsigned int uint;
+#endif
+
+const static char riffStr[] = "RIFF";
+const static char waveStr[] = "WAVE";
+const static char fmtStr[] = "fmt ";
+const static char dataStr[] = "data";
+
+typedef struct {
+ char riff_char[4];
+ int package_len;
+ char wave[4];
+} WavRiff;
+
+typedef struct {
+ char fmt[4];
+ int format_len;
+ short fixed;
+ short channel_number;
+ int sample_rate;
+ int byte_rate;
+ short byte_per_sample;
+ short bits_per_sample;
+} WavFormat;
+
+typedef struct {
+ char data_field[4];
+ uint data_len;
+} WavData;
+
+
+typedef struct {
+ WavRiff riff;
+ WavFormat format;
+ WavData data;
+} WavHeader;
+
+
+#ifdef BYTE_ORDER
+ #if BYTE_ORDER == BIG_ENDIAN
+ #define _BIG_ENDIAN_
+ #endif
+#endif
+
+#ifdef _BIG_ENDIAN_
+ static inline void _swap32 (unsigned int &dwData) {
+ dwData = ((dwData >> 24) & 0x000000FF) |
+ ((dwData >> 8) & 0x0000FF00) |
+ ((dwData << 8) & 0x00FF0000) |
+ ((dwData << 24) & 0xFF000000);
+ }
+
+ static inline void _swap16 (unsigned short &wData) {
+ wData = ((wData >> 8) & 0x00FF) |
+ ((wData << 8) & 0xFF00);
+ }
+
+ static inline void _swap16Buffer (unsigned short *pData, unsigned int dwNumWords) {
+ unsigned long i;
+
+ for (i = 0; i < dwNumWords; i ++) {
+ _swap16 (pData[i]);
+ }
+ }
+
+#else // BIG_ENDIAN
+ static inline void _swap32 (unsigned int &dwData) {
+ // do nothing
+ }
+
+ static inline void _swap16 (unsigned short &wData) {
+ // do nothing
+ }
+
+
+ static inline void _swap16Buffer (unsigned short *pData, unsigned int dwNumBytes) {
+ // do nothing
+ }
+
+#endif // BIG_ENDIAN
+
+// test if character code is between a white space ' ' and little 'z'
+static int isAlpha (char c) {
+ return (c >= ' ' && c <= 'z') ? 1 : 0;
+}
+
+
+// test if all characters are between a white space ' ' and little 'z'
+static int isAlphaStr (char *str) {
+ int c;
+
+ c = str[0];
+ while (c) {
+ if (isAlpha(c) == 0) return 0;
+ str ++;
+ c = str[0];
+ }
+
+ return 1;
+}
+
+class WavInFile {
+public:
+ WavInFile (const char *fileName) {
+ int hdrsOk;
+
+ // Try to open the file for reading
+ fptr = fopen(fileName, "rb");
+ if (fptr == NULL) {
+ // didn't succeed
+ std::string msg = "unable to open file \"";
+ msg += fileName;
+ msg += "\" for reading.";
+ throw std::runtime_error(msg);
+ }
+
+ // Read the file headers
+ hdrsOk = readWavHeaders();
+ if (hdrsOk != 0) {
+ // Something didn't match in the wav file headers
+ std::string msg = "File \"";
+ msg += fileName;
+ msg += "\" is corrupt or not a WAV file";
+ throw std::runtime_error(msg);
+ }
+
+ if (header.format.fixed != 1) {
+ std::string msg = "File \"";
+ msg += fileName;
+ msg += "\" uses unsupported encoding.";
+ throw std::runtime_error(msg);
+ }
+
+ dataRead = 0;
+ }
+
+
+
+ ~WavInFile() {
+ close();
+ }
+
+
+
+ void rewind() {
+ int hdrsOk;
+
+ fseek(fptr, 0, SEEK_SET);
+ hdrsOk = readWavHeaders();
+ assert(hdrsOk == 0);
+ dataRead = 0;
+ }
+
+
+ int checkCharTags() {
+ // header.format.fmt should equal to 'fmt '
+ if (memcmp(fmtStr, header.format.fmt, 4) != 0) return -1;
+ // header.data.data_field should equal to 'data'
+ if (memcmp(dataStr, header.data.data_field, 4) != 0) return -1;
+
+ return 0;
+ }
+
+
+ int read(char *buffer, int maxElems) {
+ int numBytes;
+ uint afterDataRead;
+
+ // ensure it's 8 bit format
+ if (header.format.bits_per_sample != 8) {
+ throw std::runtime_error("read(char*, int) works only with 8bit samples.");
+ }
+ assert(sizeof(char) == 1);
+
+ numBytes = maxElems;
+ afterDataRead = dataRead + numBytes;
+ if (afterDataRead > header.data.data_len) {
+ // Don't read more samples than are marked available in header
+ numBytes = header.data.data_len - dataRead;
+ assert(numBytes >= 0);
+ }
+
+ numBytes = fread(buffer, 1, numBytes, fptr);
+ dataRead += numBytes;
+
+ return numBytes;
+ }
+
+
+ int read(short *buffer, int maxElems) {
+ unsigned int afterDataRead;
+ int numBytes;
+ int numElems;
+
+ if (header.format.bits_per_sample == 8) {
+ // 8 bit format
+ char *temp = new char[maxElems];
+ int i;
+
+ numElems = read(temp, maxElems);
+ // convert from 8 to 16 bit
+ for (i = 0; i < numElems; i ++) {
+ buffer[i] = temp[i] << 8;
+ }
+ delete[] temp;
+ } else {
+ // 16 bit format
+ assert(header.format.bits_per_sample == 16);
+ assert(sizeof(short) == 2);
+
+ numBytes = maxElems * 2;
+ afterDataRead = dataRead + numBytes;
+ if (afterDataRead > header.data.data_len) {
+ // Don't read more samples than are marked available in header
+ numBytes = header.data.data_len - dataRead;
+ assert(numBytes >= 0);
+ }
+
+ numBytes = fread(buffer, 1, numBytes, fptr);
+ dataRead += numBytes;
+ numElems = numBytes / 2;
+
+ // 16bit samples, swap byte order if necessary
+ _swap16Buffer((unsigned short *)buffer, numElems);
+ }
+
+ return numElems;
+ }
+
+
+
+ int read(float *buffer, int maxElems) {
+ short *temp = new short[maxElems];
+ int num;
+ int i;
+ double fscale;
+
+ num = read(temp, maxElems);
+
+ fscale = 1.0 / 32768.0;
+ // convert to floats, scale to range [-1..+1[
+ for (i = 0; i < num; i ++) {
+ buffer[i] = (float)(fscale * (double)temp[i]);
+ }
+
+ delete[] temp;
+
+ return num;
+ }
+
+ int read(double *buffer, int maxElems) {
+ short *temp = new short[maxElems];
+ int num;
+ int i;
+ double fscale;
+
+ num = read(temp, maxElems);
+
+ fscale = 1.0 / 32768.0;
+ // convert to doubles, scale to range [-1..+1[
+ for (i = 0; i < num; i ++) {
+ buffer[i] = (double)(fscale * (double)temp[i]);
+ }
+
+ delete[] temp;
+
+ return num;
+ }
+
+ int eof() const {
+ // return true if all data has been read or file eof has reached
+ return (dataRead == header.data.data_len || feof(fptr));
+ }
+
+
+ void close() {
+ fclose(fptr);
+ fptr = NULL;
+ }
+
+
+ int readRIFFBlock() {
+ fread(&(header.riff), sizeof(WavRiff), 1, fptr);
+
+ // swap 32bit data byte order if necessary
+ _swap32((unsigned int &)header.riff.package_len);
+
+ // header.riff.riff_char should equal to 'RIFF');
+ if (memcmp(riffStr, header.riff.riff_char, 4) != 0) return -1;
+ // header.riff.wave should equal to 'WAVE'
+ if (memcmp(waveStr, header.riff.wave, 4) != 0) return -1;
+
+ return 0;
+ }
+
+ int readHeaderBlock() {
+ char label[5];
+ std::string sLabel;
+
+ // lead label std::string
+ fread(label, 1, 4, fptr);
+ label[4] = 0;
+
+ if (isAlphaStr(label) == 0) return -1; // not a valid label
+
+ // Decode blocks according to their label
+ if (strcmp(label, fmtStr) == 0) {
+ int nLen, nDump;
+
+ // 'fmt ' block
+ memcpy(header.format.fmt, fmtStr, 4);
+
+ // read length of the format field
+ fread(&nLen, sizeof(int), 1, fptr);
+ // swap byte order if necessary
+ _swap32((unsigned int &)nLen); // int format_len;
+ header.format.format_len = nLen;
+
+ // calculate how much length differs from expected
+ nDump = nLen - (sizeof(header.format) - 8);
+
+ // if format_len is larger than expected, read only as much data as we've space for
+ if (nDump > 0) {
+ nLen = sizeof(header.format) - 8;
+ }
+
+ // read data
+ fread(&(header.format.fixed), nLen, 1, fptr);
+
+ // swap byte order if necessary
+ _swap16((unsigned short &)header.format.fixed); // short int fixed;
+ _swap16((unsigned short &)header.format.channel_number); // short int channel_number;
+ _swap32((unsigned int &)header.format.sample_rate); // int sample_rate;
+ _swap32((unsigned int &)header.format.byte_rate); // int byte_rate;
+ _swap16((unsigned short &)header.format.byte_per_sample); // short int byte_per_sample;
+ _swap16((unsigned short &)header.format.bits_per_sample); // short int bits_per_sample;
+
+ // if format_len is larger than expected, skip the extra data
+ if (nDump > 0) {
+ fseek(fptr, nDump, SEEK_CUR);
+ }
+
+ return 0;
+ } else if (strcmp(label, dataStr) == 0) {
+ // 'data' block
+ memcpy(header.data.data_field, dataStr, 4);
+ fread(&(header.data.data_len), sizeof(uint), 1, fptr);
+
+ // swap byte order if necessary
+ _swap32((unsigned int &)header.data.data_len);
+
+ return 1;
+ } else {
+ uint len, i;
+ uint temp;
+ // unknown block
+
+ // read length
+ fread(&len, sizeof(len), 1, fptr);
+ // scan through the block
+ for (i = 0; i < len; i ++) {
+ fread(&temp, 1, 1, fptr);
+ if (feof(fptr)) return -1; // unexpected eof
+ }
+ }
+ return 0;
+ }
+
+
+ int readWavHeaders() {
+ int res;
+
+ memset(&header, 0, sizeof(header));
+
+ res = readRIFFBlock();
+ if (res) return 1;
+ // read header blocks until data block is found
+ do {
+ // read header blocks
+ res = readHeaderBlock();
+ if (res < 0) return 1; // error in file structure
+ } while (res == 0);
+ // check that all required tags are legal
+ return checkCharTags();
+ }
+
+
+ uint getNumChannels() const {
+ return header.format.channel_number;
+ }
+
+
+ uint getNumBits() const {
+ return header.format.bits_per_sample;
+ }
+
+
+ uint getBytesPerSample() const {
+ return getNumChannels() * getNumBits() / 8;
+ }
+
+
+ uint getSampleRate() const {
+ return header.format.sample_rate;
+ }
+
+
+
+ uint getDataSizeInBytes() const {
+ return header.data.data_len;
+ }
+
+
+ uint getNumSamples() const {
+ return header.data.data_len / header.format.byte_per_sample;
+ }
+
+
+ uint getLengthMS() const {
+ uint numSamples;
+ uint sampleRate;
+
+ numSamples = getNumSamples();
+ sampleRate = getSampleRate();
+
+ assert(numSamples < UINT_MAX / 1000);
+ return (1000 * numSamples / sampleRate);
+ }
+private:
+ FILE *fptr;
+ uint dataRead;
+ WavHeader header;
+};
+
+
+/// Class for writing WAV audio files.
+class WavOutFile {
+private:
+ /// Pointer to the WAV file
+ FILE *fptr;
+
+ /// WAV file header data.
+ WavHeader header;
+
+ /// Counter of how many bytes have been written to the file so far.
+ int bytesWritten;
+
+public:
+ //////////////////////////////////////////////////////////////////////////////
+ //
+ // Class WavOutFile
+ //
+
+ WavOutFile(const char *fileName, int sampleRate, int bits, int channels) {
+ bytesWritten = 0;
+ fptr = fopen(fileName, "wb");
+ if (fptr == NULL) {
+ std::string msg = "unable to open file \"";
+ msg += fileName;
+ msg += "\" for writing.";
+ //pmsg = msg.c_str;
+ throw std::runtime_error(msg);
+ }
+
+ fillInHeader(sampleRate, bits, channels);
+ writeHeader();
+ }
+
+
+
+ ~WavOutFile() {
+ close();
+ }
+
+
+
+ void fillInHeader(uint sampleRate, uint bits, uint channels) {
+ // fill in the 'riff' part..
+
+ // copy std::string 'RIFF' to riff_char
+ memcpy(&(header.riff.riff_char), riffStr, 4);
+ // package_len unknown so far
+ header.riff.package_len = 0;
+ // copy std::string 'WAVE' to wave
+ memcpy(&(header.riff.wave), waveStr, 4);
+
+
+ // fill in the 'format' part..
+
+ // copy std::string 'fmt ' to fmt
+ memcpy(&(header.format.fmt), fmtStr, 4);
+
+ header.format.format_len = 0x10;
+ header.format.fixed = 1;
+ header.format.channel_number = (short)channels;
+ header.format.sample_rate = sampleRate;
+ header.format.bits_per_sample = (short)bits;
+ header.format.byte_per_sample = (short)(bits * channels / 8);
+ header.format.byte_rate = header.format.byte_per_sample * sampleRate;
+ header.format.sample_rate = sampleRate;
+
+ // fill in the 'data' part..
+
+ // copy std::string 'data' to data_field
+ memcpy(&(header.data.data_field), dataStr, 4);
+ // data_len unknown so far
+ header.data.data_len = 0;
+ }
+
+
+ void finishHeader() {
+ // supplement the file length into the header structure
+ header.riff.package_len = bytesWritten + 36;
+ header.data.data_len = bytesWritten;
+
+ writeHeader();
+ }
+
+
+
+ void writeHeader() {
+ WavHeader hdrTemp;
+
+ // swap byte order if necessary
+ hdrTemp = header;
+ _swap32((unsigned int &)hdrTemp.riff.package_len);
+ _swap32((unsigned int &)hdrTemp.format.format_len);
+ _swap16((unsigned short &)hdrTemp.format.fixed);
+ _swap16((unsigned short &)hdrTemp.format.channel_number);
+ _swap32((unsigned int &)hdrTemp.format.sample_rate);
+ _swap32((unsigned int &)hdrTemp.format.byte_rate);
+ _swap16((unsigned short &)hdrTemp.format.byte_per_sample);
+ _swap16((unsigned short &)hdrTemp.format.bits_per_sample);
+ _swap32((unsigned int &)hdrTemp.data.data_len);
+
+ // write the supplemented header in the beginning of the file
+ fseek(fptr, 0, SEEK_SET);
+ fwrite(&hdrTemp, sizeof(hdrTemp), 1, fptr);
+ // jump back to the end of the file
+ fseek(fptr, 0, SEEK_END);
+ }
+
+
+
+ void close() {
+ finishHeader();
+ fclose(fptr);
+ fptr = NULL;
+ }
+
+
+ void write(const char *buffer, int numElems) {
+ int res;
+
+ if (header.format.bits_per_sample != 8) {
+ throw std::runtime_error("write(const char*, int) accepts only 8bit samples.");
+ }
+ assert(sizeof(char) == 1);
+
+ res = fwrite(buffer, 1, numElems, fptr);
+ if (res != numElems) {
+ throw std::runtime_error("problem while writing to a wav file.");
+ }
+
+ bytesWritten += numElems;
+ }
+
+
+ void write(const short *buffer, int numElems) {
+ int res;
+
+ // 16 bit samples
+ if (numElems < 1) return; // nothing to do
+
+ if (header.format.bits_per_sample == 8) {
+ int i;
+ char *temp = new char[numElems];
+ // convert from 16bit format to 8bit format
+ for (i = 0; i < numElems; i ++) {
+ temp[i] = buffer[i] >> 8;
+ }
+ // write in 8bit format
+ write(temp, numElems);
+ delete[] temp;
+ } else {
+ // 16bit format
+ unsigned short *pTemp = new unsigned short[numElems];
+
+ assert(header.format.bits_per_sample == 16);
+
+ // allocate temp buffer to swap byte order if necessary
+ memcpy(pTemp, buffer, numElems * 2);
+ _swap16Buffer(pTemp, numElems);
+
+ res = fwrite(pTemp, 2, numElems, fptr);
+
+ delete[] pTemp;
+
+ if (res != numElems) {
+ throw std::runtime_error("problem while writing to a wav file.");
+ }
+ bytesWritten += 2 * numElems;
+ }
+ }
+
+
+ void write(const float *buffer, int numElems) {
+ int i;
+ short *temp = new short[numElems];
+ int iTemp;
+
+ // convert to 16 bit integer
+ for (i = 0; i < numElems; i ++) {
+ // convert to integer
+ iTemp = (int)(32768.0f * buffer[i]);
+
+ // saturate
+ if (iTemp < -32768) iTemp = -32768;
+ if (iTemp > 32767) iTemp = 32767;
+ temp[i] = (short)iTemp;
+ }
+
+ write(temp, numElems);
+
+ delete[] temp;
+ }
+
+ void write(const double *buffer, int numElems) {
+ int i;
+ short *temp = new short[numElems];
+ int iTemp;
+
+ // convert to 16 bit integer
+ for (i = 0; i < numElems; i ++) {
+ // convert to integer
+ iTemp = (int)(32768.0f * buffer[i]);
+
+ // saturate
+ if (iTemp < -32768) iTemp = -32768;
+ if (iTemp > 32767) iTemp = 32767;
+ temp[i] = (short)iTemp;
+ }
+
+ write(temp, numElems);
+
+ delete[] temp;
+ }
+
+};
+
+#endif
diff --git a/source/cella/algorithms.h b/source/cella/algorithms.h
new file mode 100644
index 0000000..05621e7
--- /dev/null
+++ b/source/cella/algorithms.h
@@ -0,0 +1,299 @@
+// algorithms.h
+//
+
+#ifndef ALGORITHMS_H
+#define ALGORITHMS_H
+
+#include
+#include
+
+#define P0 (1 / (std::pow (2.0, 32.) / 2 - 1))
+#define LOG2 log (2.)
+#define LOG2OF10 3.32192809488736234787
+
+// numerical/statistic -------------------------------------------------- //
+template
+inline T logTwo (const T& x) {
+ return LOG2OF10 * std::log10 ((T) x);
+}
+
+template
+void normalize (const T* data, T* result, int N) {
+ T min = data[0];
+ T max = data[0];
+ T mean = 0;
+
+ for (int i = 0; i < N; ++i) {
+ if (max < data[i]) max = data[i];
+ if (min > data[i]) min = data[i];
+ }
+
+ T fmin = fabs (min);
+ for (int i = 0; i < N; ++i) {
+ result[i] = (data[i] + fmin) / (max + fmin);
+ }
+}
+
+template
+void scale (const T* buff, T* out, int len, T factor) {
+ for (int i = 0; i < len; ++i) {
+ out[i] = buff[i] * factor;
+ }
+}
+
+template
+void swapElem (T& a, T& b) {
+ T t = (a); (a) = (b); (b) = t;
+}
+
+template
+T kth_smallest (T a[], int n, int k) {
+ int i, j, l, m;
+ T x ;
+ l = 0 ; m = n - 1 ;
+ while (l < m) {
+ x = a[k];
+ i = l;
+ j= m;
+ do {
+ while (a[i] < x) i++;
+ while (x < a[j]) j--;
+ if (i <= j) {
+ swapElem (a[i], a[j]);
+ i++; j--;
+ }
+ } while (i <= j) ;
+ if (j < k) l = i;
+ if (k < i) m = j;
+ }
+ return a[k] ;
+}
+
+template
+T median (T a[], int n) {
+ return kth_smallest (a, n, (((n) & 1) ? ((n) / 2):(((n) / 2) -1)));
+}
+
+template
+T sum (
+ const T* values,
+ int N) {
+ if (1 > N) return 0;
+
+ T sum = 0.;
+ for (int i = 0; i < N; ++i) {
+ sum += values[i];
+ }
+
+ return sum;
+}
+
+template
+T mean (
+ const T* values,
+ int N) {
+ if (1 > N) return 0;
+
+ return sum (values, N) / N;
+}
+
+template
+T min (
+ const T* values,
+ int N) {
+ if (1 > N) return 0;
+
+ T m = values[0];
+ for (int i = 0; i < N; ++i) {
+ if (m > values[i]) {
+ m = values[i];
+ }
+ }
+
+ return m;
+}
+
+template
+T max (
+ const T* values,
+ int N) {
+ if (1 > N) return 0;
+
+ T m = values[0];
+ for (int i = 0; i < N; ++i) {
+ if (m < values[i]) {
+ m = values[i];
+ }
+ }
+
+ return m;
+}
+
+template
+T centroid (
+ const T* weights,
+ const T* values,
+ int N) {
+ T sumWeigth = 0;
+ T sumWeigthDistance = 0;
+ T distance = 0;
+ T weigth = 0;
+ T weigthDistance = 0;
+
+ T maxAmp = weights[0];
+ for (int index = 0; index < N; ++index) {
+ distance = values[index];
+ weigth = weights[index] / N;
+ weigthDistance = weigth * distance;
+ sumWeigth += weigth;
+ sumWeigthDistance += weigthDistance;
+ if (maxAmp < weights[index]) {
+ maxAmp = weights[index];
+ }
+ }
+ return (sumWeigth > .1 ? sumWeigthDistance / sumWeigth : 0);
+}
+
+template
+T stddev (const T* weights, const T* values, T mean, int N) {
+ if (N <= 0) return 0;
+
+ T max = weights[0];
+ for (int f = 0; f < N; ++f) {
+ if (weights[f] >= max) max = weights[f];
+ }
+ T scaledMax = 0.3 * max;
+
+ T sumWeight = 0.0;
+ T sumDistance = 0.0;
+ for (int index = 0; index < N; index++) {
+ T tmp = values[index] - mean;
+ T value = tmp * tmp;
+ T weight = weights[index];
+ if (weight <= scaledMax) weight = 0;
+ T distance = value * weight;
+
+ sumWeight += weight;
+ sumDistance += distance;
+ }
+
+ T sd = 0.0;
+ if (0 != sumWeight) {
+ sd = sumDistance / sumWeight;
+ } else {
+ sd = 0.;
+
+ }
+ return std::sqrt (sd);
+}
+
+template
+T moment (
+ const T* weights,
+ const T* values,
+ int N,
+ int order,
+ T centroid) {
+ T sumWeigth = 0;
+ T sumWeigthDistance = 0;
+ T distance = 0;
+ T weigth = 0;
+ T weigthDistance = 0;
+
+ for (int index = 0; index < N; ++index) {
+ distance = values[index];
+ distance -= centroid;
+ distance = std::pow (distance, order);
+ weigth = weights[index];
+ weigthDistance = weigth * distance;
+ sumWeigth += weigth;
+ sumWeigthDistance += weigthDistance;
+ }
+
+ return (sumWeigth != 0 ? sumWeigthDistance / sumWeigth : 0);
+}
+
+template
+T linreg (
+ const T* weights,
+ const T* values,
+ int N,
+ T& step) {
+ int index;
+ int nbValue = N;
+
+ T Xk = 0;
+ T Yk = 0;
+ T Xk2 = 0;
+ T XkYk = 0;
+
+ T sumXk = 0;
+ T sumYk = 0;
+ T sumXk2 = 0;
+ T sumXkYk = 0;
+
+ T slope = 0;
+
+ for (index = 0; index < nbValue; index++) {
+ Xk = values[index];
+ Yk = weights[index];
+ Xk2 = Xk * Xk;
+ XkYk = Xk * Yk;
+
+ sumXk += Xk;
+ sumYk += Yk;
+ sumXkYk += XkYk;
+ sumXk2 += Xk2;
+ }
+
+ T numSlope = ((T) nbValue) * sumXkYk - (sumXk * sumYk);
+ T numStep = sumXk2 * sumYk - sumXk * sumXkYk;
+ T denSlope = ((T) nbValue) * sumXk2 - std::pow (sumXk, 2.0);
+
+ if (0 != denSlope) {
+ slope = numSlope / denSlope;
+ step = numStep / denSlope;
+ } else {
+ slope = 0;
+ step = 0;
+ }
+ return slope;
+}
+
+// dsp ------------------------------------------------------------------ //
+
+
+template
+void conv (T* X, T* Y, T* Z, int lenx, int leny, T scale) {
+ T *zptr, s, *xp, *yp;
+ int i, n, n_lo, n_hi;
+ int lenz = lenx + leny - 1;
+
+ zptr = Z;
+ for (i = 0; i < lenz; i++) {
+ s=0.0;
+ n_lo= 0 > (i - leny + 1) ? 0 : i - leny + 1;
+ n_hi = lenx - 1 < i ? lenx - 1 : i;
+ xp = X + n_lo;
+ yp = Y + i-n_lo;
+ for (n = n_lo; n <= n_hi; n++) {
+ s += *xp * *yp;
+ xp++;
+ yp--;
+ }
+ *zptr = s * scale;
+ zptr++;
+ }
+}
+
+template
+T amp2db (const T& input) {
+ // sound pressure level (minimum 32 bits)
+ if (input < P0) return 0;
+ T f = 20. * std::log10 ((T) (input / P0));
+ return f;
+}
+
+#endif // ALGORITHMS_H
+
+// EOF
diff --git a/source/cella/analysis.h b/source/cella/analysis.h
new file mode 100644
index 0000000..93c4cf1
--- /dev/null
+++ b/source/cella/analysis.h
@@ -0,0 +1,112 @@
+// analysis.h
+//
+
+#ifndef ANALYSIS_H
+#define ANALYSIS_H
+
+#include "WavFile.h"
+#include "utilities.h"
+#include "parser.h"
+#include "FFT.h"
+#include "features.h"
+
+#include
+#include
+#include
+#include
+
+struct Descriptors {
+ std::vector highest_peak_Hz;
+ std::vector irr;
+ std::vector energy; //implemented in features.h
+ std::vector zcr; //implemented in features.h
+};
+
+void analyse (const std::string& fname, Parameters& p, Descriptors& d) {
+ WavInFile in (fname.c_str ());
+ if (in.getSampleRate () != p.sr) {
+ std::stringstream err;
+ err << "invalid sampling rate in " << fname << ": " << in.getSampleRate () << " vs" << p.sr;
+ throw std::runtime_error (err.str ());
+ }
+ int nsamps = in.getNumSamples ();
+
+ std::vector samples (nsamps);
+ if (in.getNumChannels () == 1) {
+ in.read (&samples[0], nsamps);
+ } else if (in.getNumChannels () == 2) {
+ std::vector stereo_samples (nsamps * 2);
+ in.read (&stereo_samples[0], nsamps);
+ for (unsigned i = 0; i < nsamps; ++i) {
+ samples[i] = .5 * stereo_samples[2 * i] + .5 * stereo_samples[2 * i +1];
+ }
+ } else {
+ std::stringstream err;
+ err << "invalid number of channels in << " << fname << "(mono or stereo only)";
+ throw std::runtime_error (err.str ());
+ }
+
+ std::vector buffer (p.fft_size * 2); // complex
+ std::vector window (p.fft_size);
+ makeWindow (&window[0], p.fft_size, .5, .5, 0); // hanning
+ double winEnergy = 0;
+ for (unsigned i = 0; i < p.fft_size; ++i) {
+ winEnergy += window[i] * window[i];
+ }
+ std::vector amps (p.fft_size / 2); // discard symmetry
+ std::vector freqs (p.fft_size / 2); // discard symmetry
+ for (unsigned i = 0; i < p.fft_size / 2; ++i) {
+ freqs[i] = p.sr / p.fft_size * i;
+ }
+ int ptr = 0;
+ while (ptr < nsamps) {
+ if (p.fft_size + ptr > nsamps) break; // discard incomplete frame
+
+ memset (&buffer[0], 0, sizeof (double) * 2 * p.fft_size);
+ for (unsigned i = 0; i < p.fft_size; ++i) {
+ buffer[2 * i] = window[i] * samples[i + ptr];
+ }
+ fft (&buffer[0], p.fft_size, -1);
+
+ // double max_peak = 0;
+ // int max_pos = 0;
+ // for (unsigned i = 0; i < p.fft_size; ++i) {
+ // double r = buffer[2 * i];
+ // double im = buffer[2 * i + 1];
+ // double a = sqrt (r * r + im * im);
+ // if (a > max_peak) {
+ // max_peak = a;
+ // max_pos = i;
+ // }
+ // }
+ // d.highest_peak_Hz.push_back ((p.sr / p.fft_size) * max_pos);
+
+ for (unsigned i = 0; i < p.fft_size / 2; ++i) {
+ double r = buffer[2 * i];
+ double im = buffer[2 * i + 1];
+ double a = sqrt (r * r + im * im);
+ amps[i] = a;
+ }
+ double c = speccentr (&s[0], &freqs[0], p.fft_size / 2);
+ d.highest_peak_Hz.push_back (c);
+
+ double irr = specirr (&s[0], p.fft_size / 2);
+ d.irr.push_back (irr);
+
+ double z = zcr (&samples[ptr], p.fft_size);
+ d.zcr.push_back (z);
+
+ double sum = 0;
+ for (int i = 0; i < p.fft_size; ++i) {
+ double a = samples[i + ptr];
+ sum += a * a;
+ }
+ d.energy.push_back (std::sqrt (sum / p.fft_size));
+
+ ptr += p.hop_size;
+ }
+}
+
+#endif // ANALYSIS_H
+
+// EOF
diff --git a/source/cella/features.h b/source/cella/features.h
new file mode 100644
index 0000000..9358fc6
--- /dev/null
+++ b/source/cella/features.h
@@ -0,0 +1,255 @@
+// features.h
+//
+
+#ifndef FEATURES_H
+#define FEATURES_H
+
+#include "FFT.h"
+#include "algorithms.h"
+#include
+
+// spectral descriptors ------------------------------------------------- //
+
+template
+inline T speccentr(
+ const T* amplitudes,
+ const T* frequencies,
+ int N) {
+ return centroid(amplitudes, frequencies, N);
+}
+
+template
+inline T specspread(
+ const T* amplitudes,
+ const T* frequencies,
+ int N,
+ T centroid) {
+ return std::sqrt(moment (amplitudes, frequencies, N, 2, centroid));
+}
+
+template
+inline T specskew(
+ const T* amplitudes,
+ const T* frequencies,
+ int N,
+ T centroid,
+ T spread) {
+
+ T delta = std::pow(spread, 3);
+ T tmp = moment (amplitudes, frequencies, N, 3, centroid);
+ if (delta != 0) {
+ tmp /= delta;
+ }
+
+ return tmp;
+}
+
+template
+inline T speckurt(
+ const T* amplitudes,
+ const T* frequencies,
+ int N,
+ T centroid,
+ T spread) {
+
+ T delta = std::pow(spread, 4);
+ T tmp = moment (amplitudes, frequencies, N, 4, centroid);
+ if (delta != 0) {
+ tmp /= delta;
+ }
+
+ return tmp;
+}
+
+template
+inline T specflux(T* amplitudes, T* oldAmplitudes, int N) {
+ T sf = 0; // spectral flux
+ T a = 0;
+ for (int i = 0; i < N; ++i) {
+ a = (amplitudes[i] - oldAmplitudes[i]);
+ oldAmplitudes[i] = amplitudes[i];
+ sf += a < 0 ? 0 : a; // rectification
+ //sf += a;
+ }
+
+ return sf;
+}
+
+template
+inline T specirr(T* amplitudes, int N) {
+ if (1 > N) return 0;
+ T si = 0; // spectral irregularity
+ T a = 0;
+ for (int i = 1; i < N; ++i) {
+ a = fabs(amplitudes[i] - amplitudes[i - 1]);
+ si += a;
+ }
+
+ return si;
+}
+
+template
+inline T specdecr(
+ const T* amplitudes,
+ int N) {
+ T decs = 0;
+ T den = 0;
+ if (N <= 1) return 0;
+
+ for (int index = 1; index < N; ++index) {
+ decs += (amplitudes[index] - amplitudes[0]) / (T) (index);
+ den += amplitudes[index];
+ }
+
+ if (den != 0.0) {
+ decs /= den;
+ }
+
+ return decs;
+}
+
+template
+inline T specslope(
+ const T* amplitudes,
+ const T* frequencies,
+ int N) {
+ T step = 0;
+ T sum = 0;
+ for (int i = 0; i < N; ++i) {
+ sum += amplitudes[i];
+ }
+
+ T sl = linreg(amplitudes, frequencies, N, step);
+
+ if (sum != 0.0) {
+ sl /= sum;
+ }
+
+ return sl;
+}
+
+template
+inline T specflat(
+ const T* amplitudes,
+ int N) {
+ T prod = amplitudes[0];
+ T sum = amplitudes[0];
+ for (int i = 1; i < N; ++i) {
+ sum += amplitudes[i];
+ //prod *= (amplitudes[i]);
+ prod += log(amplitudes[i]);
+ }
+
+
+ //T geom = pow (prod, 1. / N);
+ T geom = std::exp(prod / N);
+ T flat = geom;
+ if (sum != 0.0) {
+ flat /= sum;
+ } else flat = 0;
+
+ return flat;
+}
+
+template
+inline T speccrest(
+ const T* amplitudes,
+ int N) {
+ T max = amplitudes[0];
+ T sum = amplitudes[0];
+ for (int i = 1; i < N; ++i) {
+ sum += amplitudes[i];
+ if (max < amplitudes[i]) max = amplitudes[i];
+ }
+
+
+ T crest = max;
+ if (sum != 0.0) {
+ crest /= sum;
+ } else crest = 0;
+
+ return crest;
+}
+
+template
+inline T hfc(T* amplitudes, int N) {
+ T hfc = 0; // high frequency content
+ T a = 0;
+ T n = 0;
+ for (int i = 0; i < N; ++i) {
+ a = amplitudes[i];
+ hfc += (a * a) * i;
+ n += i;
+ }
+
+ return hfc / n;
+}
+
+template
+T inharmonicity(T* amp, T* freq, int size, T f0, T R, T& sumAmpl) {
+ if (f0 <= 0) return 0;
+
+ sumAmpl = 0.0;
+ T sumVariation = 0.0;
+ for (int i = 0; i < size / 2; ++i) {
+ T harmo = (i + 1) * f0;
+ T variation = std::fabs(freq[i] - harmo);
+ sumAmpl += amp[i];
+ sumVariation += variation * amp[i];
+ }
+
+ T estInharm = 0.;
+ if (sumAmpl != 0) {
+ estInharm = (2 / f0) * (sumVariation / sumAmpl) / R;
+ }
+
+ return estInharm;
+}
+
+
+// temporal descriptors ------------------------------------------------- //
+
+template
+inline T energy(const T* samples, int N, T winEnergy) {
+ T sum = 0.;
+ T en = 0.;
+ T a = 0.;
+ for (int index = 0; index < N; ++index) {
+ a = samples[index];
+ sum += a * a;
+ }
+
+ if (winEnergy > 0) {
+ en = std::sqrt(sum / winEnergy);
+ }
+
+ return en;
+}
+
+template
+static inline T zcr(const T* samples, int N) {
+ T nb_zcr = 0.0;
+ int sign1, sign2;
+
+ if (1 > N) return 0;
+
+ sign1 = (samples[0] < 0 ? -1 : 1);
+ if (samples[0] == 0) sign1 = 0;
+
+ for (int index = 1; index < N; ++index) {
+ sign2 = (samples[index] < 0 ? -1 : 1);
+ if (samples[index] == 0) sign2 = 0;
+
+ if (sign1 != sign2) {
+ ++nb_zcr;
+ }
+
+ sign1 = sign2;
+ }
+
+ return (T) (nb_zcr / N);
+}
+
+#endif // FEATURES_H
+
+// EOF
diff --git a/source/cella/filewrite.h b/source/cella/filewrite.h
new file mode 100644
index 0000000..c02bd00
--- /dev/null
+++ b/source/cella/filewrite.h
@@ -0,0 +1,10 @@
+// basic file operations
+#include
+#include
+
+void writefile(const std::string path, const std::string filename) {
+ std::ofstream myfile;
+ myfile.open (path + filename);
+ myfile << "Writing this to a file.\n";
+ myfile.close();
+}
diff --git a/source/cella/parser.h b/source/cella/parser.h
new file mode 100644
index 0000000..9511ba5
--- /dev/null
+++ b/source/cella/parser.h
@@ -0,0 +1,18 @@
+// parser.h
+//
+
+#ifndef PARSER_H
+#define PARSER_H
+
+struct Parameters {
+ double sr;
+ int fft_size;
+ int hop_size;
+ int width;
+ int height;
+ int zoom;
+};
+
+#endif // PARSER_H
+
+// EOF
diff --git a/source/cella/svg.h b/source/cella/svg.h
new file mode 100644
index 0000000..360da39
--- /dev/null
+++ b/source/cella/svg.h
@@ -0,0 +1,612 @@
+// svg_utils.h
+
+#ifndef SVG_UTILS_H
+#define SVG_UTILS_H
+
+#include
+#include
+#include
+#include
+
+#include
+
+ // Utility XML/String Functions.
+ template
+ std::string attribute(std::string const & attribute_name,
+ T const & value, std::string const & unit = "")
+ {
+ std::stringstream ss;
+ ss << attribute_name << "=\"" << value << unit << "\" ";
+ return ss.str();
+ }
+ std::string elemStart(std::string const & element_name)
+ {
+ return "\t<" + element_name + " ";
+ }
+ std::string elemEnd(std::string const & element_name)
+ {
+ return "" + element_name + ">\n";
+ }
+ std::string emptyElemEnd()
+ {
+ return "/>\n";
+ }
+
+ // Quick optional return type. This allows functions to return an invalid
+ // value if no good return is possible. The user checks for validity
+ // before using the returned value.
+ template
+ class optional
+ {
+ public:
+ optional(T const & type)
+ : valid(true), type(type) { }
+ optional() : valid(false), type(T()) { }
+ T * operator->()
+ {
+ // If we try to access an invalid value, an exception is thrown.
+ if (!valid)
+ throw std::exception();
+
+ return &type;
+ }
+ // Test for validity.
+ bool operator!() const { return !valid; }
+ private:
+ bool valid;
+ T type;
+ };
+
+ struct Dimensions
+ {
+ Dimensions(double width, double height) : width(width), height(height) { }
+ Dimensions(double combined = 0) : width(combined), height(combined) { }
+ double width;
+ double height;
+ };
+
+ struct pPoint
+ {
+ pPoint(double x = 0, double y = 0) : x(x), y(y) { }
+ double x;
+ double y;
+ };
+ optional getMinPoint(std::vector const & pPoints)
+ {
+ if (pPoints.empty())
+ return optional();
+
+ pPoint min = pPoints[0];
+ for (unsigned i = 0; i < pPoints.size(); ++i) {
+ if (pPoints[i].x < min.x)
+ min.x = pPoints[i].x;
+ if (pPoints[i].y < min.y)
+ min.y = pPoints[i].y;
+ }
+ return optional(min);
+ }
+ optional getMaxPoint(std::vector const & pPoints)
+ {
+ if (pPoints.empty())
+ return optional();
+
+ pPoint max = pPoints[0];
+ for (unsigned i = 0; i < pPoints.size(); ++i) {
+ if (pPoints[i].x > max.x)
+ max.x = pPoints[i].x;
+ if (pPoints[i].y > max.y)
+ max.y = pPoints[i].y;
+ }
+ return optional(max);
+ }
+
+ // Defines the dimensions, scale, origin, and origin offset of the document.
+ struct Layout
+ {
+ enum Origin { TopLeft, BottomLeft, TopRight, BottomRight };
+
+ Layout(Dimensions const & dimensions = Dimensions(400, 300), Origin origin = BottomLeft,
+ double scale = 1, pPoint const & origin_offset = pPoint(0, 0))
+ : dimensions(dimensions), scale(scale), origin(origin), origin_offset(origin_offset) { }
+ Dimensions dimensions;
+ double scale;
+ Origin origin;
+ pPoint origin_offset;
+ };
+
+ // Convert coordinates in user space to SVG native space.
+ double translateX(double x, Layout const & layout)
+ {
+ if (layout.origin == Layout::BottomRight || layout.origin == Layout::TopRight)
+ return layout.dimensions.width - ((x + layout.origin_offset.x) * layout.scale);
+ else
+ return (layout.origin_offset.x + x) * layout.scale;
+ }
+
+ double translateY(double y, Layout const & layout)
+ {
+ if (layout.origin == Layout::BottomLeft || layout.origin == Layout::BottomRight)
+ return layout.dimensions.height - ((y + layout.origin_offset.y) * layout.scale);
+ else
+ return (layout.origin_offset.y + y) * layout.scale;
+ }
+ double translateScale(double dimension, Layout const & layout)
+ {
+ return dimension * layout.scale;
+ }
+
+ class Serializeable
+ {
+ public:
+ Serializeable() { }
+ virtual ~Serializeable() { };
+ virtual std::string toString(Layout const & layout) const = 0;
+ };
+
+ class Color : public Serializeable
+ {
+ public:
+ enum Defaults { Transparent = -1, Aqua, Black, Blue, Brown, Cyan, Fuchsia,
+ Green, Lime, Magenta, Orange, Purple, Red, Silver, White, Yellow };
+
+ Color(int r, int g, int b) : transparent(false), red(r), green(g), blue(b) { }
+ Color(Defaults color)
+ : transparent(false), red(0), green(0), blue(0)
+ {
+ switch (color)
+ {
+ case Aqua: assign(0, 255, 255); break;
+ case Black: assign(0, 0, 0); break;
+ case Blue: assign(0, 0, 255); break;
+ case Brown: assign(165, 42, 42); break;
+ case Cyan: assign(0, 255, 255); break;
+ case Fuchsia: assign(255, 0, 255); break;
+ case Green: assign(0, 128, 0); break;
+ case Lime: assign(0, 255, 0); break;
+ case Magenta: assign(255, 0, 255); break;
+ case Orange: assign(255, 165, 0); break;
+ case Purple: assign(128, 0, 128); break;
+ case Red: assign(255, 0, 0); break;
+ case Silver: assign(192, 192, 192); break;
+ case White: assign(255, 255, 255); break;
+ case Yellow: assign(255, 255, 0); break;
+ default: transparent = true; break;
+ }
+ }
+ virtual ~Color() { }
+ std::string toString(Layout const &) const
+ {
+ std::stringstream ss;
+ if (transparent)
+ ss << "transparent";
+ else
+ ss << "rgb(" << red << "," << green << "," << blue << ")";
+ return ss.str();
+ }
+ private:
+ bool transparent;
+ int red;
+ int green;
+ int blue;
+
+ void assign(int r, int g, int b)
+ {
+ red = r;
+ green = g;
+ blue = b;
+ }
+ };
+
+ class Fill : public Serializeable
+ {
+ public:
+ Fill(Color::Defaults color) : color(color) { }
+ Fill(Color color = Color::Transparent)
+ : color(color) { }
+ std::string toString(Layout const & layout) const
+ {
+ std::stringstream ss;
+ ss << attribute("fill", color.toString(layout));
+ return ss.str();
+ }
+ private:
+ Color color;
+ };
+
+ class Stroke : public Serializeable
+ {
+ public:
+ Stroke(double width = -1, Color color = Color::Transparent)
+ : width(width), color(color) { }
+ std::string toString(Layout const & layout) const
+ {
+ // If stroke width is invalid.
+ if (width < 0)
+ return std::string();
+
+ std::stringstream ss;
+ ss << attribute("stroke-width", translateScale(width, layout)) << attribute("stroke", color.toString(layout));
+ return ss.str();
+ }
+ private:
+ double width;
+ Color color;
+ };
+
+ class Font : public Serializeable
+ {
+ public:
+ Font(double size = 12, std::string const & family = "Verdana") : size(size), family(family) { }
+ std::string toString(Layout const & layout) const
+ {
+ std::stringstream ss;
+ ss << attribute("font-size", translateScale(size, layout)) << attribute("font-family", family);
+ return ss.str();
+ }
+ private:
+ double size;
+ std::string family;
+ };
+
+ class Shape : public Serializeable
+ {
+ public:
+ Shape(Fill const & fill = Fill(), Stroke const & stroke = Stroke())
+ : fill(fill), stroke(stroke) { }
+ virtual ~Shape() { }
+ virtual std::string toString(Layout const & layout) const = 0;
+ virtual void offset(pPoint const & offset) = 0;
+ protected:
+ Fill fill;
+ Stroke stroke;
+ };
+ template
+ std::string vectorToString(std::vector collection, Layout const & layout)
+ {
+ std::string combination_str;
+ for (unsigned i = 0; i < collection.size(); ++i)
+ combination_str += collection[i].toString(layout);
+
+ return combination_str;
+ }
+
+ class Circle : public Shape
+ {
+ public:
+ Circle(pPoint const & center, double diameter, Fill const & fill,
+ Stroke const & stroke = Stroke())
+ : Shape(fill, stroke), center(center), radius(diameter / 2) { }
+ std::string toString(Layout const & layout) const
+ {
+ std::stringstream ss;
+ ss << elemStart("circle") << attribute("cx", translateX(center.x, layout))
+ << attribute("cy", translateY(center.y, layout))
+ << attribute("r", translateScale(radius, layout)) << fill.toString(layout)
+ << stroke.toString(layout) << emptyElemEnd();
+ return ss.str();
+ }
+ void offset(pPoint const & offset)
+ {
+ center.x += offset.x;
+ center.y += offset.y;
+ }
+ private:
+ pPoint center;
+ double radius;
+ };
+
+ class Elipse : public Shape
+ {
+ public:
+ Elipse(pPoint const & center, double width, double height,
+ Fill const & fill = Fill(), Stroke const & stroke = Stroke())
+ : Shape(fill, stroke), center(center), radius_width(width / 2),
+ radius_height(height / 2) { }
+ std::string toString(Layout const & layout) const
+ {
+ std::stringstream ss;
+ ss << elemStart("ellipse") << attribute("cx", translateX(center.x, layout))
+ << attribute("cy", translateY(center.y, layout))
+ << attribute("rx", translateScale(radius_width, layout))
+ << attribute("ry", translateScale(radius_height, layout))
+ << fill.toString(layout) << stroke.toString(layout) << emptyElemEnd();
+ return ss.str();
+ }
+ void offset(pPoint const & offset)
+ {
+ center.x += offset.x;
+ center.y += offset.y;
+ }
+ private:
+ pPoint center;
+ double radius_width;
+ double radius_height;
+ };
+
+ class Rectangle : public Shape
+ {
+ public:
+ Rectangle(pPoint const & edge, double width, double height,
+ Fill const & fill = Fill(), Stroke const & stroke = Stroke())
+ : Shape(fill, stroke), edge(edge), width(width),
+ height(height) { }
+ std::string toString(Layout const & layout) const
+ {
+ std::stringstream ss;
+ ss << elemStart("rect") << attribute("x", translateX(edge.x, layout))
+ << attribute("y", translateY(edge.y, layout))
+ << attribute("width", translateScale(width, layout))
+ << attribute("height", translateScale(height, layout))
+ << fill.toString(layout) << stroke.toString(layout) << emptyElemEnd();
+ return ss.str();
+ }
+ void offset(pPoint const & offset)
+ {
+ edge.x += offset.x;
+ edge.y += offset.y;
+ }
+ private:
+ pPoint edge;
+ double width;
+ double height;
+ };
+
+ class Line : public Shape
+ {
+ public:
+ Line(pPoint const & start_pPoint, pPoint const & end_pPoint,
+ Stroke const & stroke = Stroke())
+ : Shape(Fill(), stroke), start_pPoint(start_pPoint),
+ end_pPoint(end_pPoint) { }
+ std::string toString(Layout const & layout) const
+ {
+ std::stringstream ss;
+ ss << elemStart("line") << attribute("x1", translateX(start_pPoint.x, layout))
+ << attribute("y1", translateY(start_pPoint.y, layout))
+ << attribute("x2", translateX(end_pPoint.x, layout))
+ << attribute("y2", translateY(end_pPoint.y, layout))
+ << stroke.toString(layout) << emptyElemEnd();
+ return ss.str();
+ }
+ void offset(pPoint const & offset)
+ {
+ start_pPoint.x += offset.x;
+ start_pPoint.y += offset.y;
+
+ end_pPoint.x += offset.x;
+ end_pPoint.y += offset.y;
+ }
+ private:
+ pPoint start_pPoint;
+ pPoint end_pPoint;
+ };
+
+ class Polygon : public Shape
+ {
+ public:
+ Polygon(Fill const & fill = Fill(), Stroke const & stroke = Stroke())
+ : Shape(fill, stroke) { }
+ Polygon(Stroke const & stroke = Stroke()) : Shape(Color::Transparent, stroke) { }
+ Polygon & operator<<(pPoint const & pPoint)
+ {
+ pPoints.push_back(pPoint);
+ return *this;
+ }
+ std::string toString(Layout const & layout) const
+ {
+ std::stringstream ss;
+ ss << elemStart("polygon");
+
+ ss << "points=\"";
+ for (unsigned i = 0; i < pPoints.size(); ++i)
+ ss << translateX(pPoints[i].x, layout) << "," << translateY(pPoints[i].y, layout) << " ";
+ ss << "\" ";
+
+ ss << fill.toString(layout) << stroke.toString(layout) << emptyElemEnd();
+ return ss.str();
+ }
+ void offset(pPoint const & offset)
+ {
+ for (unsigned i = 0; i < pPoints.size(); ++i) {
+ pPoints[i].x += offset.x;
+ pPoints[i].y += offset.y;
+ }
+ }
+ private:
+ std::vector pPoints;
+ };
+
+ class Polyline : public Shape
+ {
+ public:
+ Polyline(Fill const & fill = Fill(), Stroke const & stroke = Stroke())
+ : Shape(fill, stroke) { }
+ Polyline(Stroke const & stroke = Stroke()) : Shape(Color::Transparent, stroke) { }
+ Polyline(std::vector const & pPoints,
+ Fill const & fill = Fill(), Stroke const & stroke = Stroke())
+ : Shape(fill, stroke), pPoints(pPoints) { }
+ Polyline & operator<<(pPoint const & pPoint)
+ {
+ pPoints.push_back(pPoint);
+ return *this;
+ }
+ std::string toString(Layout const & layout) const
+ {
+ std::stringstream ss;
+ ss << elemStart("polyline");
+
+ ss << "Points=\"";
+ for (unsigned i = 0; i < pPoints.size(); ++i)
+ ss << translateX(pPoints[i].x, layout) << "," << translateY(pPoints[i].y, layout) << " ";
+ ss << "\" ";
+
+ ss << fill.toString(layout) << stroke.toString(layout) << emptyElemEnd();
+ return ss.str();
+ }
+ void offset(pPoint const & offset)
+ {
+ for (unsigned i = 0; i < pPoints.size(); ++i) {
+ pPoints[i].x += offset.x;
+ pPoints[i].y += offset.y;
+ }
+ }
+ std::vector pPoints;
+ };
+
+ class Text : public Shape
+ {
+ public:
+ Text(pPoint const & origin, std::string const & content, Fill const & fill = Fill(),
+ Font const & font = Font(), Stroke const & stroke = Stroke())
+ : Shape(fill, stroke), origin(origin), content(content), font(font) { }
+ std::string toString(Layout const & layout) const
+ {
+ std::stringstream ss;
+ ss << elemStart("text") << attribute("x", translateX(origin.x, layout))
+ << attribute("y", translateY(origin.y, layout))
+ << fill.toString(layout) << stroke.toString(layout) << font.toString(layout)
+ << ">" << content << elemEnd("text");
+ return ss.str();
+ }
+ void offset(pPoint const & offset)
+ {
+ origin.x += offset.x;
+ origin.y += offset.y;
+ }
+ private:
+ pPoint origin;
+ std::string content;
+ Font font;
+ };
+
+ // Sample charting class.
+ class LineChart : public Shape
+ {
+ public:
+ LineChart(Dimensions margin = Dimensions(), double scale = 1,
+ Stroke const & axis_stroke = Stroke(.5, Color::Purple))
+ : axis_stroke(axis_stroke), margin(margin), scale(scale) { }
+ LineChart & operator<<(Polyline const & polyline)
+ {
+ if (polyline.pPoints.empty())
+ return *this;
+
+ polylines.push_back(polyline);
+ return *this;
+ }
+ std::string toString(Layout const & layout) const
+ {
+ if (polylines.empty())
+ return "";
+
+ std::string ret;
+ for (unsigned i = 0; i < polylines.size(); ++i)
+ ret += polylineToString(polylines[i], layout);
+
+ return ret + axisString(layout);
+ }
+ void offset(pPoint const & offset)
+ {
+ for (unsigned i = 0; i < polylines.size(); ++i)
+ polylines[i].offset(offset);
+ }
+ private:
+ Stroke axis_stroke;
+ Dimensions margin;
+ double scale;
+ std::vector polylines;
+
+ optional getDimensions() const
+ {
+ if (polylines.empty())
+ return optional();
+
+ optional min = getMinPoint(polylines[0].pPoints);
+ optional max = getMaxPoint(polylines[0].pPoints);
+ for (unsigned i = 0; i < polylines.size(); ++i) {
+ if (getMinPoint(polylines[i].pPoints)->x < min->x)
+ min->x = getMinPoint(polylines[i].pPoints)->x;
+ if (getMinPoint(polylines[i].pPoints)->y < min->y)
+ min->y = getMinPoint(polylines[i].pPoints)->y;
+ if (getMaxPoint(polylines[i].pPoints)->x > max->x)
+ max->x = getMaxPoint(polylines[i].pPoints)->x;
+ if (getMaxPoint(polylines[i].pPoints)->y > max->y)
+ max->y = getMaxPoint(polylines[i].pPoints)->y;
+ }
+
+ return optional(Dimensions(max->x - min->x, max->y - min->y));
+ }
+ std::string axisString(Layout const & layout) const
+ {
+ optional dimensions = getDimensions();
+ if (!dimensions)
+ return "";
+
+ // Make the axis 10% wider and higher than the data Points.
+ double width = dimensions->width * 1.1;
+ double height = dimensions->height * 1.1;
+
+ // Draw the axis.
+ Polyline axis(Color::Transparent, axis_stroke);
+ axis << pPoint(margin.width, margin.height + height) << pPoint(margin.width, margin.height)
+ << pPoint(margin.width + width, margin.height);
+
+ return axis.toString(layout);
+ }
+ std::string polylineToString(Polyline const & polyline, Layout const & layout) const
+ {
+ Polyline shifted_polyline = polyline;
+ shifted_polyline.offset(pPoint(margin.width, margin.height));
+
+ std::vector vertices;
+ for (unsigned i = 0; i < shifted_polyline.pPoints.size(); ++i)
+ vertices.push_back(Circle(shifted_polyline.pPoints[i], getDimensions()->height / 30.0, Color::Black));
+
+ return shifted_polyline.toString(layout) + vectorToString(vertices, layout);
+ }
+ };
+
+ class Document
+ {
+ public:
+ Document(std::string const & file_name, Layout layout = Layout())
+ : file_name(file_name), layout(layout) { }
+
+ Document & operator<<(Shape const & shape)
+ {
+ body_nodes_str += shape.toString(layout);
+ return *this;
+ }
+ std::string toString() const
+ {
+ std::stringstream ss;
+ ss << "\n\n