diff --git a/samplebrain/Makefile.in b/samplebrain/Makefile.in index 610dfa8..f84305b 100644 --- a/samplebrain/Makefile.in +++ b/samplebrain/Makefile.in @@ -5,12 +5,16 @@ SRCS := src/fft.cpp \ src/brain.cpp \ src/brain_block.cpp \ src/libmfcc.cpp \ - src/main.cpp + src/main.cpp \ + src/mfcc.cpp \ + src/aquila/filter/MelFilterBank.cpp \ + src/aquila/filter/MelFilter.cpp \ + src/aquila/transform/Dct.cpp TARGET_SRCS := src/main.cpp # for the minute, go out and up to link to the vision lib -CCFLAGS = @CFLAGS@ -ffast-math -Wno-unused -Isrc +CCFLAGS = @CFLAGS@ -std=c++11 -ffast-math -Wno-unused -Isrc LDFLAGS = @LDFLAGS@ LIBS = @LIBS@ diff --git a/samplebrain/dist.sh b/samplebrain/dist.sh index 4d12a7a..44b24b4 100755 --- a/samplebrain/dist.sh +++ b/samplebrain/dist.sh @@ -4,3 +4,4 @@ make distclean autoheader # build configure autoconf configure.ac > configure +./configure CXX=g++-4.7 diff --git a/samplebrain/src/aquila/.#global.h b/samplebrain/src/aquila/.#global.h new file mode 120000 index 0000000..9f6e7bd --- /dev/null +++ b/samplebrain/src/aquila/.#global.h @@ -0,0 +1 @@ +dave@fulmar.4739:1436183149 \ No newline at end of file diff --git a/samplebrain/src/aquila/Exceptions.h b/samplebrain/src/aquila/Exceptions.h new file mode 100644 index 0000000..974b4ee --- /dev/null +++ b/samplebrain/src/aquila/Exceptions.h @@ -0,0 +1,83 @@ +/** + * @file Exceptions.h + * + * Exception class definitions. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 2.0.0 + */ + +#ifndef EXCEPTIONS_H +#define EXCEPTIONS_H + +#include "global.h" +#include +#include + +namespace Aquila +{ + /** + * Base exception class of the library. + * + * Class clients should rather catch exceptions of specific types, such as + * Aquila::FormatException, however it is allowed to catch Aquila::Exception + * as the last resort (but catch(...)). + */ + class AQUILA_EXPORT Exception : public std::runtime_error + { + public: + /** + * Creates an exception object. + * + * @param message exception message + */ + Exception(const std::string& message): + runtime_error(message) + { + } + }; + + /** + * Data format-related exception. + */ + class AQUILA_EXPORT FormatException : public Exception + { + public: + /** + * Creates a data format exception object. + * + * @param message exception message + */ + FormatException(const std::string& message): + Exception(message) + { + } + }; + + /** + * Runtime configuration exception. + */ + class AQUILA_EXPORT ConfigurationException : public Exception + { + public: + /** + * Creates a configuration exception object. + * + * @param message exception message + */ + ConfigurationException(const std::string& message): + Exception(message) + { + } + }; +} + +#endif // EXCEPTIONS_H diff --git a/samplebrain/src/aquila/aquila.h b/samplebrain/src/aquila/aquila.h new file mode 100644 index 0000000..00e6601 --- /dev/null +++ b/samplebrain/src/aquila/aquila.h @@ -0,0 +1,48 @@ +/** + * @file aquila.h + * + * Library "master" header - includes all component headers. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + * + * @mainpage + * + * @section what-is-aquila What is Aquila? + * Aquila is an open source and cross-platform DSP (Digital Signal Processing) + * library for C++11. + * + * Aquila provides a set of classes for common DSP operations, such as FFT, DCT, + * Mel-frequency filtering, calculating spectrograms etc. It supports reading + * and writing signals in various formats, such as raw binary files, text files + * or WAVE audio recordings. + * + * @section motivation Motivation + * The initial goal of this project was to develop computer software capable + * of recognizing birds' songs. Since then the library was redesigned and + * extended with more general DSP tools. There are still a few major + * shortcomings, for example the lack of general purpose filter classes, but + * hopefully this will change soon. + */ + +#ifndef AQUILA_H +#define AQUILA_H + +#include "global.h" +#include "Exceptions.h" +#include "functions.h" +#include "source.h" +#include "transform.h" +#include "filter.h" +#include "ml.h" +#include "tools.h" + +#endif // AQUILA_H diff --git a/samplebrain/src/aquila/filter.h b/samplebrain/src/aquila/filter.h new file mode 100644 index 0000000..24ec8a9 --- /dev/null +++ b/samplebrain/src/aquila/filter.h @@ -0,0 +1,25 @@ +/** + * @file filter.h + * + * Convenience header that includes all filter headers. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + + +#ifndef AQUILA_FILTER_H +#define AQUILA_FILTER_H + +#include "filter/MelFilter.h" +#include "filter/MelFilterBank.h" + +#endif // AQUILA_FILTER_H diff --git a/samplebrain/src/aquila/filter/MelFilter.cpp b/samplebrain/src/aquila/filter/MelFilter.cpp new file mode 100644 index 0000000..5d4dd9b --- /dev/null +++ b/samplebrain/src/aquila/filter/MelFilter.cpp @@ -0,0 +1,157 @@ +/** + * @file MelFilter.cpp + * + * Triangular filters in Mel frequency scale. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 0.3.3 + */ + +#include "MelFilter.h" +#include +#include + +namespace Aquila +{ + /** + * Creates the filter and sets sample frequency. + * + * @param sampleFrequency sample frequency in Hz + */ + MelFilter::MelFilter(FrequencyType sampleFrequency): + m_sampleFrequency(sampleFrequency), m_spectrum() + { + } + + /** + * Move constructor. + * + * @param other other filter to be moved from + */ + MelFilter::MelFilter(MelFilter&& other): + m_sampleFrequency(other.m_sampleFrequency), + m_spectrum(std::move(other.m_spectrum)) + { + } + + /** + * Copy assignment operator. + * + * @param other filter to be copied from + * @return reference to assigned value + */ + MelFilter& MelFilter::operator=(const MelFilter& other) + { + m_sampleFrequency = other.m_sampleFrequency; + m_spectrum = other.m_spectrum; + return *this; + } + + /** + * Designs the Mel filter and creates triangular spectrum. + * + * @param filterNum which filter in a sequence it is + * @param melFilterWidth filter width in Mel scale (eg. 200) + * @param N filter spectrum size (must be the same as filtered spectrum) + */ + void MelFilter::createFilter(std::size_t filterNum, + FrequencyType melFilterWidth, std::size_t N) + { + // calculate frequencies in Mel scale + FrequencyType melMinFreq = filterNum * melFilterWidth / 2.0; + FrequencyType melCenterFreq = melMinFreq + melFilterWidth / 2.0; + FrequencyType melMaxFreq = melMinFreq + melFilterWidth; + + // convert frequencies to linear scale + FrequencyType minFreq = melToLinear(melMinFreq); + FrequencyType centerFreq = melToLinear(melCenterFreq); + FrequencyType maxFreq = melToLinear(melMaxFreq); + + // generate filter spectrum in linear scale + generateFilterSpectrum(minFreq, centerFreq, maxFreq, N); + } + + /** + * Returns a single value computed by multiplying signal spectrum with + * Mel filter spectrum and summing all the products. + * + * @param dataSpectrum complex signal spectrum + * @return dot product of the spectra + */ + double MelFilter::apply(const SpectrumType& dataSpectrum) const + { + double value = 0.0; + const std::size_t N = dataSpectrum.size(); + for (std::size_t i = 0; i < N; ++i) + { + value += std::abs(dataSpectrum[i]) * m_spectrum[i]; + } + return value; + } + + /** + * Generates a vector of values shaped as a triangular filter. + * + * ^ [2] + * | /\ + * | / \ + * | / \ + * |_____________________/ \__________________ + * +--------------------[1]----[3]---------------------> f + * + * @param minFreq frequency at [1] in linear scale + * @param centerFreq frequency at [2] in linear scale + * @param maxFreq frequency at [3] in linear scale + * @param N length of the spectrum + */ + void MelFilter::generateFilterSpectrum(FrequencyType minFreq, + FrequencyType centerFreq, + FrequencyType maxFreq, std::size_t N) + { + m_spectrum.clear(); + m_spectrum.resize(N, 0.0); + + // find spectral peak positions corresponding to frequencies + std::size_t minPos = static_cast(N * minFreq / m_sampleFrequency); + std::size_t maxPos = static_cast(N * maxFreq / m_sampleFrequency); + // limit maxPos not to write out of bounds of vector storage + maxPos = std::min(maxPos, N - 1); + if (maxPos <= minPos) { + return; + } + + const double max = 1.0; + + // outside the triangle spectrum values are 0, guaranteed by + // earlier call to resize + for (std::size_t k = minPos; k <= maxPos; ++k) + { + Aquila::FrequencyType currentFreq = (k * m_sampleFrequency) / N; + if (currentFreq < minFreq) + { + continue; + } + if (currentFreq < centerFreq) + { + // in the triangle on the ascending slope + m_spectrum[k] = (currentFreq - minFreq) * max / (centerFreq - minFreq); + } + else + { + if (currentFreq < maxFreq) + { + // in the triangle on the descending slope + m_spectrum[k] = (maxFreq - currentFreq) * max / (maxFreq - centerFreq); + } + } + } + } +} diff --git a/samplebrain/src/aquila/filter/MelFilter.dep b/samplebrain/src/aquila/filter/MelFilter.dep new file mode 100644 index 0000000..a8f6135 --- /dev/null +++ b/samplebrain/src/aquila/filter/MelFilter.dep @@ -0,0 +1,2 @@ +MelFilter.o: src/aquila/filter/MelFilter.cpp \ + src/aquila/filter/MelFilter.h src/aquila/filter/../global.h diff --git a/samplebrain/src/aquila/filter/MelFilter.h b/samplebrain/src/aquila/filter/MelFilter.h new file mode 100644 index 0000000..f321846 --- /dev/null +++ b/samplebrain/src/aquila/filter/MelFilter.h @@ -0,0 +1,89 @@ +/** + * @file MelFilter.h + * + * Triangular filters in Mel frequency scale. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 0.3.3 + */ + +#ifndef MELFILTER_H +#define MELFILTER_H + +#include "../global.h" +#include +#include +#include + +namespace Aquila +{ + /** + * Encapsulation of a single Mel-frequency filter. + */ + class AQUILA_EXPORT MelFilter + { + public: + explicit MelFilter(FrequencyType sampleFrequency); + MelFilter(MelFilter&& other); + MelFilter& operator=(const MelFilter& other); + + void createFilter(std::size_t filterNum, FrequencyType melFilterWidth, + std::size_t N); + + double apply(const SpectrumType& dataSpectrum) const; + + /** + * Converts frequency from linear to Mel scale. + * + * @param linearFrequency frequency in linear scale + * @return frequency in Mel scale + */ + static FrequencyType linearToMel(FrequencyType linearFrequency) + { + return 1127.01048 * std::log(1.0 + linearFrequency / 700.0); + } + + /** + * Converts frequency from Mel to linear scale. + * + * @param melFrequency frequency in Mel scale + * @return frequency in linear scale + */ + static FrequencyType melToLinear(FrequencyType melFrequency) + { + return 700.0 * (std::exp(melFrequency / 1127.01048) - 1.0); + } + + /** + * Returns sample frequency for which the filter was designed. + * + * @return sample frequency + */ + FrequencyType getSampleFrequency() const + { + return m_sampleFrequency; + } + + private: + FrequencyType m_sampleFrequency; + + /** + * Filter spectrum (real-valued). + */ + std::vector m_spectrum; + + void generateFilterSpectrum(FrequencyType minFreq, + FrequencyType centerFreq, + FrequencyType maxFreq, std::size_t N); + }; +} + +#endif // MELFILTER_H diff --git a/samplebrain/src/aquila/filter/MelFilter.o b/samplebrain/src/aquila/filter/MelFilter.o new file mode 100644 index 0000000..3527332 Binary files /dev/null and b/samplebrain/src/aquila/filter/MelFilter.o differ diff --git a/samplebrain/src/aquila/filter/MelFilterBank.cpp b/samplebrain/src/aquila/filter/MelFilterBank.cpp new file mode 100644 index 0000000..1d52f5b --- /dev/null +++ b/samplebrain/src/aquila/filter/MelFilterBank.cpp @@ -0,0 +1,59 @@ +/** + * @file MelFilterBank.cpp + * + * A bank of number of Mel frequency triangular filters. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 0.3.3 + */ + +#include "MelFilterBank.h" + +namespace Aquila +{ + /** + * Creates all the filters in the bank. + * + * @param sampleFrequency sample frequency in Hz + * @param length spectrum size of each filter + * @param melFilterWidth filter width in Mel frequency scale + * @param bankSize number of filters in the bank + */ + MelFilterBank::MelFilterBank(FrequencyType sampleFrequency, + std::size_t length, + FrequencyType melFilterWidth, + std::size_t bankSize): + m_filters(), m_sampleFrequency(sampleFrequency), N(length) + { + m_filters.reserve(bankSize); + for (std::size_t i = 0; i < bankSize; ++i) + { + m_filters.push_back(MelFilter(m_sampleFrequency)); + m_filters[i].createFilter(i, melFilterWidth, N); + } + } + + /** + * Processes frame spectrum through all filters. + * + * @param frameSpectrum frame spectrum + * @return vector of results (one value per each filter) + */ + std::vector MelFilterBank::applyAll(const SpectrumType& frameSpectrum) const + { + std::vector output(size(), 0.0); + for (std::size_t i = 0; i < size(); ++i) + { + output[i] = m_filters[i].apply(frameSpectrum); + } + return output; + } +} diff --git a/samplebrain/src/aquila/filter/MelFilterBank.dep b/samplebrain/src/aquila/filter/MelFilterBank.dep new file mode 100644 index 0000000..66f6b52 --- /dev/null +++ b/samplebrain/src/aquila/filter/MelFilterBank.dep @@ -0,0 +1,3 @@ +MelFilterBank.o: src/aquila/filter/MelFilterBank.cpp \ + src/aquila/filter/MelFilterBank.h src/aquila/filter/../global.h \ + src/aquila/filter/MelFilter.h diff --git a/samplebrain/src/aquila/filter/MelFilterBank.h b/samplebrain/src/aquila/filter/MelFilterBank.h new file mode 100644 index 0000000..4497920 --- /dev/null +++ b/samplebrain/src/aquila/filter/MelFilterBank.h @@ -0,0 +1,78 @@ +/** + * @file MelFilterBank.h + * + * A bank of number of Mel frequency triangular filters. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 0.3.3 + */ + +#ifndef MELFILTERBANK_H +#define MELFILTERBANK_H + +#include "../global.h" +#include "MelFilter.h" +#include +#include + +namespace Aquila +{ + /** + * A wrapper class for a vector of triangular filters. + */ + class AQUILA_EXPORT MelFilterBank + { + public: + MelFilterBank(FrequencyType sampleFrequency, std::size_t length, + FrequencyType melFilterWidth = 200.0, + std::size_t bankSize = 24); + + std::vector applyAll(const SpectrumType &frameSpectrum) const; + + /** + * Returns sample frequency of all filters. + * + * @return sample frequency in Hz + */ + FrequencyType getSampleFrequency() const { return m_sampleFrequency; } + + /** + * Returns spectrum size of any filter spectra. + * + * @return spectrum size + */ + std::size_t getSpectrumLength() const { return N; } + + /** + * Returns the number of filters in bank. + * + * @return number of filters + */ + std::size_t size() const { return m_filters.size(); } + + private: + /** + * Vector of Mel filters. + */ + std::vector m_filters; + + /** + * Sample frequency of the filtered signal. + */ + FrequencyType m_sampleFrequency; + + /** + * Filter spectrum size (equal to zero-padded length of signal frame). + */ + std::size_t N; + }; +} +#endif // MELFILTERBANK_H diff --git a/samplebrain/src/aquila/filter/MelFilterBank.o b/samplebrain/src/aquila/filter/MelFilterBank.o new file mode 100644 index 0000000..d1f49e1 Binary files /dev/null and b/samplebrain/src/aquila/filter/MelFilterBank.o differ diff --git a/samplebrain/src/aquila/functions.h b/samplebrain/src/aquila/functions.h new file mode 100644 index 0000000..d3991b9 --- /dev/null +++ b/samplebrain/src/aquila/functions.h @@ -0,0 +1,200 @@ +/** + * @file functions.h + * + * Miscellaneous utility functions. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef FUNCTIONS_H +#define FUNCTIONS_H + +#include "global.h" +#include +#include +#include +#include + +namespace Aquila +{ + /** + * Converts the value to decibels (assuming reference value equal to 1). + * + * @param value input value + * @return value in dB + */ + template + AQUILA_EXPORT inline Numeric dB(Numeric value) + { + return 20.0 * std::log10(value); + } + + /** + * Convert the magnitude of a complex number to decibels. + * + * @param value input value (complex number) + * @return magnitude in dB + */ + AQUILA_EXPORT inline double dB(ComplexType value) + { + return dB(std::abs(value)); + } + + /** + * Converts the value to decibels, relative to the reference value. + * + * @param value input value + * @param refValue reference value + * @return value in dB, relative to reference value + */ + template + AQUILA_EXPORT inline Numeric dB(Numeric value, Numeric refValue) + { + return 20.0 * std::log10(value / refValue); + } + + /** + * Clamps (limits) the value inside a range. + * + * @param min lower limit + * @param value numver to clamp + * @param max upper limit + * @return bounded value + */ + template + AQUILA_EXPORT inline Numeric clamp(Numeric min, Numeric value, Numeric max) + { + return std::max(min, std::min(value, max)); + } + + /** + * Returns a pseudorandom value from a range. + * + * @param from lower limit + * @param to upper limit + * @return random number + */ + AQUILA_EXPORT inline int random(int from, int to) + { + return std::rand() % (to - from) + from; + } + + /** + * Returns a pseudorandom double number from 0 to 1. + */ + AQUILA_EXPORT inline double randomDouble() + { + return std::rand() / static_cast(RAND_MAX); + } + + /** + * Checks if n is an exact power of 2. + */ + template + AQUILA_EXPORT inline bool isPowerOf2(Integer n) + { + return (n > 0) && ((n & (n - 1)) == 0); + } + + /** + * Returns the smallest power of 2 greater than n. + */ + template + AQUILA_EXPORT inline Integer nextPowerOf2(Integer n) + { + if (isPowerOf2(n)) + { + return 2 * n; + } + + #ifdef _MSC_VER + size_t size_in_bits = sizeof(Integer) * 8; + #else + constexpr size_t size_in_bits = sizeof(Integer) * 8; + #endif + + for (size_t shift = 1; shift < size_in_bits; shift *= 2) + { + n |= (n >> shift); + } + return (n + 1); + } + + /** + * Prototype of distance calculating functions. + */ + typedef std::function&, + const std::vector&)> DistanceFunctionType; + + /** + * Returns Euclidean distance between two vectors. + * + * @param v1 first vector + * @param v2 second vector + * @return Euclidean distance + */ + AQUILA_EXPORT inline double euclideanDistance(const std::vector& v1, + const std::vector& v2) + { + double distance = 0.0; + for (std::size_t i = 0, size = v1.size(); i < size; i++) + { + distance += (v1[i] - v2[i])*(v1[i] - v2[i]); + } + + return std::sqrt(distance); + } + + /** + * Returns Manhattan (taxicab) distance between two vectors. + * + * @param v1 first vector + * @param v2 second vector + * @return Manhattan distance + */ + AQUILA_EXPORT inline double manhattanDistance(const std::vector& v1, + const std::vector& v2) + { + double distance = 0.0; + for (std::size_t i = 0, size = v1.size(); i < size; i++) + { + distance += std::abs(v1[i] - v2[i]); + } + + return distance; + } + + /** + * Returns Chebyshev distance between two vectors. + * + * @param v1 first vector + * @param v2 second vector + * @return Chebyshev distance + */ + AQUILA_EXPORT inline double chebyshevDistance(const std::vector& v1, + const std::vector& v2) + { + double distance = 0.0, max = 0.0; + for (std::size_t i = 0, size = v1.size(); i < size; i++) + { + distance = std::abs(v1[i] - v2[i]); + if (distance > max) + { + max = distance; + } + } + + return max; + } +} + +#endif // FUNCTIONS_H diff --git a/samplebrain/src/aquila/global.h b/samplebrain/src/aquila/global.h new file mode 100644 index 0000000..666c417 --- /dev/null +++ b/samplebrain/src/aquila/global.h @@ -0,0 +1,70 @@ +/** + * @file global.h + * + * Global library typedefs and constants. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 2.4.1 + */ + +#ifndef GLOBAL_H +#define GLOBAL_H + +#include +#include + +#if defined (_WIN32) && defined(BUILD_SHARED_LIBS) +# if defined(Aquila_EXPORTS) +# define AQUILA_EXPORT __declspec(dllexport) +# else +# define AQUILA_EXPORT __declspec(dllimport) +# endif +#else +# define AQUILA_EXPORT +#endif + +/** + * Main library namespace. + */ +namespace Aquila +{ + /** + * Library version in an easily comparable format. + */ + const long VERSION_NUMBER = 0x300000; + + /** + * Library version as a string. + */ + const char* const VERSION_STRING = "3.0.0-dev"; + + /** + * Sample value type. + */ + typedef double SampleType; + + /** + * Sample frequency type. + */ + typedef double FrequencyType; + + /** + * Our standard complex number type, using double precision. + */ + typedef std::complex ComplexType; + + /** + * Spectrum type - a vector of complex values. + */ + typedef std::vector SpectrumType; +} + +#endif // GLOBAL_H diff --git a/samplebrain/src/aquila/ml.h b/samplebrain/src/aquila/ml.h new file mode 100644 index 0000000..9abf80c --- /dev/null +++ b/samplebrain/src/aquila/ml.h @@ -0,0 +1,25 @@ +/** + * @file ml.h + * + * Convenience header that includes all machine learning-related headers. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + + +#ifndef AQUILA_ML_H +#define AQUILA_ML_H + +#include "ml/DtwPoint.h" +#include "ml/Dtw.h" + +#endif // AQUILA_ML_H diff --git a/samplebrain/src/aquila/ml/Dtw.cpp b/samplebrain/src/aquila/ml/Dtw.cpp new file mode 100644 index 0000000..5ce2b67 --- /dev/null +++ b/samplebrain/src/aquila/ml/Dtw.cpp @@ -0,0 +1,117 @@ +/** + * @file Dtw.cpp + * + * An implementation of the Dynamic Time Warping algorithm. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 0.5.7 + */ + +#include "Dtw.h" + +namespace Aquila +{ + /** + * Computes the distance between two sets of data. + * + * @param from first vector of features + * @param to second vector of features + * @return double DTW distance + */ + double Dtw::getDistance(const DtwDataType& from, const DtwDataType& to) + { + m_fromSize = from.size(); + m_toSize = to.size(); + + // fill the local distances array + m_points.clear(); + m_points.resize(m_fromSize); + for (std::size_t i = 0; i < m_fromSize; ++i) + { + m_points[i].reserve(m_toSize); + for (std::size_t j = 0; j < m_toSize; ++j) + { + // use emplace_back, once all compilers support it correctly + m_points[i].push_back( + DtwPoint(i, j, m_distanceFunction(from[i], to[j]) + )); + } + } + + // the actual pathfinding algorithm + DtwPoint *top = nullptr, *center = nullptr, *bottom = nullptr, *previous = nullptr; + for (std::size_t i = 1; i < m_fromSize; ++i) + { + for (std::size_t j = 1; j < m_toSize; ++j) + { + center = &m_points[i - 1][j - 1]; + if (Neighbors == m_passType) + { + top = &m_points[i - 1][j]; + bottom = &m_points[i][j - 1]; + } + else // Diagonals + { + if (i > 1 && j > 1) + { + top = &m_points[i - 2][j - 1]; + bottom = &m_points[i - 1][j - 2]; + } + else + { + top = &m_points[i - 1][j]; + bottom = &m_points[i][j - 1]; + } + } + + if (top->dAccumulated < center->dAccumulated) + { + previous = top; + } + else + { + previous = center; + } + + if (bottom->dAccumulated < previous->dAccumulated) + { + previous = bottom; + } + + m_points[i][j].dAccumulated = m_points[i][j].dLocal + previous->dAccumulated; + m_points[i][j].previous = previous; + } + } + + return getFinalPoint().dAccumulated; + } + + /** + * Returns the lowest-cost path in the DTW array. + * + * @return path + */ + DtwPathType Dtw::getPath() const + { + DtwPathType path; + DtwPoint finalPoint = getFinalPoint(); + DtwPoint* point = &finalPoint; + + path.push_back(*point); + while(point->previous) + { + point = point->previous; + path.push_back(*point); + } + + return path; + } +} diff --git a/samplebrain/src/aquila/ml/Dtw.h b/samplebrain/src/aquila/ml/Dtw.h new file mode 100644 index 0000000..35e060b --- /dev/null +++ b/samplebrain/src/aquila/ml/Dtw.h @@ -0,0 +1,116 @@ +/** + * @file Dtw.h + * + * An implementation of the Dynamic Time Warping algorithm. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 0.5.7 + */ + +#ifndef DTW_H +#define DTW_H + +#include "../global.h" +#include "../functions.h" +#include "DtwPoint.h" +#include +#include + +namespace Aquila +{ + /** + * Type of compared data - vectors of features, which themselves are + * vectors of doubles. + */ + typedef std::vector> DtwDataType; + + /** + * Type of DTW point array. + */ + typedef std::vector> DtwPointsArrayType; + + /** + * Lowest-cost path is a vector of points. + */ + typedef std::vector DtwPathType; + + /** + * Dynamic Time Warping implementation. + */ + class AQUILA_EXPORT Dtw + { + public: + /** + * Type of lowest-cost passes between points. + */ + enum PassType {Neighbors, Diagonals}; + + /** + * Creates the DTW algorithm wrapper object. + * + * @param distanceFunction which function to use for calculating distance + * @param passType pass type - how to move through distance array + */ + Dtw(DistanceFunctionType distanceFunction = euclideanDistance, + PassType passType = Neighbors): + m_distanceFunction(distanceFunction), m_passType(passType), + m_points() + { + } + + double getDistance(const DtwDataType& from, const DtwDataType& to); + + /** + * Returns a const reference to the point array. + * + * @return DTW points + */ + const DtwPointsArrayType& getPoints() const + { + return m_points; + } + + /** + * Returns the final point on the DTW path (in the top right corner). + * + * @return a DTW point + */ + DtwPoint getFinalPoint() const + { + return m_points[m_fromSize - 1][m_toSize - 1]; + } + + DtwPathType getPath() const; + + private: + /** + * Distance definition used in DTW (eg. Euclidean, Manhattan etc). + */ + DistanceFunctionType m_distanceFunction; + + /** + * Type of passes between points. + */ + PassType m_passType; + + /** + * Array of DTW points. + */ + DtwPointsArrayType m_points; + + /** + * Coordinates of the top right corner of the points array. + */ + std::size_t m_fromSize, m_toSize; + }; +} + +#endif // DTW_H diff --git a/samplebrain/src/aquila/ml/DtwPoint.h b/samplebrain/src/aquila/ml/DtwPoint.h new file mode 100644 index 0000000..e11a195 --- /dev/null +++ b/samplebrain/src/aquila/ml/DtwPoint.h @@ -0,0 +1,81 @@ +/** + * @file DtwPoint.h + * + * A single point of the Dynamic Time Warping algorithm. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 0.5.7 + */ + +#ifndef DTWPOINT_H +#define DTWPOINT_H + +#include "../global.h" +#include + +namespace Aquila +{ + /** + * A struct representing a single point in the DTW array. + */ + struct AQUILA_EXPORT DtwPoint + { + /** + * Creates the point with default values. + */ + DtwPoint(): + x(0), y(0), dLocal(0.0), dAccumulated(0.0), previous(0) + { + } + + /** + * Creates the point and associates it with given coordinates. + * + * @param x_ x coordinate in DTW array + * @param y_ y coordinate in DTW array + * @param distanceLocal value of local distance at (x, y) + */ + DtwPoint(std::size_t x_, std::size_t y_, double distanceLocal = 0.0): + x(x_), y(y_), dLocal(distanceLocal), + // at the edges set accumulated distance to local. otherwise 0 + dAccumulated((0 == x || 0 == y) ? dLocal : 0.0), + previous(0) + { + } + + /** + * X coordinate of the point in the DTW array. + */ + std::size_t x; + + /** + * Y coordinate of the point in the DTW array. + */ + std::size_t y; + + /** + * Local distance at this point. + */ + double dLocal; + + /** + * Accumulated distance at this point. + */ + double dAccumulated; + + /** + * Non-owning pointer to previous point in the lowest-cost path. + */ + DtwPoint* previous; + }; +} + +#endif // DTWPOINT_H diff --git a/samplebrain/src/aquila/source.h b/samplebrain/src/aquila/source.h new file mode 100644 index 0000000..fe39dc8 --- /dev/null +++ b/samplebrain/src/aquila/source.h @@ -0,0 +1,43 @@ +/** + * @file source.h + * + * Convenience header that includes all signal source-related headers. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + + +#ifndef AQUILA_SOURCE_H +#define AQUILA_SOURCE_H + +#include "source/SignalSource.h" +#include "source/Frame.h" +#include "source/FramesCollection.h" +#include "source/PlainTextFile.h" +#include "source/RawPcmFile.h" +#include "source/WaveFile.h" +#include "source/WaveFileHandler.h" +#include "source/generator/Generator.h" +#include "source/generator/SineGenerator.h" +#include "source/generator/SquareGenerator.h" +#include "source/generator/TriangleGenerator.h" +#include "source/generator/PinkNoiseGenerator.h" +#include "source/generator/WhiteNoiseGenerator.h" +#include "source/window/BarlettWindow.h" +#include "source/window/BlackmanWindow.h" +#include "source/window/FlattopWindow.h" +#include "source/window/GaussianWindow.h" +#include "source/window/HammingWindow.h" +#include "source/window/HannWindow.h" +#include "source/window/RectangularWindow.h" + +#endif // AQUILA_SOURCE_H diff --git a/samplebrain/src/aquila/source/Frame.cpp b/samplebrain/src/aquila/source/Frame.cpp new file mode 100644 index 0000000..8aeb341 --- /dev/null +++ b/samplebrain/src/aquila/source/Frame.cpp @@ -0,0 +1,75 @@ +/** + * @file Frame.cpp + * + * Handling signal frames. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 0.2.2 + */ + +#include "Frame.h" + +namespace Aquila +{ + /** + * Creates the frame object - sets signal source and frame boundaries. + * + * Frame should not change original data, so the source is a const + * reference. + * + * @param source const reference to signal source + * @param indexBegin position of first sample of this frame in the source + * @param indexEnd position of last sample of this frame in the source + */ + Frame::Frame(const SignalSource& source, unsigned int indexBegin, + unsigned int indexEnd): + SignalSource(source.getSampleFrequency()), + m_source(&source), m_begin(indexBegin), + m_end((indexEnd > source.getSamplesCount()) ? source.getSamplesCount() : indexEnd) + { + } + + /** + * Copy constructor. + * + * @param other reference to another frame + */ + Frame::Frame(const Frame &other): + SignalSource(other.m_sampleFrequency), + m_source(other.m_source), m_begin(other.m_begin), m_end(other.m_end) + { + } + + /** + * Move constructor. + * + * @param other rvalue reference to another frame + */ + Frame::Frame(Frame&& other): + SignalSource(other.m_sampleFrequency), + m_source(other.m_source), m_begin(other.m_begin), m_end(other.m_end) + { + } + + /** + * Assignes another frame to this one using copy-and-swap idiom. + * + * @param other reference to another frame + * @return reference to the current object + */ + Frame& Frame::operator=(const Frame& other) + { + Frame temp(other); + swap(temp); + + return *this; + } +} diff --git a/samplebrain/src/aquila/source/Frame.h b/samplebrain/src/aquila/source/Frame.h new file mode 100644 index 0000000..db57f7a --- /dev/null +++ b/samplebrain/src/aquila/source/Frame.h @@ -0,0 +1,127 @@ +/** + * @file Frame.h + * + * Handling signal frames. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 0.2.2 + */ + +#ifndef FRAME_H +#define FRAME_H + +#include "../global.h" +#include "SignalSource.h" +#include +#include + +namespace Aquila +{ + /** + * An ecapsulation of a single frame of the signal. + * + * The Frame class wraps a signal frame (short fragment of a signal). + * Frame itself can be considered as a signal source, being a "window" + * over original signal data. It is a lightweight object which can be + * copied by value. No data are copied - only the pointer to source + * and frame boundaries. + * + * Frame samples are accessed by STL-compatible iterators, as is the + * case with all SignalSource-derived classes. Frame sample number N + * is the same as sample number FRAME_BEGIN+N in the original source. + * + * Example (source size = N, frame size = M, frame starts at 8th sample): + * + * @verbatim + * sample number: 0 8 N + * original source: xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx + * frame: |xxxxxxxxxxxx| + * sample number in frame: 0 M @endverbatim + */ + class AQUILA_EXPORT Frame : public SignalSource + { + public: + Frame(const SignalSource& source, unsigned int indexBegin, + unsigned int indexEnd); + Frame(const Frame& other); + Frame(Frame&& other); + Frame& operator=(const Frame& other); + + /** + * Returns the frame length. + * + * @return frame length as a number of samples + */ + virtual std::size_t getSamplesCount() const + { + return m_end - m_begin; + } + + /** + * Returns number of bits per sample. + * + * @return sample size in bits + */ + virtual unsigned short getBitsPerSample() const + { + return m_source->getBitsPerSample(); + } + + /** + * Gives access to frame samples, indexed from 0 to length()-1. + * + * @param position index of the sample in the frame + * @return sample value + */ + virtual SampleType sample(std::size_t position) const + { + return m_source->sample(m_begin + position); + } + + /** + * Returns sample data (read-only!) as a const C-style array. + * + * Calculates, using C++ pointer arithmetics, where does the frame + * start in the original source, after its convertion to an array. + * + * @return C-style array containing sample data + */ + virtual const SampleType* toArray() const + { + return m_source->toArray() + static_cast(m_begin); + } + + private: + /** + * A non-owning pointer to constant original source (eg. a WAVE file). + */ + const SignalSource* m_source; + + /** + * First and last sample of this frame in the data array/vector. + */ + unsigned int m_begin, m_end; + + /** + * Swaps the frame with another one - exception safe. + * + * @param other reference to swapped frame + */ + void swap(Frame& other) + { + std::swap(m_begin, other.m_begin); + std::swap(m_end, other.m_end); + std::swap(m_source, other.m_source); + } + }; +} + +#endif // FRAME_H diff --git a/samplebrain/src/aquila/source/FramesCollection.cpp b/samplebrain/src/aquila/source/FramesCollection.cpp new file mode 100644 index 0000000..4c61283 --- /dev/null +++ b/samplebrain/src/aquila/source/FramesCollection.cpp @@ -0,0 +1,119 @@ +/** + * @file FramesCollection.cpp + * + * A lightweight wrapper for a vector of Frames. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "FramesCollection.h" +#include "SignalSource.h" + +namespace Aquila +{ + /** + * Creates an empty frames collection. + */ + FramesCollection::FramesCollection(): + m_frames(), m_samplesPerFrame(0) + { + } + + /** + * Creates a collection and creates frames from the source. + * + * @param source a reference to source object + * @param samplesPerFrame how many samples will each frame hold + * @param samplesPerOverlap how many samples are common to adjacent frames + */ + FramesCollection::FramesCollection(const SignalSource& source, + unsigned int samplesPerFrame, + unsigned int samplesPerOverlap): + m_frames(), m_samplesPerFrame(0) + { + divideFrames(source, samplesPerFrame, samplesPerOverlap); + } + + /** + * Destroys the collection, clearing the container. + */ + FramesCollection::~FramesCollection() + { + clear(); + } + + /** + * Creates a collection when duration of each frame is known. + * + * @param source a reference to source object + * @param frameDuration frame duration in milliseconds + * @param overlap overlap as a fraction of frame length (0.0 - 1.0) + */ + FramesCollection FramesCollection::createFromDuration( + const SignalSource &source, double frameDuration, double overlap) + { + unsigned int samplesPerFrame = static_cast( + source.getSampleFrequency() * frameDuration / 1000.0 + ); + unsigned int samplesPerOverlap = static_cast( + samplesPerFrame * overlap + ); + return FramesCollection(source, samplesPerFrame, samplesPerOverlap); + } + + /** + * Performs the actual frame division. + * + * Frames are only "pointing" to the original source. There is no copying + * of sample data. Each frame can be considered as a standalone fragment + * of the source data. + * + * @param source a reference to source object + * @param samplesPerFrame how many samples will each frame hold + * @param samplesPerOverlap how many samples are common to adjacent frames + */ + void FramesCollection::divideFrames(const SignalSource& source, + unsigned int samplesPerFrame, + unsigned int samplesPerOverlap) + { + if (samplesPerFrame == 0) + { + return; + } + m_samplesPerFrame = samplesPerFrame; + const std::size_t sourceSize = source.getSamplesCount(); + const unsigned int nonOverlapped = samplesPerFrame - samplesPerOverlap; + const unsigned int framesCount = sourceSize / nonOverlapped; + + m_frames.reserve(framesCount); + unsigned int indexBegin = 0, indexEnd = 0; + for (std::size_t i = 0; i < framesCount; ++i) + { + // calculate each frame boundaries + // when frame end exceeds source size, break out + indexBegin = i * nonOverlapped; + indexEnd = indexBegin + samplesPerFrame; + if (indexEnd <= sourceSize) + m_frames.push_back(Frame(source, indexBegin, indexEnd)); + else + break; + } + } + + /** + * Deletes all contained frames and clears the collection. + */ + void FramesCollection::clear() + { + m_frames.clear(); + } +} diff --git a/samplebrain/src/aquila/source/FramesCollection.h b/samplebrain/src/aquila/source/FramesCollection.h new file mode 100644 index 0000000..2a72830 --- /dev/null +++ b/samplebrain/src/aquila/source/FramesCollection.h @@ -0,0 +1,182 @@ +/** + * @file FramesCollection.h + * + * A lightweight wrapper for a vector of Frames. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef FRAMESCOLLECTION_H +#define FRAMESCOLLECTION_H + +#include "../global.h" +#include "Frame.h" +#include +#include +#include +#include + +namespace Aquila +{ + class SignalSource; + + /** + * A lightweight wrapper for a vector of Frames. + * + * This class is neccessary to perform signal division into frames, + * which are then stored in container (currently std::vector). The frame + * objects are stored by value as they are very light and cheap to copy. + * + * The reason this wrapper was created is to create some abstraction for + * groups of frames, which can by saved or processed together. For example, + * a spectrogram is an array of spectra calculated individually for each + * frame. Sometimes the calculation doesn't involve the whole signal, + * so a part of it is divided into frames, stored in a FramesCollection + * and then processed. + * + * Individual frame objects can by accessed by iterating over the collection + * using begin() and end() methods. These calls simply return iterators + * pointing to the underlying container. + */ + class AQUILA_EXPORT FramesCollection + { + /** + * Internal storage type. + */ + typedef std::vector Container; + + public: + /** + * An iterator for the collection. + */ + typedef Container::iterator iterator; + + /** + * A const iterator for the collection. + */ + typedef Container::const_iterator const_iterator; + + FramesCollection(); + FramesCollection(const SignalSource& source, + unsigned int samplesPerFrame, + unsigned int samplesPerOverlap = 0); + ~FramesCollection(); + + static FramesCollection createFromDuration(const SignalSource& source, + double frameDuration, + double overlap = 0.0); + + void divideFrames(const SignalSource& source, + unsigned int samplesPerFrame, + unsigned int samplesPerOverlap = 0); + void clear(); + + /** + * Returns number of frames in the collection. + * + * @return frames' container size + */ + std::size_t count() const + { + return m_frames.size(); + } + + /** + * Returns number of samples in each frame. + * + * @return frame size in samples + */ + unsigned int getSamplesPerFrame() const + { + return m_samplesPerFrame; + } + + /** + * Returns nth frame in the collection. + * + * @param index index of the frame in the collection + * @return Frame instance + */ + Frame frame(std::size_t index) const + { + return m_frames[index]; + } + + /** + * Returns an iterator pointing to the first frame. + * + * @return iterator + */ + iterator begin() + { + return m_frames.begin(); + } + + /** + * Returns a const iterator pointing to the first frame. + * + * @return iterator + */ + const_iterator begin() const + { + return m_frames.begin(); + } + + /** + * Returns an iterator pointing one-past-last frame. + * + * @return iterator + */ + iterator end() + { + return m_frames.end(); + } + + /** + * Returns a const iterator pointing one-past-last frame. + * + * @return iterator + */ + const_iterator end() const + { + return m_frames.end(); + } + + /** + * Applies the calculation f to all frames in the collection. + * + * @param f a function whose single argument is a SignalSource + * @return vector of return values of f - one for each frame + */ + template + std::vector apply( + std::function f) const + { + std::vector results; + std::transform(begin(), end(), std::back_inserter(results), f); + return results; + } + + private: + /** + * Frames container. + */ + Container m_frames; + + /** + * Number of samples in each frame. + */ + unsigned int m_samplesPerFrame; + }; +} + +#endif // FRAMESCOLLECTION_H diff --git a/samplebrain/src/aquila/source/PlainTextFile.cpp b/samplebrain/src/aquila/source/PlainTextFile.cpp new file mode 100644 index 0000000..08ae67e --- /dev/null +++ b/samplebrain/src/aquila/source/PlainTextFile.cpp @@ -0,0 +1,59 @@ +/** + * @file PlainTextFile.cpp + * + * Reading samples from a plain text file. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "PlainTextFile.h" +#include +#include +#include + +namespace Aquila +{ + /** + * Creates the data source. + * + * @param filename full path to .txt file + * @param sampleFrequency sample frequency of the data in file + */ + PlainTextFile::PlainTextFile(std::string filename, + FrequencyType sampleFrequency): + SignalSource(sampleFrequency) + { + std::fstream fs; + fs.open(filename.c_str(), std::ios::in); + std::copy(std::istream_iterator(fs), + std::istream_iterator(), + std::back_inserter(m_data)); + fs.close(); + } + + /** + * Saves the given signal source as a plain text file. + * + * @param source source of the data to save + * @param filename destination file + */ + void PlainTextFile::save(const SignalSource& source, + const std::string& filename) + { + std::fstream fs; + fs.open(filename.c_str(), std::ios::out); + std::copy(std::begin(source), + std::end(source), + std::ostream_iterator(fs, "\n")); + fs.close(); + } +} diff --git a/samplebrain/src/aquila/source/PlainTextFile.h b/samplebrain/src/aquila/source/PlainTextFile.h new file mode 100644 index 0000000..1cabb65 --- /dev/null +++ b/samplebrain/src/aquila/source/PlainTextFile.h @@ -0,0 +1,46 @@ +/** + * @file PlainTextFile.h + * + * Reading samples from a plain text file. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef PLAINTEXTFILE_H +#define PLAINTEXTFILE_H + +#include "../global.h" +#include "SignalSource.h" +#include + +namespace Aquila +{ + /** + * Plain text file, where each sample is in new line. + * + * No headers are allowed in the file, only a simple list of numbers + * will work at the moment. + * + * Any numeric type will be converted on the fly to SampleType. Sample + * rate must be known prior to opening the file as the constructor expects + * sample frequency as its second argument. + */ + class AQUILA_EXPORT PlainTextFile : public SignalSource + { + public: + PlainTextFile(std::string filename, FrequencyType sampleFrequency); + + static void save(const SignalSource& source, const std::string& file); + }; +} + +#endif // PLAINTEXTFILE_H diff --git a/samplebrain/src/aquila/source/RawPcmFile.h b/samplebrain/src/aquila/source/RawPcmFile.h new file mode 100644 index 0000000..d679e30 --- /dev/null +++ b/samplebrain/src/aquila/source/RawPcmFile.h @@ -0,0 +1,90 @@ +/** + * @file RawPcmFile.h + * + * Reading raw PCM binary data from file. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef RAWPCMFILE_H +#define RAWPCMFILE_H + +#include "../global.h" +#include "SignalSource.h" +#include +#include +#include +#include + +namespace Aquila +{ + /** + * A class to read raw PCM binary data from file. + + * No headers are allowed in the file. + * + * Any numeric type will be converted on the fly to SampleType. Sample + * rate must be known prior to opening the file as the constructor expects + * sample frequency as its second argument. + */ + template + class AQUILA_EXPORT RawPcmFile : public SignalSource + { + public: + /** + * Creates the data source. + * + * @param filename full path to data file + * @param sampleFrequency sample frequency of the data in file + */ + RawPcmFile(std::string filename, FrequencyType sampleFrequency): + SignalSource(sampleFrequency) + { + std::fstream fs; + fs.open(filename.c_str(), std::ios::in | std::ios::binary); + // get file size by seeking to the end and telling current position + fs.seekg(0, std::ios::end); + std::streamsize fileSize = fs.tellg(); + // seek back to the beginning so read() can access all content + fs.seekg(0, std::ios::beg); + std::size_t samplesCount = fileSize / sizeof(Numeric); + // read raw data into a temporary buffer + Numeric* buffer = new Numeric[samplesCount]; + fs.read((char*)buffer, fileSize); + // copy and implicit conversion to SampleType + m_data.assign(buffer, buffer + samplesCount); + delete [] buffer; + fs.close(); + } + + /** + * Saves the given signal source as a raw PCM file. + * + * @param source source of the data to save + * @param filename destination file + */ + static void save(const SignalSource& source, const std::string& filename) + { + std::fstream fs; + fs.open(filename.c_str(), std::ios::out | std::ios::binary); + std::size_t samplesCount = source.getSamplesCount(); + Numeric* buffer = new Numeric[samplesCount]; + // copy and convert from SampleType to target type + std::copy(std::begin(source), std::end(source), buffer); + fs.write((char*)buffer, samplesCount * sizeof(Numeric)); + delete [] buffer; + fs.close(); + } + }; +} + +#endif // RAWPCMFILE_H diff --git a/samplebrain/src/aquila/source/SignalSource.cpp b/samplebrain/src/aquila/source/SignalSource.cpp new file mode 100644 index 0000000..6c20f51 --- /dev/null +++ b/samplebrain/src/aquila/source/SignalSource.cpp @@ -0,0 +1,243 @@ +/** + * @file SignalSource.cpp + * + * A base signal source class. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "SignalSource.h" +#include +#include +#include +#include + + +namespace Aquila +{ + /** + * Add a constant value to each sample. + * + * @param x value to add + * @return updated source + */ + SignalSource& SignalSource::operator+=(SampleType x) + { + std::transform( + std::begin(m_data), + std::end(m_data), + std::begin(m_data), + [x] (SampleType y) { return x + y; } + ); + return *this; + } + + /** + * Per-sample addition of other signal source. + * + * @param rhs source on the right-hand side of the operator + * @return sum of two sources + */ + SignalSource& SignalSource::operator+=(const SignalSource& rhs) + { + std::transform( + std::begin(m_data), + std::end(m_data), + std::begin(rhs.m_data), + std::begin(m_data), + [] (SampleType x, SampleType y) { return x + y; } + ); + return *this; + } + + /** + * Multiply each sample by a constant value. + * + * @param x multiplier + * @return updated source + */ + SignalSource& SignalSource::operator*=(SampleType x) + { + std::transform( + std::begin(m_data), + std::end(m_data), + std::begin(m_data), + [x] (SampleType y) { return x * y; } + ); + return *this; + } + + /** + * Per-sample multiplication with other signal source. + * + * @param rhs source on the right-hand side of the operator + * @return product of two sources + */ + SignalSource& SignalSource::operator*=(const SignalSource& rhs) + { + std::transform( + std::begin(m_data), + std::end(m_data), + std::begin(rhs.m_data), + std::begin(m_data), + [] (SampleType x, SampleType y) { return x * y; } + ); + return *this; + } + + SignalSource operator+(const SignalSource& lhs, SampleType x) + { + SignalSource result(lhs); + return result += x; + } + + SignalSource operator+(SignalSource&& lhs, SampleType x) + { + lhs += x; + return std::move(lhs); + } + + SignalSource operator+(SampleType x, const SignalSource& rhs) + { + SignalSource result(rhs); + return result += x; + } + + SignalSource operator+(SampleType x, SignalSource&& rhs) + { + rhs += x; + return std::move(rhs); + } + + SignalSource operator+(const SignalSource& lhs, const SignalSource& rhs) + { + SignalSource result(lhs); + return result += rhs; + } + + SignalSource operator+(SignalSource&& lhs, const SignalSource& rhs) + { + lhs += rhs; + return std::move(lhs); + } + + SignalSource operator+(const SignalSource& lhs, SignalSource&& rhs) + { + rhs += lhs; + return std::move(rhs); + } + + SignalSource operator*(const SignalSource& lhs, SampleType x) + { + SignalSource result(lhs); + return result *= x; + } + + SignalSource operator*(SignalSource&& lhs, SampleType x) + { + lhs *= x; + return std::move(lhs); + } + + SignalSource operator*(SampleType x, const SignalSource& rhs) + { + SignalSource result(rhs); + return result *= x; + } + + SignalSource operator*(SampleType x, SignalSource&& rhs) + { + rhs *= x; + return std::move(rhs); + } + + SignalSource operator*(const SignalSource& lhs, const SignalSource& rhs) + { + SignalSource result(lhs); + return result *= rhs; + } + + SignalSource operator*(SignalSource&& lhs, const SignalSource& rhs) + { + lhs *= rhs; + return std::move(lhs); + } + + SignalSource operator*(const SignalSource& lhs, SignalSource&& rhs) + { + rhs *= lhs; + return std::move(rhs); + } + + /** + * Calculates mean value of the signal. + * + * @param source signal source + * @return signal mean + */ + double mean(const SignalSource& source) + { + double sum = std::accumulate(std::begin(source), std::end(source), 0.0); + return sum / source.getSamplesCount(); + } + + /** + * Calculates energy of the signal. + * + * @param source signal source + * @return signal energy + */ + double energy(const SignalSource& source) + { + return std::accumulate( + std::begin(source), + std::end(source), + 0.0, + [] (double acc, SampleType value) { + return acc + value * value; + } + ); + } + + /** + * Calculates power of the signal. + * + * @param source signal source + * @return signal power + */ + double power(const SignalSource& source) + { + return energy(source) / source.getSamplesCount(); + } + + /** + * Calculates Euclidean (L2) norm of the signal. + * + * @param source signal source + * @return norm + */ + double norm(const SignalSource& source) + { + return std::sqrt(energy(source)); + } + + /** + * Calculates root mean square level of the signal. + * + * @param source signal source + * @return RMS level + */ + double rms(const SignalSource& source) + { + return std::sqrt(power(source)); + } +} diff --git a/samplebrain/src/aquila/source/SignalSource.h b/samplebrain/src/aquila/source/SignalSource.h new file mode 100644 index 0000000..3a12230 --- /dev/null +++ b/samplebrain/src/aquila/source/SignalSource.h @@ -0,0 +1,374 @@ +/** + * @file SignalSource.h + * + * A base signal source class. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef SIGNALSOURCE_H +#define SIGNALSOURCE_H + +#include "../global.h" +#include +#include +#include +#include + +namespace Aquila +{ + /** + * An abstraction of any signal source. + * + * This is the base class for all signal sources cooperating with the + * library, be it an array, a text file or a WAVE binary file. Most of + * algorithms defined in the library expect a pointer or a reference + * to SignalSource. The library ships with a few derived + * classes for a quick start, including WaveFile, generators etc. + * + * Signal sources support the concept of iteration. Use + * SignalSource::begin() and SignalSource::end() to obtain iterator objects, + * which allow per-sample data access. The iterators work well with + * C++ standard library algorithms, so feel free to use them instead of + * manually looping and calling SignalSource::sample(). + */ + class AQUILA_EXPORT SignalSource + { + public: + /** + * The default constructor sets sample frequency to 0. + */ + SignalSource(): + m_data(), m_sampleFrequency(0) + { + } + + /** + * CThis overload initializes sample frequency of the source. + * + * @param sampleFrequency sample frequency in Hz + */ + SignalSource(FrequencyType sampleFrequency): + m_data(), m_sampleFrequency(sampleFrequency) + { + } + + /** + * Initialize the source from data sampled at a given frequency. + * + * @param data pointer to a C-style array of samples (numeric values) + * @param dataLength length of the array + * @param sampleFrequency sample frequency in Hz + */ + template + SignalSource(Numeric* data, std::size_t dataLength, + FrequencyType sampleFrequency = 0): + m_data(data, data + dataLength), m_sampleFrequency(sampleFrequency) + { + } + + /** + * Create the source from a vector of samples. + * + * @param data vector of samples + * @param sampleFrequency sample frequency in Hz + */ + SignalSource(const std::vector& data, + FrequencyType sampleFrequency = 0): + m_data(data), m_sampleFrequency(sampleFrequency) + { + } + + /** + * Create the source from a (temporary) vector of samples. + * + * @param data vector of samples + * @param sampleFrequency sample frequency in Hz + */ + SignalSource(std::vector&& data, + FrequencyType sampleFrequency = 0): + m_data(std::move(data)), m_sampleFrequency(sampleFrequency) + { + } + + /** + * The destructor does nothing, but must be defined as virtual. + */ + virtual ~SignalSource() {} + + /** + * Returns sample frequency of the signal. + * + * @return sample frequency in Hz + */ + virtual FrequencyType getSampleFrequency() const + { + return m_sampleFrequency; + } + + /** + * Sets sample frequency of the signal. + * + * @param sampleFrequency sample frequency in Hz + */ + virtual void setSampleFrequency(FrequencyType sampleFrequency) + { + m_sampleFrequency = sampleFrequency; + } + + /** + * Returns number of bits per signal sample. + * + * @return sample size in bits + */ + virtual unsigned short getBitsPerSample() const + { + return 8 * sizeof(SampleType); + } + + /** + * Returns number of samples in the source. + * + * @return samples count + */ + virtual std::size_t getSamplesCount() const + { + return m_data.size(); + } + + /** + * Returns sample located at the "position" in the signal. + * + * @param position sample index in the signal + * @return sample value + */ + virtual SampleType sample(std::size_t position) const + { + return m_data[position]; + } + + /** + * Returns sample data (read-only!) as a const C-style array. + * + * Because vector guarantees to be contiguous in memory, we can + * return the address of the first element in the vector. + * It is valid only until next operation which modifies the vector. + * + * @return C-style array containing sample data + */ + virtual const SampleType* toArray() const + { + return &m_data[0]; + } + + /** + * Returns number of samples in the source. + * + * This method is an alias to getSamplesCount() and it should not be + * implemented in derived classes. + * + * @return samples count + */ + std::size_t length() const + { + return getSamplesCount(); + } + + class iterator; + + /** + * Returns an iterator pointing to the first sample in the source. + * + * @return iterator + */ + iterator begin() const + { + return iterator(this, 0); + } + + /** + * Returns an iterator pointing to the "one past last" sample. + * + * @return iterator + */ + iterator end() const + { + return iterator(this, getSamplesCount()); + } + + /** + * Iterator class enabling sequential data access. + * + * It is a forward iterator with a range from the first sample in the + * source to "one past last" sample. + */ + class AQUILA_EXPORT iterator : + public std::iterator + { + public: + /** + * Creates an iterator associated with a given source. + * + * @param source pointer to a source on which the iterator will work + * @param i index of the sample in the source + */ + explicit iterator(const SignalSource* source = nullptr, unsigned int i = 0): + m_source(source), idx(i) + { + } + + /** + * Assigns another iterator in a memberwise fashion. + * + * @param other right-hand value iterator + * @return reference to self + */ + iterator& operator=(const iterator& other) + { + m_source = other.m_source; + idx = other.idx; + return (*this); + } + + /** + * Compares two iterators for equality. + * + * Iterators are equal only when they belong to the same source + * and point to the same sample in the source. + * + * @param other right-hand value iterator + * @return true, if the iterators are equal + */ + bool operator==(const iterator& other) const + { + return m_source == other.m_source && idx == other.idx; + } + + /** + * Compares two iterators for inequality. + * + * Negates the equality operator. + * + * @param other right-hand value iterator + * @return true only when the iterators are not equal + */ + bool operator!=(const iterator& other) const + { + return !operator==(other); + } + + /** + * Moves the iterator one sample to the right (prefix version). + * + * @return reference to self + */ + iterator& operator++() + { + ++idx; + return (*this); + } + + /** + * Moves the iterator one sample to the right (postfix version). + * + * @return a copy of self before incrementing + */ + iterator operator++(int) + { + iterator tmp(*this); + ++(*this); + return tmp; + } + + /** + * Dereferences the iterator. + * + * @return signal sample value. + */ + SampleType operator*() const + { + return m_source->sample(idx); + } + + /** + * Returns the distance from the beginning of the source. + * + * @return number of samples between beginning and current position + */ + std::size_t getPosition() const + { + return idx; + } + + private: + /** + * Signal source - as a pointer - the iterators must be copyable. + */ + const SignalSource* m_source; + + /** + * Iterator's position in the source. + */ + std::size_t idx; + }; + + SignalSource& operator+=(SampleType x); + SignalSource& operator+=(const SignalSource& rhs); + SignalSource& operator*=(SampleType x); + SignalSource& operator*=(const SignalSource& rhs); + + protected: + /** + * Actual sample data. + */ + std::vector m_data; + + /** + * Sample frequency of the data. + */ + FrequencyType m_sampleFrequency; + }; + + SignalSource operator+(const SignalSource& lhs, SampleType x); + SignalSource operator+(SignalSource&& lhs, SampleType x); + SignalSource operator+(SampleType x, const SignalSource& rhs); + SignalSource operator+(SampleType x, SignalSource&& rhs); + SignalSource operator+(const SignalSource& lhs, const SignalSource& rhs); + SignalSource operator+(SignalSource&& lhs, const SignalSource& rhs); + SignalSource operator+(const SignalSource& lhs, SignalSource&& rhs); + + SignalSource operator*(const SignalSource& lhs, SampleType x); + SignalSource operator*(SignalSource&& lhs, SampleType x); + SignalSource operator*(SampleType x, const SignalSource& rhs); + SignalSource operator*(SampleType x, SignalSource&& rhs); + SignalSource operator*(const SignalSource& lhs, const SignalSource& rhs); + SignalSource operator*(SignalSource&& lhs, const SignalSource& rhs); + SignalSource operator*(const SignalSource& lhs, SignalSource&& rhs); + + /*************************************************************************** + * + * Free-standing functions closely related to signals. + * + **************************************************************************/ + + double mean(const SignalSource& source); + + double energy(const SignalSource& source); + + double power(const SignalSource& source); + + double norm(const SignalSource& source); + + double rms(const SignalSource& source); +} + +#endif // SIGNALSOURCE_H diff --git a/samplebrain/src/aquila/source/WaveFile.cpp b/samplebrain/src/aquila/source/WaveFile.cpp new file mode 100644 index 0000000..a12c7a2 --- /dev/null +++ b/samplebrain/src/aquila/source/WaveFile.cpp @@ -0,0 +1,93 @@ +/** + * @file WaveFile.cpp + * + * WAVE file handling as a signal source. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 0.0.7 + */ + +#include "WaveFile.h" +#include "WaveFileHandler.h" + +namespace Aquila +{ + /** + * Creates a wave file object and immediately loads data from file. + * + * @param filename full path to .wav file + * @param channel LEFT or RIGHT (the default setting is LEFT) + */ + WaveFile::WaveFile(const std::string& filename, StereoChannel channel): + SignalSource(), m_filename(filename) + { + load(m_filename, channel); + } + + /** + * Deletes the WaveFile object. + */ + WaveFile::~WaveFile() + { + } + + /** + * Reads the header and channel data from given .wav file. + * + * Data are read into a temporary buffer and then converted to + * channel sample vectors. If source is a mono recording, samples + * are written to left channel. + * + * To improve performance, no format checking is performed. + * + * @param filename full path to .wav file + * @param channel which audio channel to read (for formats other than mono) + */ + void WaveFile::load(const std::string& filename, StereoChannel channel) + { + m_filename = filename; + m_data.clear(); + ChannelType dummy; + WaveFileHandler handler(m_filename); + if (LEFT == channel) + { + handler.readHeaderAndChannels(m_header, m_data, dummy); + } + else + { + handler.readHeaderAndChannels(m_header, dummy, m_data); + } + m_sampleFrequency = m_header.SampFreq; + } + + /** + * Saves the given signal source as a .wav file. + * + * @param source source of the data to save + * @param filename destination file + */ + void WaveFile::save(const SignalSource& source, const std::string& filename) + { + WaveFileHandler handler(filename); + handler.save(source); + } + + /** + * Returns the audio recording length + * + * @return recording length in milliseconds + */ + unsigned int WaveFile::getAudioLength() const + { + return static_cast(m_header.WaveSize / + static_cast(m_header.BytesPerSec) * 1000); + } +} diff --git a/samplebrain/src/aquila/source/WaveFile.h b/samplebrain/src/aquila/source/WaveFile.h new file mode 100644 index 0000000..67eb195 --- /dev/null +++ b/samplebrain/src/aquila/source/WaveFile.h @@ -0,0 +1,186 @@ +/** + * @file WaveFile.h + * + * WAVE file handling as a signal source. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 0.0.7 + */ + +#ifndef WAVEFILE_H +#define WAVEFILE_H + +#include "../global.h" +#include "SignalSource.h" +#include +#include +#include +#include + +namespace Aquila +{ + /** + * Which channel to use when reading stereo recordings. + */ + enum StereoChannel { LEFT, RIGHT }; + + /** + * .wav file header structure. + */ + struct WaveHeader + { + char RIFF[4]; + std::uint32_t DataLength; + char WAVE[4]; + char fmt_[4]; + std::uint32_t SubBlockLength; + std::uint16_t formatTag; + std::uint16_t Channels; + std::uint32_t SampFreq; + std::uint32_t BytesPerSec; + std::uint16_t BytesPerSamp; + std::uint16_t BitsPerSamp; + char data[4]; + std::uint32_t WaveSize; + }; + + /** + * Wave file data access. + * + * Binary files in WAVE format (.wav extension) can serve as data input for + * Aquila. With this class, you can read the metadata and the actual + * waveform data from the file. The supported formats are: + * + * - 8-bit mono + * - 8-bit stereo* + * - 16-bit mono + * - 16-bit stereo* + * + * For stereo data, only only one of the channels is loaded from file. + * By default this is the left channel, but you can control this from the + * constructor parameter. + * + * There are no requirements for sample frequency of the data. + */ + class AQUILA_EXPORT WaveFile : public SignalSource + { + public: + /** + * Audio channel representation. + */ + typedef decltype(m_data) ChannelType; + + explicit WaveFile(const std::string& filename, + StereoChannel channel = LEFT); + ~WaveFile(); + + void load(const std::string& filename, StereoChannel channel); + static void save(const SignalSource& source, const std::string& file); + + /** + * Returns the filename. + * + * @return full path to currently loaded file + */ + std::string getFilename() const + { + return m_filename; + } + + /** + * Returns number of channels. + * + * @return 1 for mono, 2 for stereo, other types are not recognized + */ + unsigned short getChannelsNum() const + { + return m_header.Channels; + } + + /** + * Checks if this is a mono recording. + * + * @return true if there is only one channel + */ + bool isMono() const + { + return 1 == getChannelsNum(); + } + + /** + * Checks if this is a stereo recording. + * + * @return true if there are two channels + */ + bool isStereo() const + { + return 2 == getChannelsNum(); + } + + /** + * Returns the number of bytes per second. + * + * @return product of sample frequency and bytes per sample + */ + unsigned int getBytesPerSec() const + { + return m_header.BytesPerSec; + } + + /** + * Returns number of bytes per sample. + * + * @return 1 for 8b-mono, 2 for 8b-stereo or 16b-mono, 4 for 16b-stereo + */ + unsigned int getBytesPerSample() const + { + return m_header.BytesPerSamp; + } + + /** + * Returns number of bits per sample + * + * @return 8 or 16 + */ + virtual unsigned short getBitsPerSample() const + { + return m_header.BitsPerSamp; + } + + /** + * Returns the recording size (without header). + * + * The return value is a raw byte count. To know the real sample count, + * it must be divided by bytes per sample. + * + * @return byte count + */ + unsigned int getWaveSize() const + { + return m_header.WaveSize; + } + + unsigned int getAudioLength() const; + + private: + /** + * Full path of the .wav file. + */ + std::string m_filename; + + /** + * Header structure. + */ + WaveHeader m_header; + }; +} + +#endif // WAVEFILE_H diff --git a/samplebrain/src/aquila/source/WaveFileHandler.cpp b/samplebrain/src/aquila/source/WaveFileHandler.cpp new file mode 100644 index 0000000..0fe6fb1 --- /dev/null +++ b/samplebrain/src/aquila/source/WaveFileHandler.cpp @@ -0,0 +1,269 @@ +/** + * @file WaveFileHandler.cpp + * + * A utility class to handle loading and saving of .wav files. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "WaveFileHandler.h" +#include "WaveFile.h" +#include +#include +#include + +namespace Aquila +{ + /** + * Create the handler and tell it which file to read later. + * + * @param filename .wav file name + */ + WaveFileHandler::WaveFileHandler(const std::string& filename): + m_filename(filename) + { + } + + /** + * Reads WAVE header and audio channel data from file. + * + * @param header reference to header instance which will be filled + * @param leftChannel reference to left audio channel + * @param rightChannel reference to right audio channel + */ + void WaveFileHandler::readHeaderAndChannels(WaveHeader &header, + WaveFile::ChannelType& leftChannel, WaveFile::ChannelType& rightChannel) + { + // first we read header from the stream + // then as we know now the data size, we create a temporary + // buffer and read raw data into that buffer + std::fstream fs; + fs.open(m_filename.c_str(), std::ios::in | std::ios::binary); + fs.read((char*)(&header), sizeof(WaveHeader)); + short* data = new short[header.WaveSize/2]; + fs.read((char*)data, header.WaveSize); + fs.close(); + + // initialize data channels (using right channel only in stereo mode) + unsigned int channelSize = header.WaveSize/header.BytesPerSamp; + leftChannel.resize(channelSize); + if (2 == header.Channels) + rightChannel.resize(channelSize); + + // most important conversion happens right here + if (16 == header.BitsPerSamp) + { + if (2 == header.Channels) + decode16bitStereo(leftChannel, rightChannel, data, channelSize); + else + decode16bit(leftChannel, data, channelSize); + } + else + { + if (2 == header.Channels) + decode8bitStereo(leftChannel, rightChannel, data, channelSize); + else + decode8bit(leftChannel, data, channelSize); + } + + // clear the buffer + delete [] data; + } + + /** + * Saves the given signal source as a .wav file. + * + * @param source source of the data to save + */ + void WaveFileHandler::save(const SignalSource& source) + { + WaveHeader header; + createHeader(source, header); + std::ofstream fs; + fs.open(m_filename.c_str(), std::ios::out | std::ios::binary); + fs.write((const char*)(&header), sizeof(WaveHeader)); + + std::size_t waveSize = header.WaveSize; + short* data = new short[waveSize/2]; + if (16 == header.BitsPerSamp) + { + encode16bit(source, data, waveSize/2); + } + else + { + encode8bit(source, data, waveSize/2); + } + fs.write((const char*)data, waveSize); + + delete [] data; + fs.close(); + } + + + /** + * Populates a .wav file header with values obtained from the source. + * + * @param source data source + * @param header the header which will be filled with correct parameters + */ + void WaveFileHandler::createHeader(const SignalSource &source, WaveHeader &header) + { + std::uint32_t frequency = static_cast(source.getSampleFrequency()); + // saving only mono files at the moment + std::uint16_t channels = 1; + std::uint16_t bitsPerSample = source.getBitsPerSample(); + // higher dynamic sources will be converted down to 16 bits per sample + if (bitsPerSample > 16) + bitsPerSample = 16; + std::uint32_t bytesPerSec = frequency * channels * bitsPerSample / 8; + std::uint32_t waveSize = source.getSamplesCount() * channels * bitsPerSample / 8; + + strncpy(header.RIFF, "RIFF", 4); + // DataLength is the file size excluding first two header fields - + // - RIFF and DataLength itself, which together take 8 bytes to store + header.DataLength = waveSize + sizeof(WaveHeader) - 8; + strncpy(header.WAVE, "WAVE", 4); + strncpy(header.fmt_, "fmt ", 4); + header.SubBlockLength = 16; + header.formatTag = 1; + header.Channels = channels; + header.SampFreq = frequency; + header.BytesPerSec = bytesPerSec; + header.BytesPerSamp = header.Channels * bitsPerSample / 8; + header.BitsPerSamp = bitsPerSample; + strncpy(header.data, "data", 4); + header.WaveSize = waveSize; + } + + /** + * Decodes 16 bit mono data into a suitable audio channel format. + * + * @param channel a reference to audio channel + * @param data raw data buffer + * @param channelSize expected number of samples in channel + */ + void WaveFileHandler::decode16bit(WaveFile::ChannelType& channel, short* data, std::size_t channelSize) + { + for (std::size_t i = 0; i < channelSize; ++i) + { + channel[i] = data[i]; + } + } + + /** + * Decodes 16 bit stereo data into two audio channels. + * + * @param leftChannel a reference to left audio channel + * @param rightChannel a reference to right audio channel + * @param data raw data buffer + * @param channelSize expected number of samples (same for both channels) + */ + void WaveFileHandler::decode16bitStereo(WaveFile::ChannelType& leftChannel, + WaveFile::ChannelType& rightChannel, short* data, std::size_t channelSize) + { + for (std::size_t i = 0; i < channelSize; ++i) + { + leftChannel[i] = data[2*i]; + rightChannel[i] = data[2*i+1]; + } + } + + /** + * Decodes 8 bit mono data into a suitable audio channel format. + * + * @param channel a reference to audio channel + * @param data raw data buffer + * @param channelSize expected number of samples in channel + */ + void WaveFileHandler::decode8bit(WaveFile::ChannelType& channel, short* data, std::size_t channelSize) + { + // low byte and high byte of a 16b word + unsigned char lb, hb; + for (std::size_t i = 0; i < channelSize; ++i) + { + splitBytes(data[i / 2], lb, hb); + // only one channel collects samples + channel[i] = lb - 128; + } + } + + /** + * Decodes 8 bit stereo data into two audio channels. + * + * @param leftChannel a reference to left audio channel + * @param rightChannel a reference to right audio channel + * @param data raw data buffer + * @param channelSize expected number of samples (same for both channels) + */ + void WaveFileHandler::decode8bitStereo(WaveFile::ChannelType& leftChannel, + WaveFile::ChannelType& rightChannel, short* data, std::size_t channelSize) + { + // low byte and high byte of a 16b word + unsigned char lb, hb; + for (std::size_t i = 0; i < channelSize; ++i) + { + splitBytes(data[i / 2], lb, hb); + // left channel is in low byte, right in high + // values are unipolar, so we move them by half + // of the dynamic range + leftChannel[i] = lb - 128; + rightChannel[i] = hb - 128; + } + } + + /** + * Encodes the source data as an array of 16-bit values. + * + * @param source original signal source + * @param data the data buffer to be written + * @param dataSize size of the buffer + */ + void WaveFileHandler::encode16bit(const SignalSource& source, short* data, std::size_t dataSize) + { + for (std::size_t i = 0; i < dataSize; ++i) + { + short sample = static_cast(source.sample(i)); + data[i] = sample; + } + } + + /** + * Encodes the source data as an array of 8-bit values stored in shorts. + * + * @param source original signal source + * @param data the data buffer to be written + * @param dataSize size of the buffer + */ + void WaveFileHandler::encode8bit(const SignalSource& source, short* data, std::size_t dataSize) + { + for (std::size_t i = 0; i < dataSize; ++i) + { + unsigned char sample1 = static_cast(source.sample(2 * i) + 128); + unsigned char sample2 = static_cast(source.sample(2 * i + 1) + 128); + short hb = sample1, lb = sample2; + data[i] = ((hb << 8) & 0xFF00) | (lb & 0x00FF); + } + } + + /** + * Splits a 16-b number to lower and upper byte. + * + * @param twoBytes number to split + * @param lb lower byte (by reference) + * @param hb upper byte (by reference) + */ + void WaveFileHandler::splitBytes(short twoBytes, unsigned char& lb, unsigned char& hb) + { + lb = twoBytes & 0x00FF; + hb = (twoBytes >> 8) & 0x00FF; + } +} diff --git a/samplebrain/src/aquila/source/WaveFileHandler.h b/samplebrain/src/aquila/source/WaveFileHandler.h new file mode 100644 index 0000000..4fadd1a --- /dev/null +++ b/samplebrain/src/aquila/source/WaveFileHandler.h @@ -0,0 +1,71 @@ +/** + * @file WaveFileHandler.h + * + * A utility class to handle loading and saving of .wav files. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef WAVEFILEHANDLER_H +#define WAVEFILEHANDLER_H + +#include "../global.h" +#include "SignalSource.h" +#include "WaveFile.h" +#include +#include + +namespace Aquila +{ + /** + * Forward declaration to avoid including WaveFile.h header. + */ + struct WaveHeader; + + /** + * A utility class to handle loading and saving of .wav files. + */ + class AQUILA_EXPORT WaveFileHandler + { + public: + WaveFileHandler(const std::string& filename); + + void readHeaderAndChannels(WaveHeader& header, + WaveFile::ChannelType& leftChannel, WaveFile::ChannelType& rightChannel); + + void save(const SignalSource& source); + + static void decode16bit(WaveFile::ChannelType& channel, + short* data, std::size_t channelSize); + static void decode16bitStereo(WaveFile::ChannelType& leftChannel, + WaveFile::ChannelType& rightChannel, short* data, std::size_t channelSize); + + static void decode8bit(WaveFile::ChannelType& channel, + short* data, std::size_t channelSize); + static void decode8bitStereo(WaveFile::ChannelType& leftChannel, + WaveFile::ChannelType& rightChannel, short* data, std::size_t channelSize); + + static void encode16bit(const SignalSource& source, short* data, std::size_t dataSize); + static void encode8bit(const SignalSource& source, short* data, std::size_t dataSize); + + private: + void createHeader(const SignalSource& source, WaveHeader& header); + static void splitBytes(short twoBytes, unsigned char& lb, unsigned char& hb); + + /** + * Destination or source file. + */ + std::string m_filename; + }; +} + +#endif // WAVEFILEHANDLER_H diff --git a/samplebrain/src/aquila/source/generator/Generator.cpp b/samplebrain/src/aquila/source/generator/Generator.cpp new file mode 100644 index 0000000..bd550c0 --- /dev/null +++ b/samplebrain/src/aquila/source/generator/Generator.cpp @@ -0,0 +1,32 @@ +/** + * @file Generator.cpp + * + * An interface for signal generators. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "Generator.h" + +namespace Aquila +{ + /** + * Creates the generator object. + * + * @param sampleFrequency sample frequency of the signal + */ + Generator::Generator(FrequencyType sampleFrequency): + SignalSource(sampleFrequency), m_frequency(0), m_amplitude(0), + m_phase(0.0) + { + } +} diff --git a/samplebrain/src/aquila/source/generator/Generator.h b/samplebrain/src/aquila/source/generator/Generator.h new file mode 100644 index 0000000..897325a --- /dev/null +++ b/samplebrain/src/aquila/source/generator/Generator.h @@ -0,0 +1,104 @@ +/** + * @file Generator.h + * + * An interface for signal generators. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef GENERATOR_H +#define GENERATOR_H + +#include "../SignalSource.h" +#include "../../global.h" +#include + +namespace Aquila +{ + /** + * The base interface for signal generators. + */ + class AQUILA_EXPORT Generator : public SignalSource + { + public: + /** + * Creates the generator object. + * + * @param sampleFrequency sample frequency of the data in array + */ + Generator(FrequencyType sampleFrequency); + + /** + * Sets frequency of the generated signal. + * + * @param frequency signal frequency + * @return the current object for fluent interface + */ + Generator& setFrequency(FrequencyType frequency) + { + m_frequency = frequency; + + return *this; + } + + /** + * Sets amplitude of the generated signal. + * + * @param amplitude signal amplitude + * @return the current object for fluent interface + */ + Generator& setAmplitude(SampleType amplitude) + { + m_amplitude = amplitude; + + return *this; + } + + /** + * Sets phase shift of the generated wave. + * + * @param phase phase shift (0 < phase <= 1) + * @return the current object for fluent interface + */ + Generator& setPhase(double phase) + { + m_phase = phase; + + return *this; + } + + /** + * Generates a given number of samples. + * + * @param samplesCount how many samples to generate + */ + virtual void generate(std::size_t samplesCount) = 0; + + protected: + /** + * Frequency of the generated signal (not always used). + */ + FrequencyType m_frequency; + + /** + * Amplitude of the generated signal (not always used). + */ + SampleType m_amplitude; + + /** + * Phase shift as a fraction of whole period (default = 0.0). + */ + double m_phase; + }; +} + +#endif // GENERATOR_H diff --git a/samplebrain/src/aquila/source/generator/PinkNoiseGenerator.cpp b/samplebrain/src/aquila/source/generator/PinkNoiseGenerator.cpp new file mode 100644 index 0000000..51e0b28 --- /dev/null +++ b/samplebrain/src/aquila/source/generator/PinkNoiseGenerator.cpp @@ -0,0 +1,80 @@ +/** + * @file PinkNoiseGenerator.cpp + * + * Pink noise generator. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "PinkNoiseGenerator.h" +#include "../../functions.h" + +namespace Aquila +{ + /** + * Creates the generator object. + * + * @param sampleFrequency sample frequency of the signal + */ + PinkNoiseGenerator::PinkNoiseGenerator(FrequencyType sampleFrequency): + Generator(sampleFrequency), key(0), maxKey(0xFFFF) + { + } + + /** + * Fills the buffer with generated pink noise samples. + * + * @param samplesCount how many samples to generate + */ + void PinkNoiseGenerator::generate(std::size_t samplesCount) + { + m_data.resize(samplesCount); + + // Voss algorithm initialization + maxKey = 0xFFFF; + key = 0; + for (std::size_t i = 0; i < whiteSamplesNum; ++i) + { + whiteSamples[i] = randomDouble() - 0.5; + } + + for (std::size_t i = 0; i < samplesCount; ++i) + { + m_data[i] = m_amplitude * pinkSample(); + } + } + + /** + * Generates a single pink noise sample using Voss algorithm. + * + * @return pink noise sample + */ + double PinkNoiseGenerator::pinkSample() + { + int lastKey = key; + double sum = 0; + + key++; + if (key > maxKey) + key = 0; + + int diff = lastKey ^ key; + for (std::size_t i = 0; i < whiteSamplesNum; ++i) + { + if (diff & (1 << i)) + whiteSamples[i] = randomDouble() - 0.5; + sum += whiteSamples[i]; + } + + return sum / whiteSamplesNum; + } +} diff --git a/samplebrain/src/aquila/source/generator/PinkNoiseGenerator.h b/samplebrain/src/aquila/source/generator/PinkNoiseGenerator.h new file mode 100644 index 0000000..cb087b6 --- /dev/null +++ b/samplebrain/src/aquila/source/generator/PinkNoiseGenerator.h @@ -0,0 +1,60 @@ +/** + * @file PinkNoiseGenerator.h + * + * Pink noise generator. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef PINKNOISEGENERATOR_H +#define PINKNOISEGENERATOR_H + +#include "Generator.h" + +namespace Aquila +{ + /** + * Pink noise generator using Voss algorithm. + */ + class AQUILA_EXPORT PinkNoiseGenerator : public Generator + { + public: + PinkNoiseGenerator(FrequencyType sampleFrequency); + + virtual void generate(std::size_t samplesCount); + + private: + double pinkSample(); + + /** + * Number of white noise samples to use in Voss algorithm. + */ + enum { whiteSamplesNum = 20 }; + + /** + * An internal buffer for white noise samples. + */ + double whiteSamples[whiteSamplesNum]; + + /** + * A key marking which of the white noise samples should change. + */ + int key; + + /** + * Maximum key value. + */ + int maxKey; + }; +} + +#endif // PINKNOISEGENERATOR_H diff --git a/samplebrain/src/aquila/source/generator/SineGenerator.cpp b/samplebrain/src/aquila/source/generator/SineGenerator.cpp new file mode 100644 index 0000000..87c540b --- /dev/null +++ b/samplebrain/src/aquila/source/generator/SineGenerator.cpp @@ -0,0 +1,51 @@ +/** + * @file SineGenerator.cpp + * + * Sine wave generator. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "SineGenerator.h" +#include + +namespace Aquila +{ + /** + * Creates the generator object. + * + * @param sampleFrequency sample frequency of the signal + */ + SineGenerator::SineGenerator(FrequencyType sampleFrequency): + Generator(sampleFrequency) + { + } + + /** + * Fills the buffer with generated sine samples. + * + * @param samplesCount how many samples to generate + */ + void SineGenerator::generate(std::size_t samplesCount) + { + m_data.resize(samplesCount); + double normalizedFrequency = m_frequency / + static_cast(m_sampleFrequency); + for (std::size_t i = 0; i < samplesCount; ++i) + { + m_data[i] = m_amplitude * std::sin( + 2.0 * M_PI * normalizedFrequency * i + + m_phase * 2.0 * M_PI + ); + } + } +} diff --git a/samplebrain/src/aquila/source/generator/SineGenerator.h b/samplebrain/src/aquila/source/generator/SineGenerator.h new file mode 100644 index 0000000..2bf6d81 --- /dev/null +++ b/samplebrain/src/aquila/source/generator/SineGenerator.h @@ -0,0 +1,37 @@ +/** + * @file SineGenerator.h + * + * Sine wave generator. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef SINEGENERATOR_H +#define SINEGENERATOR_H + +#include "Generator.h" + +namespace Aquila +{ + /** + * Sine wave generator. + */ + class AQUILA_EXPORT SineGenerator : public Generator + { + public: + SineGenerator(FrequencyType sampleFrequency); + + virtual void generate(std::size_t samplesCount); + }; +} + +#endif // SINEGENERATOR_H diff --git a/samplebrain/src/aquila/source/generator/SquareGenerator.cpp b/samplebrain/src/aquila/source/generator/SquareGenerator.cpp new file mode 100644 index 0000000..9f82ce3 --- /dev/null +++ b/samplebrain/src/aquila/source/generator/SquareGenerator.cpp @@ -0,0 +1,52 @@ +/** + * @file SquareGenerator.cpp + * + * Square wave generator. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "SquareGenerator.h" +#include + +namespace Aquila +{ + /** + * Creates the generator object. + * + * @param sampleFrequency sample frequency of the signal + */ + SquareGenerator::SquareGenerator(FrequencyType sampleFrequency): + Generator(sampleFrequency), m_duty(0.5) + { + } + + /** + * Fills the buffer with generated square samples. + *` + * @param samplesCount how many samples to generate + */ + void SquareGenerator::generate(std::size_t samplesCount) + { + m_data.resize(samplesCount); + std::size_t samplesPerPeriod = static_cast( + m_sampleFrequency / static_cast(m_frequency)); + std::size_t positiveLength = static_cast(m_duty * + samplesPerPeriod); + + for (std::size_t i = 0; i < samplesCount; ++i) + { + std::size_t t = i % samplesPerPeriod; + m_data[i] = m_amplitude * (t < positiveLength ? 1 : -1); + } + } +} diff --git a/samplebrain/src/aquila/source/generator/SquareGenerator.h b/samplebrain/src/aquila/source/generator/SquareGenerator.h new file mode 100644 index 0000000..077842a --- /dev/null +++ b/samplebrain/src/aquila/source/generator/SquareGenerator.h @@ -0,0 +1,58 @@ +/** + * @file SquareGenerator.h + * + * Square wave generator. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef SQUAREGENERATOR_H +#define SQUAREGENERATOR_H + +#include "Generator.h" + +namespace Aquila +{ + /** + * Square wave generator. + */ + class AQUILA_EXPORT SquareGenerator : public Generator + { + public: + SquareGenerator(FrequencyType sampleFrequency); + + virtual void generate(std::size_t samplesCount); + + /** + * Sets duty cycle of the generated square wave. + * + * Duty cycle is a fraction of period in which the signal is positive. + * + * @param duty duty cycle (0 < duty <= 1) + * @return the current object for fluent interface + */ + SquareGenerator& setDuty(double duty) + { + m_duty = duty; + + return *this; + } + + private: + /** + * Duty cycle, default = 0.5. + */ + double m_duty; + }; +} + +#endif // SQUAREGENERATOR_H diff --git a/samplebrain/src/aquila/source/generator/TriangleGenerator.cpp b/samplebrain/src/aquila/source/generator/TriangleGenerator.cpp new file mode 100644 index 0000000..17c1ff7 --- /dev/null +++ b/samplebrain/src/aquila/source/generator/TriangleGenerator.cpp @@ -0,0 +1,71 @@ +/** + * @file TriangleGenerator.cpp + * + * Triangle (and sawtooth) wave generator. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "TriangleGenerator.h" +#include + +namespace Aquila +{ + /** + * Creates the generator object. + * + * @param sampleFrequency sample frequency of the signal + */ + TriangleGenerator::TriangleGenerator(FrequencyType sampleFrequency): + Generator(sampleFrequency), m_width(1.0) + { + } + + /** + * Fills the buffer with generated triangle wave samples. + * + * The default behaviour is to generate a sawtooth wave. To change that + * to a triangle wave, call setWidth() with some value between 0 and 1. + * + * @param samplesCount how many samples to generate + */ + void TriangleGenerator::generate(std::size_t samplesCount) + { + m_data.resize(samplesCount); + + double dt = 1.0 / m_sampleFrequency, period = 1.0 / m_frequency; + double risingLength = m_width * period; + double fallingLength = period - risingLength; + double risingIncrement = + (risingLength != 0) ? (2.0 * m_amplitude / risingLength) : 0; + double fallingDecrement = + (fallingLength != 0) ? (2.0 * m_amplitude / fallingLength) : 0; + + double t = 0; + for (std::size_t i = 0; i < samplesCount; ++i) + { + if (t > period) + { + t -= period; + } + if (t < risingLength) + { + m_data[i] = -m_amplitude + t * risingIncrement; + } + else + { + m_data[i] = m_amplitude - (t - risingLength) * fallingDecrement; + } + t += dt; + } + } +} diff --git a/samplebrain/src/aquila/source/generator/TriangleGenerator.h b/samplebrain/src/aquila/source/generator/TriangleGenerator.h new file mode 100644 index 0000000..8382e0a --- /dev/null +++ b/samplebrain/src/aquila/source/generator/TriangleGenerator.h @@ -0,0 +1,58 @@ +/** + * @file TriangleGenerator.h + * + * Triangle (and sawtooth) wave generator. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef TRIANGLEGENERATOR_H +#define TRIANGLEGENERATOR_H + +#include "Generator.h" + +namespace Aquila +{ + /** + * Triangle (and sawtooth) wave generator. + */ + class AQUILA_EXPORT TriangleGenerator : public Generator + { + public: + TriangleGenerator(FrequencyType sampleFrequency); + + virtual void generate(std::size_t samplesCount); + + /** + * Sets slope width of the generated triangle wave. + * + * Slope width is a fraction of period in which signal is rising. + * + * @param width slope width (0 < width <= 1) + * @return the current object for fluent interface + */ + TriangleGenerator& setWidth(double width) + { + m_width = width; + + return *this; + } + + private: + /** + * Slope width, default = 1.0 (generates sawtooth wave). + */ + double m_width; + }; +} + +#endif // TRIANGLEGENERATOR_H diff --git a/samplebrain/src/aquila/source/generator/WhiteNoiseGenerator.cpp b/samplebrain/src/aquila/source/generator/WhiteNoiseGenerator.cpp new file mode 100644 index 0000000..bcce665 --- /dev/null +++ b/samplebrain/src/aquila/source/generator/WhiteNoiseGenerator.cpp @@ -0,0 +1,46 @@ +/** + * @file WhiteNoiseGenerator.cpp + * + * White noise generator. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "WhiteNoiseGenerator.h" +#include "../../functions.h" + +namespace Aquila +{ + /** + * Creates the generator object. + * + * @param sampleFrequency sample frequency of the signal + */ + WhiteNoiseGenerator::WhiteNoiseGenerator(FrequencyType sampleFrequency): + Generator(sampleFrequency) + { + } + + /** + * Fills the buffer with generated white noise samples. + * + * @param samplesCount how many samples to generate + */ + void WhiteNoiseGenerator::generate(std::size_t samplesCount) + { + m_data.resize(samplesCount); + for (std::size_t i = 0; i < samplesCount; ++i) + { + m_data[i] = m_amplitude * (randomDouble() - 0.5); + } + } +} diff --git a/samplebrain/src/aquila/source/generator/WhiteNoiseGenerator.h b/samplebrain/src/aquila/source/generator/WhiteNoiseGenerator.h new file mode 100644 index 0000000..04d15e8 --- /dev/null +++ b/samplebrain/src/aquila/source/generator/WhiteNoiseGenerator.h @@ -0,0 +1,37 @@ +/** + * @file WhiteNoiseGenerator.h + * + * White noise generator. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef WHITENOISEGENERATOR_H +#define WHITENOISEGENERATOR_H + +#include "Generator.h" + +namespace Aquila +{ + /** + * White noise generator. + */ + class AQUILA_EXPORT WhiteNoiseGenerator : public Generator + { + public: + WhiteNoiseGenerator(FrequencyType sampleFrequency); + + virtual void generate(std::size_t samplesCount); + }; +} + +#endif // WHITENOISEGENERATOR_H diff --git a/samplebrain/src/aquila/source/window/BarlettWindow.cpp b/samplebrain/src/aquila/source/window/BarlettWindow.cpp new file mode 100644 index 0000000..020c0c1 --- /dev/null +++ b/samplebrain/src/aquila/source/window/BarlettWindow.cpp @@ -0,0 +1,39 @@ +/** + * @file BarlettWindow.cpp + * + * Barlett (triangular) window. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "BarlettWindow.h" +#include + +namespace Aquila +{ + /** + * Creates Barlett (triangular) window of given size. + * + * @param size window length + */ + BarlettWindow::BarlettWindow(std::size_t size): + SignalSource() + { + m_data.reserve(size); + for (std::size_t n = 0; n < size; ++n) + { + m_data.push_back( + 1.0 - (2.0 * std::fabs(n - (size - 1) / 2.0)) / (double(size - 1)) + ); + } + } +} diff --git a/samplebrain/src/aquila/source/window/BarlettWindow.h b/samplebrain/src/aquila/source/window/BarlettWindow.h new file mode 100644 index 0000000..8cc62e0 --- /dev/null +++ b/samplebrain/src/aquila/source/window/BarlettWindow.h @@ -0,0 +1,37 @@ +/** + * @file BarlettWindow.h + * + * Barlett (triangular) window. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef BARLETTWINDOW_H +#define BARLETTWINDOW_H + +#include "../../global.h" +#include "../SignalSource.h" +#include + +namespace Aquila +{ + /** + * Barlett (triangular) window. + */ + class AQUILA_EXPORT BarlettWindow : public SignalSource + { + public: + BarlettWindow(std::size_t size); + }; +} + +#endif // BARLETTWINDOW_H diff --git a/samplebrain/src/aquila/source/window/BlackmanWindow.cpp b/samplebrain/src/aquila/source/window/BlackmanWindow.cpp new file mode 100644 index 0000000..7b55906 --- /dev/null +++ b/samplebrain/src/aquila/source/window/BlackmanWindow.cpp @@ -0,0 +1,40 @@ +/** + * @file BlackmanWindow.cpp + * + * Blackman window. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "BlackmanWindow.h" +#include + +namespace Aquila +{ + /** + * Creates Blackman window of given size. + * + * @param size window length + */ + BlackmanWindow::BlackmanWindow(std::size_t size): + SignalSource() + { + m_data.reserve(size); + for (std::size_t n = 0; n < size; ++n) + { + m_data.push_back( + 0.42 - 0.5 * std::cos(2.0 * M_PI * n / double(size - 1)) + + 0.08 * std::cos(4.0 * M_PI * n / double(size - 1)) + ); + } + } +} diff --git a/samplebrain/src/aquila/source/window/BlackmanWindow.h b/samplebrain/src/aquila/source/window/BlackmanWindow.h new file mode 100644 index 0000000..8d99b84 --- /dev/null +++ b/samplebrain/src/aquila/source/window/BlackmanWindow.h @@ -0,0 +1,37 @@ +/** + * @file BlackmanWindow.h + * + * Blackman window. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef BLACKMANWINDOW_H +#define BLACKMANWINDOW_H + +#include "../../global.h" +#include "../SignalSource.h" +#include + +namespace Aquila +{ + /** + * Blackman window. + */ + class AQUILA_EXPORT BlackmanWindow : public SignalSource + { + public: + BlackmanWindow(std::size_t size); + }; +} + +#endif // BLACKMANWINDOW_H diff --git a/samplebrain/src/aquila/source/window/FlattopWindow.cpp b/samplebrain/src/aquila/source/window/FlattopWindow.cpp new file mode 100644 index 0000000..c8a84e1 --- /dev/null +++ b/samplebrain/src/aquila/source/window/FlattopWindow.cpp @@ -0,0 +1,42 @@ +/** + * @file FlattopWindow.cpp + * + * Flat-top window. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "FlattopWindow.h" +#include + +namespace Aquila +{ + /** + * Creates flat-top window of given size. + * + * @param size window length + */ + FlattopWindow::FlattopWindow(std::size_t size): + SignalSource() + { + m_data.reserve(size); + for (std::size_t n = 0; n < size; ++n) + { + m_data.push_back( + 1.0 - 1.93 * std::cos(2.0 * M_PI * n / double(size - 1)) + + 1.29 * std::cos(4.0 * M_PI * n / double(size - 1)) - + 0.388 * std::cos(6.0 * M_PI * n / double(size - 1)) + + 0.0322 * std::cos(8.0 * M_PI * n / double(size - 1)) + ); + } + } +} diff --git a/samplebrain/src/aquila/source/window/FlattopWindow.h b/samplebrain/src/aquila/source/window/FlattopWindow.h new file mode 100644 index 0000000..825d4c6 --- /dev/null +++ b/samplebrain/src/aquila/source/window/FlattopWindow.h @@ -0,0 +1,37 @@ +/** + * @file FlattopWindow.h + * + * Flat-top window. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef FLATTOPWINDOW_H +#define FLATTOPWINDOW_H + +#include "../../global.h" +#include "../SignalSource.h" +#include + +namespace Aquila +{ + /** + * Flat-top window. + */ + class AQUILA_EXPORT FlattopWindow : public SignalSource + { + public: + FlattopWindow(std::size_t size); + }; +} + +#endif // FLATTOPWINDOW_H diff --git a/samplebrain/src/aquila/source/window/GaussianWindow.cpp b/samplebrain/src/aquila/source/window/GaussianWindow.cpp new file mode 100644 index 0000000..58999db --- /dev/null +++ b/samplebrain/src/aquila/source/window/GaussianWindow.cpp @@ -0,0 +1,41 @@ +/** + * @file GaussianWindow.cpp + * + * Gaussian (triangular) window. Based on the formula given at: + * http://en.wikipedia.org/wiki/Window_function#Gaussian_window + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Chris Vandevelde + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "GaussianWindow.h" +#include + +namespace Aquila +{ + /** + * Creates Gaussian window of given size, with optional sigma parameter. + * + * @param size window length + * @param sigma standard deviation + */ + GaussianWindow::GaussianWindow(std::size_t size, double sigma/* = 0.5*/): + SignalSource() + { + m_data.reserve(size); + for (std::size_t n = 0; n < size; ++n) + { + m_data.push_back( + std::exp((-0.5) * std::pow((n - (size - 1.0) / 2.0)/(sigma * (size - 1.0) / 2.0), 2.0)) + ); + } + } +} diff --git a/samplebrain/src/aquila/source/window/GaussianWindow.h b/samplebrain/src/aquila/source/window/GaussianWindow.h new file mode 100644 index 0000000..748a193 --- /dev/null +++ b/samplebrain/src/aquila/source/window/GaussianWindow.h @@ -0,0 +1,38 @@ +/** + * @file GaussianWindow.h + * + * Gaussian (triangular) window. Based on the formula given at: + * http://en.wikipedia.org/wiki/Window_function#Gaussian_window + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Chris Vandevelde + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef GAUSSIANWINDOW_H +#define GAUSSIANWINDOW_H + +#include "../../global.h" +#include "../SignalSource.h" +#include + +namespace Aquila +{ + /** + * Creates Gaussian window of given size, with optional sigma parameter. + */ + class AQUILA_EXPORT GaussianWindow : public SignalSource + { + public: + GaussianWindow(std::size_t size, double sigma = 0.5); + }; +} + +#endif // GAUSSIANWINDOW_H diff --git a/samplebrain/src/aquila/source/window/HammingWindow.cpp b/samplebrain/src/aquila/source/window/HammingWindow.cpp new file mode 100644 index 0000000..a666573 --- /dev/null +++ b/samplebrain/src/aquila/source/window/HammingWindow.cpp @@ -0,0 +1,39 @@ +/** + * @file HammingWindow.cpp + * + * Hamming window. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "HammingWindow.h" +#include + +namespace Aquila +{ + /** + * Creates Hamming window of given size. + * + * @param size window length + */ + HammingWindow::HammingWindow(std::size_t size): + SignalSource() + { + m_data.reserve(size); + for (std::size_t n = 0; n < size; ++n) + { + m_data.push_back( + 0.53836 - 0.46164 * std::cos(2.0 * M_PI * n / double(size - 1)) + ); + } + } +} diff --git a/samplebrain/src/aquila/source/window/HammingWindow.h b/samplebrain/src/aquila/source/window/HammingWindow.h new file mode 100644 index 0000000..25a4201 --- /dev/null +++ b/samplebrain/src/aquila/source/window/HammingWindow.h @@ -0,0 +1,37 @@ +/** + * @file HammingWindow.h + * + * Hamming window. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef HAMMINGWINDOW_H +#define HAMMINGWINDOW_H + +#include "../../global.h" +#include "../SignalSource.h" +#include + +namespace Aquila +{ + /** + * Hamming window. + */ + class AQUILA_EXPORT HammingWindow : public SignalSource + { + public: + HammingWindow(std::size_t size); + }; +} + +#endif // HAMMINGWINDOW_H diff --git a/samplebrain/src/aquila/source/window/HannWindow.cpp b/samplebrain/src/aquila/source/window/HannWindow.cpp new file mode 100644 index 0000000..d77ec68 --- /dev/null +++ b/samplebrain/src/aquila/source/window/HannWindow.cpp @@ -0,0 +1,39 @@ +/** + * @file HannWindow.cpp + * + * Hann window. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "HannWindow.h" +#include + +namespace Aquila +{ + /** + * Creates Hann window of given size. + * + * @param size window length + */ + HannWindow::HannWindow(std::size_t size): + SignalSource() + { + m_data.reserve(size); + for (std::size_t n = 0; n < size; ++n) + { + m_data.push_back( + 0.5 * (1.0 - std::cos(2.0 * M_PI * n / double(size - 1))) + ); + } + } +} diff --git a/samplebrain/src/aquila/source/window/HannWindow.h b/samplebrain/src/aquila/source/window/HannWindow.h new file mode 100644 index 0000000..aa8bf8a --- /dev/null +++ b/samplebrain/src/aquila/source/window/HannWindow.h @@ -0,0 +1,37 @@ +/** + * @file HannWindow.h + * + * Hann window. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef HANNWINDOW_H +#define HANNWINDOW_H + +#include "../../global.h" +#include "../SignalSource.h" +#include + +namespace Aquila +{ + /** + * Hann window. + */ + class AQUILA_EXPORT HannWindow : public SignalSource + { + public: + HannWindow(std::size_t size); + }; +} + +#endif // HANNWINDOW_H diff --git a/samplebrain/src/aquila/source/window/RectangularWindow.h b/samplebrain/src/aquila/source/window/RectangularWindow.h new file mode 100644 index 0000000..86e6213 --- /dev/null +++ b/samplebrain/src/aquila/source/window/RectangularWindow.h @@ -0,0 +1,46 @@ +/** + * @file RectangularWindow.h + * + * A rectangular window (equivalent to multiplication by one). + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef RECTANGULARWINDOW_H +#define RECTANGULARWINDOW_H + +#include "../../global.h" +#include "../SignalSource.h" +#include + +namespace Aquila +{ + /** + * Rectangular window. + */ + class AQUILA_EXPORT RectangularWindow : public SignalSource + { + public: + /** + * Creates rectangular window of given size. + * + * @param size window length + */ + RectangularWindow(std::size_t size): + SignalSource() + { + m_data.assign(size, 1.0); + } + }; +} + +#endif // RECTANGULARWINDOW_H diff --git a/samplebrain/src/aquila/synth.h b/samplebrain/src/aquila/synth.h new file mode 100644 index 0000000..e5d42e8 --- /dev/null +++ b/samplebrain/src/aquila/synth.h @@ -0,0 +1,26 @@ +/** + * @file synth.h + * + * Convenience header that includes all synthesizer headers. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + + +#ifndef AQUILA_SYNTH_H +#define AQUILA_SYNTH_H + +#include "synth/Synthesizer.h" +#include "synth/KarplusStrongSynthesizer.h" +#include "synth/SineSynthesizer.h" + +#endif // AQUILA_SYNTH_H diff --git a/samplebrain/src/aquila/synth/KarplusStrongSynthesizer.cpp b/samplebrain/src/aquila/synth/KarplusStrongSynthesizer.cpp new file mode 100644 index 0000000..d0f1a46 --- /dev/null +++ b/samplebrain/src/aquila/synth/KarplusStrongSynthesizer.cpp @@ -0,0 +1,63 @@ +/** + * @file KarplusStrongSynthesizer.cpp + * + * Plucked string synthesis using Karplus-Strong algorithm. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "KarplusStrongSynthesizer.h" +#include +#include +#include +#include + +namespace Aquila +{ + /** + * Plays a single note by "plucking" a string. + * + * A short noise burst is fed through a feedback loop including delay and + * a first-order lowpass filter (in this case a simple moving average). + * Resulting waveform is played using SFML - the sound is similar to a + * plucked guitar string. + * + * @param frequency base frequency of the guitar note + * @param duration tone duration in milliseconds + */ + void KarplusStrongSynthesizer::playFrequency(FrequencyType frequency, + unsigned int duration) + { + std::size_t delay = static_cast(m_sampleFrequency / frequency); + std::size_t totalSamples = static_cast(m_sampleFrequency * duration / 1000.0); + m_generator.setAmplitude(8192).generate(delay); + + // copy initial noise burst at the beginning of output array + sf::Int16* arr = new sf::Int16[totalSamples]; + std::copy(std::begin(m_generator), std::end(m_generator), arr); + // first sample that goes into feedback loop; + // cannot be averaged with previous + arr[delay] = m_alpha * arr[0]; + for (std::size_t i = delay + 1; i < totalSamples; ++i) + { + // average two consecutive delayed samples and dampen by alpha + arr[i] = m_alpha * (0.5 * (arr[i - delay] + arr[i - delay - 1])); + } + + m_buffer.loadFromSamples(arr, totalSamples, 1, m_sampleFrequency); + sf::Sound sound(m_buffer); + sound.play(); + sf::sleep(sf::milliseconds(duration)); + + delete [] arr; + } +} diff --git a/samplebrain/src/aquila/synth/KarplusStrongSynthesizer.h b/samplebrain/src/aquila/synth/KarplusStrongSynthesizer.h new file mode 100644 index 0000000..4487aa9 --- /dev/null +++ b/samplebrain/src/aquila/synth/KarplusStrongSynthesizer.h @@ -0,0 +1,68 @@ +/** + * @file KarplusStrongSynthesizer.h + * + * Plucked string synthesis using Karplus-Strong algorithm. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef KARPLUSSTRONGSYNTHESIZER_H +#define KARPLUSSTRONGSYNTHESIZER_H + +#include "../global.h" +#include "Synthesizer.h" +#include "../source/generator/WhiteNoiseGenerator.h" + +namespace Aquila +{ + /** + * Very simple guitar synthesizer using Karplus-Strong algorithm. + */ + class AQUILA_EXPORT KarplusStrongSynthesizer : public Synthesizer + { + public: + /** + * Creates the synthesizer object. + * + * @param sampleFrequency sample frequency of the audio signal + */ + KarplusStrongSynthesizer(FrequencyType sampleFrequency): + Synthesizer(sampleFrequency), m_generator(sampleFrequency), + m_alpha(0.99) + { + } + + /** + * Sets feedback loop parameter. + */ + void setAlpha(double alpha) + { + m_alpha = alpha; + } + + protected: + virtual void playFrequency(FrequencyType frequency, unsigned int duration); + + private: + /** + * Generator used to provide initial noise burst. + */ + WhiteNoiseGenerator m_generator; + + /** + * Feedback loop parameter. + */ + double m_alpha; + }; +} + +#endif // KARPLUSSTRONGSYNTHESIZER_H diff --git a/samplebrain/src/aquila/synth/SineSynthesizer.cpp b/samplebrain/src/aquila/synth/SineSynthesizer.cpp new file mode 100644 index 0000000..3d6900b --- /dev/null +++ b/samplebrain/src/aquila/synth/SineSynthesizer.cpp @@ -0,0 +1,42 @@ +/** + * @file SineSynthesizer.cpp + * + * Simple sine tone synthesis. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "SineSynthesizer.h" +#include +#include +#include + +namespace Aquila +{ + /** + * Plays a tone at given frequency. + * + * @param frequency frequency of the generated sound + * @param duration beep duration in milliseconds + */ + void SineSynthesizer::playFrequency(FrequencyType frequency, + unsigned int duration) + { + std::size_t numSamples = static_cast(m_sampleFrequency * duration / 1000); + m_generator.setFrequency(frequency).generate(numSamples); + m_buffer.loadFromSignalSource(m_generator); + sf::Sound sound(m_buffer); + sound.play(); + // the additional 50 ms is an intentional pause between tones + sf::sleep(sf::milliseconds(duration + 50)); + } +} diff --git a/samplebrain/src/aquila/synth/SineSynthesizer.h b/samplebrain/src/aquila/synth/SineSynthesizer.h new file mode 100644 index 0000000..472ff16 --- /dev/null +++ b/samplebrain/src/aquila/synth/SineSynthesizer.h @@ -0,0 +1,55 @@ +/** + * @file SineSynthesizer.h + * + * Simple sine tone synthesis. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef SINESYNTHESIZER_H +#define SINESYNTHESIZER_H + +#include "../global.h" +#include "Synthesizer.h" +#include "../source/generator/SineGenerator.h" + +namespace Aquila +{ + /** + * Sine wave synthesizer. + */ + class AQUILA_EXPORT SineSynthesizer : public Synthesizer + { + public: + /** + * Creates the synthesizer object. + * + * @param sampleFrequency sample frequency of the audio signal + */ + SineSynthesizer(FrequencyType sampleFrequency): + Synthesizer(sampleFrequency), m_generator(sampleFrequency) + { + m_generator.setAmplitude(8192); + } + + protected: + void playFrequency(FrequencyType note, unsigned int duration); + + private: + /** + * Underlying sine generator. + */ + SineGenerator m_generator; + }; +} + +#endif // SINESYNTHESIZER_H diff --git a/samplebrain/src/aquila/synth/Synthesizer.cpp b/samplebrain/src/aquila/synth/Synthesizer.cpp new file mode 100644 index 0000000..99026fd --- /dev/null +++ b/samplebrain/src/aquila/synth/Synthesizer.cpp @@ -0,0 +1,108 @@ +/** + * @file Synthesizer.cpp + * + * Base class for SFML-based audio synthesizers. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "Synthesizer.h" +#include + +namespace Aquila +{ + /** + * Creates the mapping between note names and frequencies. + * + * @return initialized note map + */ + NoteMapType initNoteMap() + { + NoteMapType notes; + notes["C2"] = 65.406; + notes["C2S"] = 69.296; + notes["D2"] = 73.416; + notes["D2S"] = 77.782; + notes["E2"] = 82.407; + notes["F2"] = 87.307; + notes["F2S"] = 92.499; + notes["G2"] = 97.999; + notes["G2S"] = 103.83; + notes["A2"] = 110.0; + notes["A2S"] = 116.54; + notes["B2"] = 123.47; + notes["C3"] = 130.81; + notes["C3S"] = 138.59; + notes["D3"] = 146.83 ; + notes["D3S"] = 155.56; + notes["E3"] = 164.81; + notes["F3"] = 174.61; + notes["F3S"] = 185.0; + notes["G3"] = 196.0; + notes["G3S"] = 207.65; + notes["A3"] = 220.00; + notes["A3S"] = 233.08; + notes["B3"] = 246.94; + notes["C4"] = 261.626; + notes["C4S"] = 277.18; + notes["D4"] = 293.665; + notes["D4S"] = 311.13; + notes["E4"] = 329.628; + notes["F4"] = 349.228; + notes["F4S"] = 369.99; + notes["G4"] = 391.995; + notes["G4S"] = 415.305; + notes["A4"] = 440.0; + notes["A4S"] = 466.164; + notes["B4"] = 493.883; + notes["C5"] = 523.251; + notes["C5S"] = 554.365; + notes["D5"] = 587.33; + notes["D5S"] = 622.254; + notes["E5"] = 659.255; + notes["F5"] = 698.456; + notes["F5S"] = 739.989; + notes["G5"] = 783.991; + notes["G5S"] = 830.609; + notes["A5"] = 880.0; + notes["A5S"] = 932.33; + notes["B5"] = 987.77; + notes["C6"] = 1046.5; + return notes; + } + + NoteMapType Synthesizer::m_notes = initNoteMap(); + + /** + * Plays a single note. + * + * This method only does the lookup from note name to frequency and + * delegates the actual playing to pure virtual method playFrequency. + * + * Unrecognized note names are silent for the given duration. + * + * @param note note name (@see m_notes) + * @param duration duration in milliseconds + */ + void Synthesizer::playNote(std::string note, unsigned int duration) + { + if (m_notes.count(note) == 0) + { + sf::sleep(sf::milliseconds(duration)); + } + else + { + FrequencyType frequency = m_notes[note]; + playFrequency(frequency, duration); + } + } +} diff --git a/samplebrain/src/aquila/synth/Synthesizer.h b/samplebrain/src/aquila/synth/Synthesizer.h new file mode 100644 index 0000000..fd25f03 --- /dev/null +++ b/samplebrain/src/aquila/synth/Synthesizer.h @@ -0,0 +1,88 @@ +/** + * @file Synthesizer.h + * + * Base class for SFML-based audio synthesizers. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef SYNTHESIZER_H +#define SYNTHESIZER_H + +#include "../global.h" +#include "../wrappers/SoundBufferAdapter.h" +#include +#include +#include +#include + +namespace Aquila +{ + /** + * Type of the mapping from note names to frequencies. + */ + typedef std::map NoteMapType; + + NoteMapType initNoteMap(); + + /** + * An abstract class from which sound synthesizers should be derived. + */ + class AQUILA_EXPORT Synthesizer + { + public: + /** + * Creates the synthesizer object. + * + * @param sampleFrequency sample frequency of the audio signal + */ + Synthesizer(FrequencyType sampleFrequency): + m_sampleFrequency(sampleFrequency), m_buffer() + { + } + + /** + * No-op virtual destructor. + */ + virtual ~Synthesizer() {} + + void playNote(std::string note, unsigned int duration = 500); + + protected: + /** + * Plays a tone at given frequency. + * + * Must be overriden in child classes. + * + * @param frequency base frequency of the generated sound + * @param duration tone duration in milliseconds + */ + virtual void playFrequency(FrequencyType frequency, unsigned int duration) = 0; + + /** + * Sample frequency of the generated signal. + */ + const FrequencyType m_sampleFrequency; + + /** + * Audio buffer for playback. + */ + SoundBufferAdapter m_buffer; + + /** + * A mapping from note names to frequencies. + */ + static NoteMapType m_notes; + }; +} + +#endif // SYNTHESIZER_H diff --git a/samplebrain/src/aquila/tools.h b/samplebrain/src/aquila/tools.h new file mode 100644 index 0000000..4f84561 --- /dev/null +++ b/samplebrain/src/aquila/tools.h @@ -0,0 +1,24 @@ +/** + * @file tools.h + * + * Convenience header that includes all utility headers. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + + +#ifndef AQUILA_TOOLS_H +#define AQUILA_TOOLS_H + +#include "tools/TextPlot.h" + +#endif // AQUILA_TOOLS_H diff --git a/samplebrain/src/aquila/tools/TextPlot.cpp b/samplebrain/src/aquila/tools/TextPlot.cpp new file mode 100644 index 0000000..25fec2c --- /dev/null +++ b/samplebrain/src/aquila/tools/TextPlot.cpp @@ -0,0 +1,70 @@ +/** + * @file TextPlot.cpp + * + * A simple console-based data plotter. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "TextPlot.h" + +namespace Aquila +{ + /** + * Creates the plot object. + * + * @param title plot title (optional, default is "Data plot") + * @param out where to output the plot data (default is std::cout) + */ + TextPlot::TextPlot(const std::string &title, std::ostream &out): + m_title(title), m_out(out), m_width(64), m_height(16) + { + } + + /** + * "Draws" the plot to the output stream. + * + * @param plot internal plot data + */ + void TextPlot::drawPlotMatrix(const PlotMatrixType &plotData) + { + const std::size_t length = plotData.size(); + + m_out << "\n" << m_title << "\n"; + // output the plot data, flushing only at the end + for (unsigned int y = 0; y < m_height; ++y) + { + for (unsigned int x = 0; x < length; ++x) + { + m_out << plotData[x][y]; + } + m_out << "\n"; + } + m_out.flush(); + } + + /** + * A shorthand for plotting only the first half of magnitude spectrum. + * + * @param spectrum a vector of complex numbers + */ + void TextPlot::plotSpectrum(Aquila::SpectrumType spectrum) + { + std::size_t halfLength = spectrum.size() / 2; + std::vector absSpectrum(halfLength); + for (std::size_t i = 0; i < halfLength; ++i) + { + absSpectrum[i] = std::abs(spectrum[i]); + } + plot(absSpectrum); + } +} diff --git a/samplebrain/src/aquila/tools/TextPlot.h b/samplebrain/src/aquila/tools/TextPlot.h new file mode 100644 index 0000000..80110d8 --- /dev/null +++ b/samplebrain/src/aquila/tools/TextPlot.h @@ -0,0 +1,203 @@ +/** + * @file TextPlot.h + * + * A simple console-based data plotter. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef TEXTPLOT_H +#define TEXTPLOT_H + +#include "../global.h" +#include "../source/SignalSource.h" +#include +#include +#include +#include +#include +#include + +namespace Aquila +{ + /** + * A utility class to "draw" data plots in the console applications. + * + * TextPlot is used for quick-and-dirty data plotting in console programs. + * It also serves as a way to take a peek into signal data. The class is + * very simple and - intentionally - lacks a lot of configuration options. + * A more serious tool will be neccessary to create pretty images. + * For example, a frontend to gnuplot is taken under consideration. + */ + class AQUILA_EXPORT TextPlot + { + public: + TextPlot(const std::string& title = "Data plot", + std::ostream& out = std::cout); + + /** + * Returns the title of the plot. + * + * @return plot title + */ + std::string getTitle() const + { + return m_title; + } + + /** + * Sets the title of the plot. + * + * @param title plot title + */ + void setTitle(const std::string& title) + { + m_title = title; + } + + /** + * Returns plot width. + * + * @return plot width + */ + std::size_t getWidth() const + { + return m_width; + } + + /** + * Returns plot height. + * + * @return plot height + */ + std::size_t getHeight() const + { + return m_height; + } + + /** + * Sets plot dimensions. + * + * @param width plot width, currently ignored + * @param height plot height + */ + void setSize(std::size_t width, std::size_t height) + { + m_width = width; + m_height = height; + } + + /** + * Plots all data from a given source. + * + * @param source data source + */ + void plot(const SignalSource& source) + { + PlotMatrixType plotData(source.length()); + doPlot(plotData, source.begin(), source.end()); + } + + /** + * An overload for plotting vectors. + * + * @param data a numeric vector + */ + template + void plot(const std::vector& data) + { + PlotMatrixType plotData(data.size()); + doPlot(plotData, data.begin(), data.end()); + } + + /** + * An overload for plotting regular C-style arrays. + * + * In theory, we could use the template argument trick to omit the + * `length` argument, however, that is possible only for static arrays. + * + * @param data the array + * @param length array size + */ + template + void plot(Numeric data[], std::size_t length) + { + PlotMatrixType plotData(length); + doPlot(plotData, data, data + length); + } + + void plotSpectrum(SpectrumType spectrum); + + private: + /** + * The internal representation of the plot. + */ + typedef std::vector> PlotMatrixType; + + /** + * Prepares an internal "pixel" matrix and calls drawPlotMatrix(). + * + * The template parameter is an iterator type, therefore the plotter + * can be more generic. + * + * @param plot reference to the plot matrix + * @param begin an iterator pointing to the beginning of plotted data + * @param end an iterator pointing "one past end" of plotted data + */ + template + void doPlot(PlotMatrixType& plotData, Iterator begin, Iterator end) + { + const double max = *std::max_element(begin, end); + const double min = *std::min_element(begin, end); + const double range = max - min; + + for (std::size_t xPos = 0; xPos < plotData.size(); ++xPos) + { + plotData[xPos].resize(m_height, ' '); + double normalizedValue = (*begin++ - min) / range; + std::size_t yPos = m_height - static_cast( + std::ceil(m_height * normalizedValue)); + + // bound the value so it stays within vector dimension + if (yPos >= m_height) + yPos = m_height - 1; + plotData[xPos][yPos] = '*'; + } + + drawPlotMatrix(plotData); + } + + void drawPlotMatrix(const PlotMatrixType& plotData); + + /** + * Plot title. + */ + std::string m_title; + + /** + * Output stream. + */ + std::ostream& m_out; + + /** + * Plot width - currently ignored. + */ + std::size_t m_width; + + /** + * Plot height. + */ + std::size_t m_height; + }; +} + +#endif // TEXTPLOT_H diff --git a/samplebrain/src/aquila/transform.h b/samplebrain/src/aquila/transform.h new file mode 100644 index 0000000..d4b13d6 --- /dev/null +++ b/samplebrain/src/aquila/transform.h @@ -0,0 +1,31 @@ +/** + * @file transform.h + * + * Convenience header that includes all transformation headers. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + + +#ifndef AQUILA_TRANSFORM_H +#define AQUILA_TRANSFORM_H + +#include "transform/Fft.h" +#include "transform/Dft.h" +#include "transform/AquilaFft.h" +#include "transform/OouraFft.h" +#include "transform/FftFactory.h" +#include "transform/Dct.h" +#include "transform/Mfcc.h" +#include "transform/Spectrogram.h" + +#endif // AQUILA_TRANSFORM_H diff --git a/samplebrain/src/aquila/transform/AquilaFft.cpp b/samplebrain/src/aquila/transform/AquilaFft.cpp new file mode 100644 index 0000000..c4252c0 --- /dev/null +++ b/samplebrain/src/aquila/transform/AquilaFft.cpp @@ -0,0 +1,214 @@ +/** + * @file AquilaFft.cpp + * + * A custom implementation of FFT radix-2 algorithm. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "AquilaFft.h" +#include +#include +#include + +namespace Aquila +{ + /** + * Complex unit. + */ + const ComplexType AquilaFft::j(0, 1); + + /** + * Applies the transformation to the signal. + * + * @param x input signal + * @return calculated spectrum + */ + SpectrumType AquilaFft::fft(const SampleType x[]) + { + SpectrumType spectrum(N); + + // bit-reversing the samples - a requirement of radix-2 + // instead of reversing in place, put the samples to result array + unsigned int a = 1, b = 0, c = 0; + std::copy(x, x + N, std::begin(spectrum)); + for (b = 1; b < N; ++b) + { + if (b < a) + { + spectrum[a - 1] = x[b - 1]; + spectrum[b - 1] = x[a - 1]; + } + c = N / 2; + while (c < a) + { + a -= c; + c /= 2; + } + a += c; + } + + // FFT calculation using "butterflies" + // code ported from Matlab, based on book by Tomasz P. Zieliński + + // FFT stages count + unsigned int numStages = static_cast( + std::log(static_cast(N)) / LN_2); + + // L = 2^k - DFT block length and offset + // M = 2^(k-1) - butterflies per block, butterfly width + // p - butterfly index + // q - block index + // r - index of sample in butterfly + // Wi - starting value of Fourier base coefficient + unsigned int L = 0, M = 0, p = 0, q = 0, r = 0; + ComplexType Wi(0, 0), Temp(0, 0); + + ComplexType** Wi_cache = getCachedFftWi(numStages); + + // iterate over the stages + for (unsigned int k = 1; k <= numStages; ++k) + { + L = 1 << k; + M = 1 << (k - 1); + Wi = Wi_cache[k][0]; + + // iterate over butterflies + for (p = 1; p <= M; ++p) + { + // iterate over blocks + for (q = p; q <= N; q += L) + { + r = q + M; + Temp = spectrum[r - 1] * Wi; + spectrum[r - 1] = spectrum[q - 1] - Temp; + spectrum[q - 1] = spectrum[q - 1] + Temp; + } + Wi = Wi_cache[k][p]; + } + } + + return spectrum; + } + + /** + * Applies the inverse transform to the spectrum. + * + * @param spectrum input spectrum + * @param x output signal + */ + void AquilaFft::ifft(SpectrumType spectrum, double x[]) + { + SpectrumType spectrumCopy(spectrum); + unsigned int a = 1, b = 0, c = 0; + for (b = 1; b < N; ++b) + { + if (b < a) + { + spectrumCopy[a - 1] = spectrum[b - 1]; + spectrumCopy[b - 1] = spectrum[a - 1]; + } + c = N / 2; + while (c < a) + { + a -= c; + c /= 2; + } + a += c; + } + + unsigned int numStages = static_cast( + std::log(static_cast(N)) / LN_2); + unsigned int L = 0, M = 0, p = 0, q = 0, r = 0; + ComplexType Wi(0, 0), Temp(0, 0); + ComplexType** Wi_cache = getCachedFftWi(numStages); + for (unsigned int k = 1; k <= numStages; ++k) + { + L = 1 << k; + M = 1 << (k - 1); + Wi = -Wi_cache[k][0]; + for (p = 1; p <= M; ++p) + { + for (q = p; q <= N; q += L) + { + r = q + M; + Temp = spectrumCopy[r - 1] * Wi; + spectrumCopy[r - 1] = spectrumCopy[q - 1] - Temp; + spectrumCopy[q - 1] = spectrumCopy[q - 1] + Temp; + } + Wi = -Wi_cache[k][p]; + } + } + + for (unsigned int k = 0; k < N; ++k) + { + x[k] = std::abs(spectrumCopy[k]) / static_cast(N); + } + } + + /** + * Returns a table of Wi (twiddle factors) stored in cache. + * + * @param numStages the FFT stages count + * @return pointer to an array of pointers to arrays of complex numbers + */ + ComplexType** AquilaFft::getCachedFftWi(unsigned int numStages) + { + fftWiCacheKeyType key = numStages; + // cache hit, return immediately + if (fftWiCache.find(key) != fftWiCache.end()) + { + return fftWiCache[key]; + } + + // nothing in cache, calculate twiddle factors + ComplexType** Wi = new ComplexType*[numStages + 1]; + for (unsigned int k = 1; k <= numStages; ++k) + { + // L = 2^k - DFT block length and offset + // M = 2^(k-1) - butterflies per block, butterfly width + // W - Fourier base multiplying factor + unsigned int L = 1 << k; + unsigned int M = 1 << (k - 1); + ComplexType W = exp((-j) * 2.0 * M_PI / static_cast(L)); + Wi[k] = new ComplexType[M + 1]; + Wi[k][0] = ComplexType(1.0); + for (unsigned int p = 1; p <= M; ++p) + { + Wi[k][p] = Wi[k][p - 1] * W; + } + } + + // store in cache and return + fftWiCache[key] = Wi; + + return Wi; + } + + /** + * Clears the twiddle factor cache. + */ + void AquilaFft::clearFftWiCache() + { + for (auto it = fftWiCache.begin(); it != fftWiCache.end(); it++) + { + ComplexType** c = it->second; + unsigned int numStages = it->first; + for (unsigned int i = 1; i <= numStages; ++i) + { + delete [] c[i]; + } + + delete [] c; + } + } +} diff --git a/samplebrain/src/aquila/transform/AquilaFft.h b/samplebrain/src/aquila/transform/AquilaFft.h new file mode 100644 index 0000000..ff89504 --- /dev/null +++ b/samplebrain/src/aquila/transform/AquilaFft.h @@ -0,0 +1,88 @@ +/** + * @file AquilaFft.h + * + * A custom implementation of FFT radix-2 algorithm. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef AQUILAFFT_H +#define AQUILAFFT_H + +#include "Fft.h" +#include +#include + +namespace Aquila +{ + /** + * ln(2) - needed for calculating number of stages in FFT. + */ + const double LN_2 = 0.69314718055994530941723212145818; + + /** + * A custom implementation of FFT radix-2 algorithm. + */ + class AQUILA_EXPORT AquilaFft : public Fft + { + public: + /** + * Initializes the transform for a given input length. + * + * @param length input signal size (usually a power of 2) + */ + AquilaFft(std::size_t length): + Fft(length), fftWiCache() + { + } + + /** + * Destroys the transform object. + * + * @todo clear the cache somewhere else, not here! + */ + ~AquilaFft() + { + clearFftWiCache(); + } + + virtual SpectrumType fft(const SampleType x[]); + virtual void ifft(SpectrumType spectrum, double x[]); + + private: + /** + * Complex unit (0.0 + 1.0j). + */ + static const ComplexType j; + + /** + * Type of the twiddle factor cache key. + */ + typedef unsigned int fftWiCacheKeyType; + + /** + * Type of the twiddle factor cache. + */ + typedef std::map fftWiCacheType; + + /** + * Twiddle factor cache - implemented as a map. + */ + fftWiCacheType fftWiCache; + + ComplexType** getCachedFftWi(unsigned int numStages); + + void clearFftWiCache(); + }; +} + +#endif // AQUILAFFT_H diff --git a/samplebrain/src/aquila/transform/Dct.cpp b/samplebrain/src/aquila/transform/Dct.cpp new file mode 100644 index 0000000..ac3098f --- /dev/null +++ b/samplebrain/src/aquila/transform/Dct.cpp @@ -0,0 +1,114 @@ +/** + * @file Dct.cpp + * + * Discrete Cosine Transform (DCT) calculation. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "Dct.h" +#include +#include +#include + +namespace Aquila +{ + /** + * Calculates the Discrete Cosine Transform, type II. + * + * See http://en.wikipedia.org/wiki/Discrete_cosine_transform for + * explanation what DCT-II is. + * + * Uses cosine value caching in order to speed up computations. + * + * @param data input data vector + * @param outputLength how many coefficients to return + * @return vector of DCT coefficients + */ + std::vector Dct::dct(const std::vector& data, std::size_t outputLength) + { + // zero-initialize output vector + std::vector output(outputLength, 0.0); + std::size_t inputLength = data.size(); + + // DCT scaling factor + double c0 = std::sqrt(1.0 / inputLength); + double cn = std::sqrt(2.0 / inputLength); + // cached cosine values + double** cosines = getCachedCosines(inputLength, outputLength); + + for (std::size_t n = 0; n < outputLength; ++n) + { + for (std::size_t k = 0; k < inputLength; ++k) + { + output[n] += data[k] * cosines[n][k]; + } + output[n] *= (0 == n) ? c0 : cn; + } + + return output; + } + + /** + * Returns a table of DCT cosine values stored in memory cache. + * + * The two params unambigiously identify which cache to use. + * + * @param inputLength length of the input vector + * @param outputLength length of the output vector + * @return pointer to array of pointers to arrays of doubles + */ + double** Dct::getCachedCosines(std::size_t inputLength, std::size_t outputLength) + { + auto key = std::make_pair(inputLength, outputLength); + + // if we have that key cached, return immediately + if (cosineCache.find(key) != cosineCache.end()) + { + return cosineCache[key]; + } + + // nothing in cache for that pair, calculate cosines + double** cosines = new double*[outputLength]; + for (std::size_t n = 0; n < outputLength; ++n) + { + cosines[n] = new double[inputLength]; + for (std::size_t k = 0; k < inputLength; ++k) + { + // from the definition of DCT-II + cosines[n][k] = std::cos((M_PI * (2 * k + 1) * n) / + (2.0 * inputLength)); + } + } + cosineCache[key] = cosines; + + return cosines; + } + + /** + * Deletes all the memory used by cache. + */ + void Dct::clearCosineCache() + { + for (auto it = std::begin(cosineCache); it != std::end(cosineCache); it++) + { + auto key = it->first; + double** cosines = it->second; + std::size_t outputLength = key.second; + for (std::size_t i = 0; i < outputLength; ++i) + { + delete [] cosines[i]; + } + delete [] cosines; + } + } +} diff --git a/samplebrain/src/aquila/transform/Dct.dep b/samplebrain/src/aquila/transform/Dct.dep new file mode 100644 index 0000000..42a89f3 --- /dev/null +++ b/samplebrain/src/aquila/transform/Dct.dep @@ -0,0 +1,2 @@ +Dct.o: src/aquila/transform/Dct.cpp src/aquila/transform/Dct.h \ + src/aquila/transform/../global.h diff --git a/samplebrain/src/aquila/transform/Dct.h b/samplebrain/src/aquila/transform/Dct.h new file mode 100644 index 0000000..16d4203 --- /dev/null +++ b/samplebrain/src/aquila/transform/Dct.h @@ -0,0 +1,75 @@ +/** + * @file Dct.h + * + * Discrete Cosine Transform (DCT) calculation. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef DCT_H +#define DCT_H + +#include "../global.h" +#include +#include +#include +#include + +namespace Aquila +{ + /** + * An implementation of the Discrete Cosine Transform. + */ + class AQUILA_EXPORT Dct + { + public: + /** + * Initializes the transform. + */ + Dct(): + cosineCache() + { + } + + /** + * Destroys the transform object. + */ + ~Dct() + { + clearCosineCache(); + } + + std::vector dct(const std::vector& data, std::size_t outputLength); + + private: + /** + * Key type for the cache, using input and output length. + */ + typedef std::pair cosineCacheKeyType; + + /** + * Cache type. + */ + typedef std::map cosineCacheType; + + /** + * Cache object, implemented as a map. + */ + cosineCacheType cosineCache; + + double** getCachedCosines(std::size_t inputLength, std::size_t outputLength); + + void clearCosineCache(); + }; +} + +#endif // DCT_H diff --git a/samplebrain/src/aquila/transform/Dct.o b/samplebrain/src/aquila/transform/Dct.o new file mode 100644 index 0000000..39e8642 Binary files /dev/null and b/samplebrain/src/aquila/transform/Dct.o differ diff --git a/samplebrain/src/aquila/transform/Dft.cpp b/samplebrain/src/aquila/transform/Dft.cpp new file mode 100644 index 0000000..322edd3 --- /dev/null +++ b/samplebrain/src/aquila/transform/Dft.cpp @@ -0,0 +1,77 @@ +/** + * @file Dft.cpp + * + * A reference implementation of the Discrete Fourier Transform. + * + * Note - this is a direct application of the DFT equations and as such it + * surely isn't optimal. The implementation here serves only as a reference + * model to compare with. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "Dft.h" +#include +#include +#include + +namespace Aquila +{ + /** + * Complex unit. + */ + const ComplexType Dft::j(0, 1); + + /** + * Applies the transformation to the signal. + * + * @param x input signal + * @return calculated spectrum + */ + SpectrumType Dft::fft(const SampleType x[]) + { + SpectrumType spectrum(N); + ComplexType WN = std::exp((-j) * 2.0 * M_PI / static_cast(N)); + + for (unsigned int k = 0; k < N; ++k) + { + ComplexType sum(0, 0); + for (unsigned int n = 0; n < N; ++n) + { + sum += x[n] * std::pow(WN, n * k); + } + spectrum[k] = sum; + } + + return spectrum; + } + + /** + * Applies the inverse transform to the spectrum. + * + * @param spectrum input spectrum + * @param x output signal + */ + void Dft::ifft(SpectrumType spectrum, double x[]) + { + ComplexType WN = std::exp((-j) * 2.0 * M_PI / static_cast(N)); + for (unsigned int k = 0; k < N; ++k) + { + ComplexType sum(0, 0); + for (unsigned int n = 0; n < N; ++n) + { + sum += spectrum[n] * std::pow(WN, -static_cast(n * k)); + } + x[k] = std::abs(sum) / static_cast(N); + } + } +} diff --git a/samplebrain/src/aquila/transform/Dft.dep b/samplebrain/src/aquila/transform/Dft.dep new file mode 100644 index 0000000..bc07d34 --- /dev/null +++ b/samplebrain/src/aquila/transform/Dft.dep @@ -0,0 +1,2 @@ +Dft.o: src/aquila/transform/Dft.cpp src/aquila/transform/Dft.h \ + src/aquila/transform/Fft.h src/aquila/transform/../global.h diff --git a/samplebrain/src/aquila/transform/Dft.h b/samplebrain/src/aquila/transform/Dft.h new file mode 100644 index 0000000..f448180 --- /dev/null +++ b/samplebrain/src/aquila/transform/Dft.h @@ -0,0 +1,63 @@ +/** + * @file Dft.h + * + * A reference implementation of the Discrete Fourier Transform. + * + * Note - this is a direct application of the DFT equation and as such it + * surely isn't optimal. The implementation here serves only as a reference + * model to compare with. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef DFT_H +#define DFT_H + +#include "Fft.h" + +namespace Aquila +{ + /** + * A straightforward implementation of the Discrete Fourier Transform. + */ + class AQUILA_EXPORT Dft : public Fft + { + public: + /** + * Initializes the transform for a given input length. + * + * @param length input signal size + */ + Dft(std::size_t length): + Fft(length) + { + } + + /** + * Destroys the transform object. + */ + ~Dft() + { + } + + virtual SpectrumType fft(const SampleType x[]); + virtual void ifft(SpectrumType spectrum, double x[]); + + private: + /** + * Complex unit (0.0 + 1.0j). + */ + static const ComplexType j; + }; +} + +#endif // DFT_H diff --git a/samplebrain/src/aquila/transform/Dft.o b/samplebrain/src/aquila/transform/Dft.o new file mode 100644 index 0000000..5c5ad38 Binary files /dev/null and b/samplebrain/src/aquila/transform/Dft.o differ diff --git a/samplebrain/src/aquila/transform/Fft.h b/samplebrain/src/aquila/transform/Fft.h new file mode 100644 index 0000000..12e1b9b --- /dev/null +++ b/samplebrain/src/aquila/transform/Fft.h @@ -0,0 +1,92 @@ +/** + * @file Fft.h + * + * An interface for FFT calculation classes. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef FFT_H +#define FFT_H + +#include "../global.h" +#include + +namespace Aquila +{ + /** + * An interface for FFT calculation classes. + * + * The abstract interface for FFT algorithm allows switching between + * implementations at runtime, or selecting a most effective implementation + * for the current platform. + * + * The FFT objects are not intended to be copied. + * + * Some of FFT implementations prepare a "plan" or create a coefficient + * cache only once, and then run the transform on multiple sets of data. + * Aquila expresses this approach by defining a meaningful constructor + * for the base FFT interface. A derived class should calculate the + * plan once - in the constructor (based on FFT length). Later calls + * to fft() / ifft() should reuse the already created plan/cache. + */ + class AQUILA_EXPORT Fft + { + public: + /** + * Initializes the transform for a given input length. + * + * @param length input signal size (usually a power of 2) + */ + Fft(std::size_t length): N(length) + { + } + + /** + * Destroys the transform object - does nothing. + * + * As the derived classes may perform some deinitialization in + * their destructors, it must be declared as virtual. + */ + virtual ~Fft() + { + } + + /** + * Applies the forward FFT transform to the signal. + * + * @param x input signal + * @return calculated spectrum + */ + virtual SpectrumType fft(const SampleType x[]) = 0; + + /** + * Applies the inverse FFT transform to the spectrum. + * + * @param spectrum input spectrum + * @param x output signal + */ + virtual void ifft(SpectrumType spectrum, double x[]) = 0; + + protected: + /** + * Signal and spectrum length. + */ + std::size_t N; + + private: + Fft( const Fft& ); + const Fft& operator=( const Fft& ); + }; +} + +#endif // FFT_H diff --git a/samplebrain/src/aquila/transform/FftFactory.cpp b/samplebrain/src/aquila/transform/FftFactory.cpp new file mode 100644 index 0000000..9f91848 --- /dev/null +++ b/samplebrain/src/aquila/transform/FftFactory.cpp @@ -0,0 +1,42 @@ +/** + * @file FftFactory.cpp + * + * A factory class to manage the creation of FFT calculation objects. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "FftFactory.h" +#include "OouraFft.h" + +namespace Aquila +{ + /** + * Returns "the best possible" FFT object. + * + * And now for some clarification about what is "the best possible": + * This method will perhaps take into consideration some optimisation + * hints, as optimizing for size or for speed. These hints will affect + * the choice of FFT implementation - hidden from the caller who gets + * only a pointer to the base abstract Fft class. + * + * As of now, the fastest implementation in Aquila is using Ooura's + * mathematical packages, so this one is always returned. + * + * @param length FFT length (number of samples) + * @return the FFT object (wrapped in a shared_ptr) + */ + std::shared_ptr FftFactory::getFft(std::size_t length) + { + return std::shared_ptr(new OouraFft(length)); + } +} diff --git a/samplebrain/src/aquila/transform/FftFactory.dep b/samplebrain/src/aquila/transform/FftFactory.dep new file mode 100644 index 0000000..f873d25 --- /dev/null +++ b/samplebrain/src/aquila/transform/FftFactory.dep @@ -0,0 +1,3 @@ +FftFactory.o: src/aquila/transform/FftFactory.cpp \ + src/aquila/transform/FftFactory.h src/aquila/transform/../global.h \ + src/aquila/transform/Fft.h src/aquila/transform/OouraFft.h diff --git a/samplebrain/src/aquila/transform/FftFactory.h b/samplebrain/src/aquila/transform/FftFactory.h new file mode 100644 index 0000000..64cb127 --- /dev/null +++ b/samplebrain/src/aquila/transform/FftFactory.h @@ -0,0 +1,38 @@ +/** + * @file FftFactory.h + * + * A factory class to manage the creation of FFT calculation objects. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef FFTFACTORY_H +#define FFTFACTORY_H + +#include "../global.h" +#include "Fft.h" +#include +#include + +namespace Aquila +{ + /** + * A factory class to manage the creation of FFT calculation objects. + */ + class AQUILA_EXPORT FftFactory + { + public: + static std::shared_ptr getFft(std::size_t length); + }; +} + +#endif // FFTFACTORY_H diff --git a/samplebrain/src/aquila/transform/FftFactory.o b/samplebrain/src/aquila/transform/FftFactory.o new file mode 100644 index 0000000..f77921e Binary files /dev/null and b/samplebrain/src/aquila/transform/FftFactory.o differ diff --git a/samplebrain/src/aquila/transform/Mfcc.cpp b/samplebrain/src/aquila/transform/Mfcc.cpp new file mode 100644 index 0000000..ac112f7 --- /dev/null +++ b/samplebrain/src/aquila/transform/Mfcc.cpp @@ -0,0 +1,43 @@ +/** + * @file Mfcc.cpp + * + * Calculation of MFCC signal features. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "Mfcc.h" +#include "Dct.h" +#include "../source/SignalSource.h" +#include "../filter/MelFilterBank.h" + +namespace Aquila +{ + /** + * Calculates a set of MFCC features from a given source. + * + * @param source input signal + * @param numFeatures how many features to calculate + * @return vector of MFCC features of length numFeatures + */ + std::vector Mfcc::calculate(const SignalSource &source, + std::size_t numFeatures) + { + auto spectrum = m_fft->fft(source.toArray()); + + Aquila::MelFilterBank bank(source.getSampleFrequency(), m_inputSize); + auto filterOutput = bank.applyAll(spectrum); + + Aquila::Dct dct; + return dct.dct(filterOutput, numFeatures); + } +} diff --git a/samplebrain/src/aquila/transform/Mfcc.h b/samplebrain/src/aquila/transform/Mfcc.h new file mode 100644 index 0000000..09ab47b --- /dev/null +++ b/samplebrain/src/aquila/transform/Mfcc.h @@ -0,0 +1,83 @@ +/** + * @file Mfcc.h + * + * Calculation of MFCC signal features. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef MFCC_H +#define MFCC_H + +#include "../global.h" +#include "FftFactory.h" +#include +#include + +namespace Aquila +{ + class SignalSource; + + /** + * The Mfcc class implements calculation of MFCC features from input signal. + * + * MFCC coefficients are commonly used in speech recognition. The common + * workflow is to split input signal in frames of equal length and apply + * MFCC calculation to each frame individually. Hence a few assumptions + * were made here: + * + * - a single Mfcc instance can be used only to process signals of equal + * length, for example consecutive frames + * - if you need to handle signals of various lengths, just create new + * Mfcc object per each signal source + * + * The code below is a simplest possible example of how to calculate MFCC + * for each frame of input signal. + * + * FramesCollection frames(data, FRAME_SIZE); + * Mfcc mfcc(FRAME_SIZE); + * for (Frame& frame : frames) { + * auto mfccValues = mfcc.calculate(frame); + * // do something with the calculated values + * } + * + */ + class AQUILA_EXPORT Mfcc + { + public: + /** + * Constructor creates the FFT object to reuse between calculations. + * + * @param inputSize input length (common to all inputs) + */ + Mfcc(std::size_t inputSize): + m_inputSize(inputSize), m_fft(FftFactory::getFft(inputSize)) + { + } + + std::vector calculate(const SignalSource& source, + std::size_t numFeatures = 12); + + private: + /** + * Number of samples in each processed input. + */ + const std::size_t m_inputSize; + + /** + * FFT calculator. + */ + std::shared_ptr m_fft; + }; +} + +#endif // MFCC_H diff --git a/samplebrain/src/aquila/transform/OouraFft.cpp b/samplebrain/src/aquila/transform/OouraFft.cpp new file mode 100644 index 0000000..e149494 --- /dev/null +++ b/samplebrain/src/aquila/transform/OouraFft.cpp @@ -0,0 +1,110 @@ +/** + * @file OouraFft.cpp + * + * A wrapper for the FFT algorithm found in Ooura mathematical packages. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "OouraFft.h" +#include +#include +#include + +namespace Aquila +{ + /** + * Initializes the transform for a given input length. + * + * Prepares the work area for Ooura's algorithm. + * + * @param length input signal size (usually a power of 2) + */ + OouraFft::OouraFft(std::size_t length): + Fft(length), + // according to the description: "length of ip >= 2+sqrt(n)" + ip(new int[static_cast(2 + std::sqrt(static_cast(N)))]), + w(new double[N / 2]) + { + ip[0] = 0; + } + + /** + * Destroys the transform object and cleans up working area. + */ + OouraFft::~OouraFft() + { + delete [] w; + delete [] ip; + } + + /** + * Applies the transformation to the signal. + * + * @param x input signal + * @return calculated spectrum + */ + SpectrumType OouraFft::fft(const SampleType x[]) + { + static_assert( + sizeof(ComplexType[2]) == sizeof(double[4]), + "complex has the same memory layout as two consecutive doubles" + ); + // create a temporary storage array and copy input to even elements + // of the array (real values), leaving imaginary components at 0 + double* a = new double[2 * N]; + for (std::size_t i = 0; i < N; ++i) + { + a[2 * i] = x[i]; + a[2 * i + 1] = 0.0; + } + + // let's call the C function from Ooura's package + cdft(2*N, -1, a, ip, w); + + // convert the array back to complex values and return as vector + ComplexType* tmpPtr = reinterpret_cast(a); + SpectrumType spectrum(tmpPtr, tmpPtr + N); + delete [] a; + + return spectrum; + } + + /** + * Applies the inverse transform to the spectrum. + * + * @param spectrum input spectrum + * @param x output signal + */ + void OouraFft::ifft(SpectrumType spectrum, double x[]) + { + static_assert( + sizeof(ComplexType[2]) == sizeof(double[4]), + "complex has the same memory layout as two consecutive doubles" + ); + // interpret the vector as consecutive pairs of doubles (re,im) + // and copy to input/output array + double* tmpPtr = reinterpret_cast(&spectrum[0]); + double* a = new double[2 * N]; + std::copy(tmpPtr, tmpPtr + 2 * N, a); + + // Ooura's function + cdft(2*N, 1, a, ip, w); + + // copy the data to the double array and scale it + for (std::size_t i = 0; i < N; ++i) + { + x[i] = a[2 * i] / static_cast(N); + } + delete [] a; + } +} diff --git a/samplebrain/src/aquila/transform/OouraFft.dep b/samplebrain/src/aquila/transform/OouraFft.dep new file mode 100644 index 0000000..c6d1758 --- /dev/null +++ b/samplebrain/src/aquila/transform/OouraFft.dep @@ -0,0 +1,3 @@ +OouraFft.o: src/aquila/transform/OouraFft.cpp \ + src/aquila/transform/OouraFft.h src/aquila/transform/Fft.h \ + src/aquila/transform/../global.h diff --git a/samplebrain/src/aquila/transform/OouraFft.h b/samplebrain/src/aquila/transform/OouraFft.h new file mode 100644 index 0000000..23406a3 --- /dev/null +++ b/samplebrain/src/aquila/transform/OouraFft.h @@ -0,0 +1,55 @@ +/** + * @file OouraFft.h + * + * A wrapper for the FFT algorithm found in Ooura mathematical packages. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef OOURAFFT_H +#define OOURAFFT_H + +#include "Fft.h" + +extern "C" { + void cdft(int, int, double *, int *, double *); + void rdft(int, int, double *, int *, double *); +} + +namespace Aquila +{ + /** + * A wrapper for the FFT algorithm found in Ooura mathematical packages. + */ + class AQUILA_EXPORT OouraFft : public Fft + { + public: + OouraFft(std::size_t length); + ~OouraFft(); + + virtual SpectrumType fft(const SampleType x[]); + virtual void ifft(SpectrumType spectrum, double x[]); + + private: + /** + * Work area for bit reversal. + */ + int* ip; + + /** + * Cos/sin table. + */ + double* w; + }; +} + +#endif // OOURAFFT_H diff --git a/samplebrain/src/aquila/transform/OouraFft.o b/samplebrain/src/aquila/transform/OouraFft.o new file mode 100644 index 0000000..98fc938 Binary files /dev/null and b/samplebrain/src/aquila/transform/OouraFft.o differ diff --git a/samplebrain/src/aquila/transform/Spectrogram.cpp b/samplebrain/src/aquila/transform/Spectrogram.cpp new file mode 100644 index 0000000..3298a80 --- /dev/null +++ b/samplebrain/src/aquila/transform/Spectrogram.cpp @@ -0,0 +1,44 @@ +/** + * @file Spectrogram.cpp + * + * Spectrogram calculation. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "Spectrogram.h" +#include "FftFactory.h" +#include "../source/Frame.h" +#include "../source/FramesCollection.h" + +namespace Aquila +{ + /** + * Creates the spectrogram from a collection of signal frames. + * + * Calculates frame spectra immediately after initialization. + * + * @param frames input frames + */ + Spectrogram::Spectrogram(FramesCollection& frames): + m_frameCount(frames.count()), + m_spectrumSize(frames.getSamplesPerFrame()), + m_fft(FftFactory::getFft(m_spectrumSize)), + m_data(new SpectrogramDataType(m_frameCount)) + { + std::size_t i = 0; + for (auto it = frames.begin(); it != frames.end(); ++it, ++i) + { + (*m_data)[i] = m_fft->fft(it->toArray()); + } + } +} diff --git a/samplebrain/src/aquila/transform/Spectrogram.h b/samplebrain/src/aquila/transform/Spectrogram.h new file mode 100644 index 0000000..84fcb4f --- /dev/null +++ b/samplebrain/src/aquila/transform/Spectrogram.h @@ -0,0 +1,102 @@ +/** + * @file Spectrogram.h + * + * Spectrogram calculation. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef SPECTROGRAM_H +#define SPECTROGRAM_H + +#include "../global.h" +#include +#include +#include + + +namespace Aquila +{ + class Fft; + class FramesCollection; + + /** + * Spectrogram class. + * + * @todo copy constructor, safe point access, code cleanup + */ + class AQUILA_EXPORT Spectrogram + { + public: + Spectrogram(FramesCollection& frames); + + /** + * Returns number of frames (spectrogram width). + * + * @return frame count + */ + std::size_t getFrameCount() const + { + return m_frameCount; + } + + /** + * Returns spectrum size (spectrogram height). + * + * @return spectrum size + */ + std::size_t getSpectrumSize() const + { + return m_spectrumSize; + } + + /** + * Returns a complex data point value at given coordinates. + * + * @param frame frame number (the x coordinate) + * @param peak spectral peak number (the y coordinate) + * @return spectrum value at given frame and peak number + */ + ComplexType getPoint(std::size_t frame, std::size_t peak) const + { + return (*m_data)[frame][peak]; + } + + private: + /** + * Spectrogram data type used for internal storage. + */ + typedef std::vector SpectrogramDataType; + + /** + * Frame count (width of the spectrogram). + */ + std::size_t m_frameCount; + + /** + * Spectrum size (height of the spectrogram). + */ + std::size_t m_spectrumSize; + + /** + * A shared pointer to FFT algorithm class. + */ + std::shared_ptr m_fft; + + /** + * A shared pointer to spectrogram data. + */ + std::shared_ptr m_data; + }; +} + +#endif // SPECTROGRAM_H diff --git a/samplebrain/src/aquila/wrappers/SoundBufferAdapter.cpp b/samplebrain/src/aquila/wrappers/SoundBufferAdapter.cpp new file mode 100644 index 0000000..93103bd --- /dev/null +++ b/samplebrain/src/aquila/wrappers/SoundBufferAdapter.cpp @@ -0,0 +1,88 @@ +/** + * @file SoundBufferAdapter.cpp + * + * A wrapper around SignalSource to use as a sound buffer in SFML. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#include "SoundBufferAdapter.h" +#include "../source/SignalSource.h" +#include +#include + +namespace Aquila +{ + /** + * Creates the buffer with no initial data. + */ + SoundBufferAdapter::SoundBufferAdapter(): + SoundBuffer() + { + } + + /** + * Copy constructor. + * + * @param other buffer instance to copy from + */ + SoundBufferAdapter::SoundBufferAdapter(const SoundBufferAdapter &other): + SoundBuffer(other) + { + } + + /** + * Creates the buffer with initial data provided by signal source. + * + * @param source signal source + */ + SoundBufferAdapter::SoundBufferAdapter(const SignalSource &source): + SoundBuffer() + { + loadFromSignalSource(source); + } + + /** + * Destructor - does nothing by itself. + * + * Relies on virtual call to the destructor of the parent class. + */ + SoundBufferAdapter::~SoundBufferAdapter() + { + } + + /** + * Loads sound data from an instance of SignalSource-subclass. + * + * Data read from source are converted to SFML-compatible sample array + * and loaded into the buffer. + * + * Name capitalized for consistency with SFML coding style. + * + * @todo get rid of copying data around, let's come up with some better way + * + * @param source signal source + * @return true if successfully loaded + */ + bool SoundBufferAdapter::loadFromSignalSource(const SignalSource &source) + { + sf::Int16* samples = new sf::Int16[source.getSamplesCount()]; + std::copy(source.begin(), source.end(), samples); + bool result = loadFromSamples(samples, + source.getSamplesCount(), + 1, + static_cast(source.getSampleFrequency())); + delete [] samples; + + return result; + } +} diff --git a/samplebrain/src/aquila/wrappers/SoundBufferAdapter.h b/samplebrain/src/aquila/wrappers/SoundBufferAdapter.h new file mode 100644 index 0000000..eeab70c --- /dev/null +++ b/samplebrain/src/aquila/wrappers/SoundBufferAdapter.h @@ -0,0 +1,43 @@ +/** + * @file SoundBufferAdapter.h + * + * A wrapper around SignalSource to use as a sound buffer in SFML. + * + * This file is part of the Aquila DSP library. + * Aquila is free software, licensed under the MIT/X11 License. A copy of + * the license is provided with the library in the LICENSE file. + * + * @package Aquila + * @version 3.0.0-dev + * @author Zbigniew Siciarz + * @date 2007-2014 + * @license http://www.opensource.org/licenses/mit-license.php MIT + * @since 3.0.0 + */ + +#ifndef SOUNDBUFFERADAPTER_H +#define SOUNDBUFFERADAPTER_H + +#include "../global.h" +#include + +namespace Aquila +{ + class SignalSource; + + /** + * A wrapper around SignalSource to use as a sound buffer in SFML. + */ + class AQUILA_EXPORT SoundBufferAdapter : public sf::SoundBuffer + { + public: + SoundBufferAdapter(); + SoundBufferAdapter(const SoundBufferAdapter& other); + SoundBufferAdapter(const SignalSource& source); + ~SoundBufferAdapter(); + + bool loadFromSignalSource(const SignalSource& source); + }; +} + +#endif // SOUNDBUFFERADAPTER_H diff --git a/samplebrain/src/brain.cpp b/samplebrain/src/brain.cpp index 629c6cf..46a6022 100644 --- a/samplebrain/src/brain.cpp +++ b/samplebrain/src/brain.cpp @@ -50,7 +50,8 @@ void brain::init(u32 block_size, u32 overlap, bool ditchpcm) { void brain::chop_and_add(const sample &s, u32 block_size, u32 overlap, bool ditchpcm) { u32 pos=0; while (pos+block_size-1::iterator i=m_blocks.begin(); i!=m_blocks.end(); ++i) { - cerr< m_blocks; @@ -36,6 +41,7 @@ private: u32 m_block_size; u32 m_overlap; + }; #endif diff --git a/samplebrain/src/brain_block.cpp b/samplebrain/src/brain_block.cpp index ea7cfa4..4235c96 100644 --- a/samplebrain/src/brain_block.cpp +++ b/samplebrain/src/brain_block.cpp @@ -7,8 +7,9 @@ using namespace spiralcore; FFT *brain_block::m_fftw; +Aquila::Mfcc *brain_block::m_mfcc_proc; -static const int MFCC_FILTERS=48; +static const int MFCC_FILTERS=12; void enveloper(sample &s, u32 start, u32 end) { for(u32 i=0; iimpulse2freq(m_pcm.get_non_const_buffer(), - m_fft.get_non_const_buffer()); + m_fftw->impulse2freq(m_pcm.get_non_const_buffer()); - if (m_block_size>30) m_fft.crop_to(30); + std::vector> mfspec; + + for (u32 i=0; im_spectrum[i][0]; + + mfspec.push_back(std::complex(m_fftw->m_spectrum[i][0], + m_fftw->m_spectrum[i][1])); + } + + if (m_block_size>100) m_fft.crop_to(100); if (ditchpcm) m_pcm.clear(); -// for (u32 i=0; i m = m_mfcc_proc->calculate(mfspec,MFCC_FILTERS); + + for (u32 i=0; im_length!=block_size) { if (m_fftw == NULL) delete m_fftw; m_fftw = new FFT(block_size); + if (m_mfcc_proc == NULL) delete m_mfcc_proc; + m_mfcc_proc = new Aquila::Mfcc(block_size); } } double brain_block::compare(const brain_block &other, float ratio) const { - double acc=0; - // just mfcc - //if (ratio==1) - { - for (u32 i=0; i +#include using namespace spiralcore; @@ -6,35 +7,19 @@ static const int MAX_FFT_LENGTH = 4096; FFT::FFT(int length) : m_length(length), -#ifndef __FFTWFLOAT__ m_in(new double[length]), -#else - m_in(new float[length]), -#endif - -#ifndef __FFTWFLOAT__ m_spectrum(new fftw_complex[length]) { m_plan = fftw_plan_dft_r2c_1d(m_length, m_in, m_spectrum, FFTW_ESTIMATE); } -#else -m_spectrum(new fftwf_complex[length]) -{ - m_plan = fftwf_plan_dft_r2c_1d(m_length, m_in, m_spectrum, FFTW_ESTIMATE); -} -#endif FFT::~FFT() { delete[] m_in; -#ifndef __FFTWFLOAT__ fftw_destroy_plan(m_plan); -#else - fftwf_destroy_plan(m_plan); -#endif } -void FFT::impulse2freq(float *imp, float *out) +void FFT::impulse2freq(float *imp) { unsigned int i; @@ -43,19 +28,5 @@ void FFT::impulse2freq(float *imp, float *out) m_in[i] = imp[i]; } -#ifndef __FFTWFLOAT__ - fftw_execute(m_plan); -#else - fftwf_execute(m_plan); -#endif - - for (i=0; i