Some initial changes that are long overdue part 1

1 files changed, 290 insertions, 0 deletions

diff --git a/libs/maxiFFT.cpp b/libs/maxiFFT.cpp
new file mode 100644
index 0000000..c648992
--- /dev/null
+++ b/libs/maxiFFT.cpp
@@ -0,0 +1,290 @@
+/*
+ *  maximilian
+ *  platform independent synthesis library using portaudio or rtaudio
+ *
+ *  Created by Mick Grierson on 29/12/2009.
+ *  Copyright 2009 Mick Grierson & Strangeloop Limited. All rights reserved.
+ *	Thanks to the Goldsmiths Creative Computing Team.
+ *	Special thanks to Arturo Castro for the PortAudio implementation.
+ * 
+ *	Permission is hereby granted, free of charge, to any person
+ *	obtaining a copy of this software and associated documentation
+ *	files (the "Software"), to deal in the Software without
+ *	restriction, including without limitation the rights to use,
+ *	copy, modify, merge, publish, distribute, sublicense, and/or sell
+ *	copies of the Software, and to permit persons to whom the
+ *	Software is furnished to do so, subject to the following
+ *	conditions:
+ *	
+ *	The above copyright notice and this permission notice shall be
+ *	included in all copies or substantial portions of the Software.
+ *
+ *	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,	
+ *	EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+ *	OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+ *	NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+ *	HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+ *	WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ *	FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+ *	OTHER DEALINGS IN THE SOFTWARE.
+ *
+ */
+
+#include "maxiFFT.h"
+#include "maximilian.h"
+#include <iostream>
+#include "math.h"
+
+using namespace std;
+
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+//F F T
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+void maxiFFT::setup(int _fftSize, int _windowSize, int _hopSize) {
+	_fft = new fft(_fftSize);
+	fftSize = _fftSize;
+	windowSize = _windowSize;
+	bins = fftSize / 2;
+	hopSize = _hopSize;
+	buffer = (float *) malloc(fftSize * sizeof(float));
+	magnitudes = (float *) malloc(bins * sizeof(float));
+	magnitudesDB = (float *) malloc(bins * sizeof(float));
+	phases = (float *) malloc(bins * sizeof(float));
+	avgPower = new float;
+	memset(buffer, 0, fftSize * sizeof(float));
+	memset(magnitudes, 0, bins * sizeof(float));
+	memset(magnitudesDB, 0, bins * sizeof(float));
+	memset(phases, 0, bins * sizeof(float));
+	*avgPower = 0;
+	pos =windowSize - hopSize;
+	newFFT = 0;
+	window = (float*) malloc(fftSize * sizeof(float));
+	memset(window, 0, fftSize * sizeof(float));
+	fft::genWindow(3, windowSize, window);
+}
+
+bool maxiFFT::process(float value) {
+	//add value to buffer at current pos
+	buffer[pos++] = value;
+	//if buffer full, run fft
+	newFFT = pos == windowSize;
+	if (newFFT) {
+#if defined(__APPLE_CC__) && !defined(_NO_VDSP)
+		_fft->powerSpectrum_vdsp(0, buffer, window, magnitudes, phases);
+#else
+		_fft->powerSpectrum(0, buffer, window, magnitudes, phases);		
+#endif
+		//shift buffer back by one hop size
+		memcpy(buffer, buffer + hopSize, (windowSize - hopSize) * sizeof(float));
+		//set the new data to zero (is this necessary?, it will get overwritten)
+		memset(buffer + windowSize - hopSize, 0, hopSize * sizeof(float));
+		//reset pos to start of hop
+		pos= windowSize - hopSize;
+	}
+	return newFFT;
+}
+
+float* maxiFFT::magsToDB() {
+#if defined(__APPLE_CC__) && !defined(_NO_VDSP)
+	_fft->convToDB_vdsp(magnitudes, magnitudesDB);
+#else
+	_fft->convToDB(magnitudes, magnitudesDB);
+#endif	
+	return magnitudesDB;
+}
+
+float maxiFFT::spectralFlatness() {
+	float geometricMean=0, arithmaticMean=0;
+	for(int i=0; i < bins; i++) {
+		if (magnitudes[i] != 0)
+			geometricMean += logf(magnitudes[i]);
+		arithmaticMean += magnitudes[i];
+	}
+	geometricMean = expf(geometricMean / (float)bins);
+	arithmaticMean /= (float)bins;
+	return arithmaticMean !=0 ?  geometricMean / arithmaticMean : 0;
+}
+
+float maxiFFT::spectralCentroid() {
+	float x=0, y=0;
+	for(int i=0; i < bins; i++) {
+		x += fabs(magnitudes[i]) * i;
+		y += fabs(magnitudes[i]);
+	}
+	return y != 0 ? x / y * ((float) maxiSettings::sampleRate / fftSize) : 0;
+}
+
+
+
+maxiFFT::~maxiFFT() {
+	delete _fft;
+	if (buffer)
+		delete[] buffer,magnitudes,phases,window, avgPower, magnitudesDB;
+}
+
+
+
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+//I N V E R S E  F F T
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+void maxiIFFT::setup(int _fftSize, int _windowSize, int _hopSize) {
+	_fft = new fft(_fftSize);	
+	fftSize = _fftSize;
+	windowSize = _windowSize;
+	bins = fftSize / 2;
+	hopSize = _hopSize;
+	buffer = (float *) malloc(fftSize * sizeof(float));
+	ifftOut = (float *) malloc(fftSize * sizeof(float));
+	memset(buffer, 0, windowSize * sizeof(float));
+	pos =0;
+	window = (float*) malloc(fftSize * sizeof(float));
+	memset(window, 0, fftSize * sizeof(float));
+	fft::genWindow(3, windowSize, window);
+}
+
+float maxiIFFT::process(float *magnitudes, float *phases) {
+	if (0==pos) {
+		//do ifft
+		memset(ifftOut, 0, fftSize * sizeof(float));
+#if defined(__APPLE_CC__) && !defined(_NO_VDSP)
+		_fft->inversePowerSpectrum_vdsp(0, ifftOut, window, magnitudes, phases);
+#else
+		_fft->inversePowerSpectrum(0, ifftOut, window, magnitudes, phases);
+#endif
+		//add to output
+		//shift back by one hop
+		memcpy(buffer, buffer+hopSize, (fftSize - hopSize) * sizeof(float));
+		//clear the end chunk
+		memset(buffer + (fftSize - hopSize), 0, hopSize * sizeof(float)); 
+		//merge new output
+		for(int i=0; i < fftSize; i++) {
+			buffer[i] += ifftOut[i];
+		}
+	}
+	
+	nextValue = buffer[pos];
+	//limit the values, this alg seems to spike occasionally (and break the audio drivers)
+	if (nextValue > 0.99999f) nextValue = 0.99999f;
+	if (nextValue < -0.99999f) nextValue = -0.99999f;
+	if (hopSize == ++pos ) {
+		pos=0;
+	}
+	
+	return nextValue;
+}
+
+maxiIFFT::~maxiIFFT() {
+	delete _fft;
+	delete[] ifftOut, buffer, window;
+};
+
+
+
+
+
+
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+//O C T A V E  A N A L Y S E R
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+
+void maxiFFTOctaveAnalyzer::setup(float samplingRate, int nBandsInTheFFT, int nAveragesPerOctave){
+	
+    samplingRate = samplingRate;
+    nSpectrum = nBandsInTheFFT;
+    spectrumFrequencySpan = (samplingRate / 2.0f) / (float)(nSpectrum);
+    nAverages = nBandsInTheFFT;
+    // fe:  2f for octave bands, sqrt(2) for half-octave bands, cuberoot(2) for third-octave bands, etc
+    if (nAveragesPerOctave==0) // um, wtf?
+		nAveragesPerOctave = 1;
+//    nAveragesPerOctave = nAveragesPerOctave;
+    averageFrequencyIncrement = pow(2.0f, 1.0f/(float)(nAveragesPerOctave));
+    // this isn't currently configurable (used once here then no effect), but here's some reasoning:
+    // 43 is a good value if you want to approximate "computer" octaves: 44100/2/2/2/2/2/2/2/2/2/2
+    // 55 is a good value if you'd rather approximate A-440 octaves: 440/2/2/2
+    // 65 is a good value if you'd rather approximate "middle-C" octaves:  ~262/2/2
+    // you could easily double it if you felt the lowest band was just rumble noise (as it probably is)
+    // but don't go much smaller unless you have a huge fft window size (see below for more why)
+    // keep in mind, if you change it, that the number of actual bands may change +/-1, and
+    // for some values, the last averaging band may not be very useful (may extend above nyquist)
+    firstOctaveFrequency = 55.0f;
+    // for each spectrum[] bin, calculate the mapping into the appropriate average[] bin.
+    // this gives us roughly log-sized averaging bins, subject to how "fine" the spectrum bins are.
+    // with more spectrum bins, you can better map into the averaging bins (especially at low
+    // frequencies) or use more averaging bins per octave.  with an fft window size of 2048,
+    // sampling rate of 44100, and first octave around 55, that's about enough to do half-octave
+    // analysis.  if you don't have enough spectrum bins to map adequately into averaging bins
+    // at the requested number per octave then you'll end up with "empty" averaging bins, where
+    // there is no spectrum available to map into it.  (so... if you have "nonreactive" averages,
+    // either increase fft buffer size, or decrease number of averages per octave, etc)
+    spe2avg = new int[nSpectrum];
+    int avgidx = 0;
+    float averageFreq = firstOctaveFrequency; // the "top" of the first averaging bin
+    // we're looking for the "top" of the first spectrum bin, and i'm just sort of
+    // guessing that this is where it is (or possibly spectrumFrequencySpan/2?)
+    // ... either way it's probably close enough for these purposes
+    float spectrumFreq = spectrumFrequencySpan;
+    for (int speidx=0; speidx < nSpectrum; speidx++) {
+		while (spectrumFreq > averageFreq) {
+			avgidx++;
+			averageFreq *= averageFrequencyIncrement;
+		}
+		spe2avg[speidx] = avgidx;
+		spectrumFreq += spectrumFrequencySpan;
+    }
+    nAverages = avgidx;
+    averages = new float[nAverages];
+    peaks = new float[nAverages];
+    peakHoldTimes = new int[nAverages];
+    peakHoldTime = 0; // arbitrary
+    peakDecayRate = 0.9f; // arbitrary
+    linearEQIntercept = 1.0f; // unity -- no eq by default
+    linearEQSlope = 0.0f; // unity -- no eq by default
+}
+
+void maxiFFTOctaveAnalyzer::calculate(float * fftData){
+	
+	int last_avgidx = 0; // tracks when we've crossed into a new averaging bin, so store current average
+    float sum = 0.0f; // running total of spectrum data
+    int count = 0; // count of spectrums accumulated (for averaging)
+    for (int speidx=0; speidx < nSpectrum; speidx++) {
+		count++;
+		sum += fftData[speidx] * (linearEQIntercept + (float)(speidx) * linearEQSlope);
+		int avgidx = spe2avg[speidx];
+		if (avgidx != last_avgidx) {
+			
+			for (int j = last_avgidx; j < avgidx; j++){
+				averages[j] = sum / (float)(count);
+			}
+			count = 0;
+			sum = 0.0f;
+		}
+		last_avgidx = avgidx;
+    }
+    // the last average was probably not calculated...
+    if ((count > 0) && (last_avgidx < nAverages)){
+		averages[last_avgidx] = sum / (float)(count);
+	}
+	
+    // update the peaks separately
+    for (int i=0; i < nAverages; i++) {
+		if (averages[i] >= peaks[i]) {
+			// save new peak level, also reset the hold timer
+			peaks[i] = averages[i];
+			peakHoldTimes[i] = peakHoldTime;
+		} else {
+			// current average does not exceed peak, so hold or decay the peak
+			if (peakHoldTimes[i] > 0) {
+				peakHoldTimes[i]--;
+			} else {
+				peaks[i] *= peakDecayRate;
+			}
+		}
+    }
+}