mfccs from aquila

This commit is contained in:
Dave Griffiths 2015-07-09 16:16:36 +01:00
parent 299af9442d
commit c21e3f4d21
102 changed files with 6587 additions and 81 deletions

View File

@ -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@

View File

@ -4,3 +4,4 @@ make distclean
autoheader
# build configure
autoconf configure.ac > configure
./configure CXX=g++-4.7

View File

@ -0,0 +1 @@
dave@fulmar.4739:1436183149

View File

@ -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 <stdexcept>
#include <string>
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

View File

@ -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

View File

@ -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

View File

@ -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 <algorithm>
#include <utility>
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<std::size_t>(N * minFreq / m_sampleFrequency);
std::size_t maxPos = static_cast<std::size_t>(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);
}
}
}
}
}

View File

@ -0,0 +1,2 @@
MelFilter.o: src/aquila/filter/MelFilter.cpp \
src/aquila/filter/MelFilter.h src/aquila/filter/../global.h

View File

@ -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 <cstddef>
#include <cmath>
#include <vector>
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<double> m_spectrum;
void generateFilterSpectrum(FrequencyType minFreq,
FrequencyType centerFreq,
FrequencyType maxFreq, std::size_t N);
};
}
#endif // MELFILTER_H

Binary file not shown.

View File

@ -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<double> MelFilterBank::applyAll(const SpectrumType& frameSpectrum) const
{
std::vector<double> output(size(), 0.0);
for (std::size_t i = 0; i < size(); ++i)
{
output[i] = m_filters[i].apply(frameSpectrum);
}
return output;
}
}

View File

@ -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

View File

@ -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 <cstddef>
#include <vector>
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<double> 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<MelFilter> 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

Binary file not shown.

View File

@ -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 <algorithm>
#include <cmath>
#include <cstdlib>
#include <functional>
namespace Aquila
{
/**
* Converts the value to decibels (assuming reference value equal to 1).
*
* @param value input value
* @return value in dB
*/
template<typename Numeric>
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<typename Numeric>
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<typename Numeric>
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<double>(RAND_MAX);
}
/**
* Checks if n is an exact power of 2.
*/
template<typename Integer>
AQUILA_EXPORT inline bool isPowerOf2(Integer n)
{
return (n > 0) && ((n & (n - 1)) == 0);
}
/**
* Returns the smallest power of 2 greater than n.
*/
template<typename Integer>
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<double(const std::vector<double>&,
const std::vector<double>&)> 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<double>& v1,
const std::vector<double>& 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<double>& v1,
const std::vector<double>& 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<double>& v1,
const std::vector<double>& 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

View File

@ -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 <complex>
#include <vector>
#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<double> ComplexType;
/**
* Spectrum type - a vector of complex values.
*/
typedef std::vector<ComplexType> SpectrumType;
}
#endif // GLOBAL_H

View File

@ -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

View File

@ -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;
}
}

View File

@ -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 <cstddef>
#include <vector>
namespace Aquila
{
/**
* Type of compared data - vectors of features, which themselves are
* vectors of doubles.
*/
typedef std::vector<std::vector<double>> DtwDataType;
/**
* Type of DTW point array.
*/
typedef std::vector<std::vector<DtwPoint>> DtwPointsArrayType;
/**
* Lowest-cost path is a vector of points.
*/
typedef std::vector<DtwPoint> 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

View File

@ -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 <cstddef>
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

View File

@ -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

View File

@ -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;
}
}

View File

@ -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 <algorithm>
#include <cstddef>
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<std::ptrdiff_t>(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

View File

@ -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<unsigned int>(
source.getSampleFrequency() * frameDuration / 1000.0
);
unsigned int samplesPerOverlap = static_cast<unsigned int>(
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();
}
}

View File

@ -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 <algorithm>
#include <cstddef>
#include <functional>
#include <vector>
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<Frame> 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 <typename ResultType>
std::vector<ResultType> apply(
std::function<ResultType (const SignalSource&)> f) const
{
std::vector<ResultType> 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

View File

@ -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 <algorithm>
#include <fstream>
#include <iterator>
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<SampleType>(fs),
std::istream_iterator<SampleType>(),
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<SampleType>(fs, "\n"));
fs.close();
}
}

View File

@ -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 <string>
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

View File

@ -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 <algorithm>
#include <cstddef>
#include <fstream>
#include <string>
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 <typename Numeric = SampleType>
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

View File

@ -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 <algorithm>
#include <cmath>
#include <numeric>
#include <utility>
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));
}
}

View File

@ -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 <cstddef>
#include <iterator>
#include <utility>
#include <vector>
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 <typename Numeric>
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<SampleType>& 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<SampleType>&& 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<std::forward_iterator_tag, int>
{
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<SampleType> 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

View File

@ -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<unsigned int>(m_header.WaveSize /
static_cast<double>(m_header.BytesPerSec) * 1000);
}
}

View File

@ -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 <cstddef>
#include <cstdint>
#include <string>
#include <vector>
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

View File

@ -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 <cstdint>
#include <cstring>
#include <fstream>
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<std::uint32_t>(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<short>(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<unsigned char>(source.sample(2 * i) + 128);
unsigned char sample2 = static_cast<unsigned char>(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;
}
}

View File

@ -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 <cstddef>
#include <string>
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

View File

@ -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)
{
}
}

View File

@ -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 <cstddef>
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

View File

@ -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;
}
}

View File

@ -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

View File

@ -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 <cmath>
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<double>(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
);
}
}
}

View File

@ -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

View File

@ -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 <cmath>
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<std::size_t>(
m_sampleFrequency / static_cast<double>(m_frequency));
std::size_t positiveLength = static_cast<std::size_t>(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);
}
}
}

View File

@ -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

View File

@ -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 <cmath>
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;
}
}
}

View File

@ -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

View File

@ -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);
}
}
}

View File

@ -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

View File

@ -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 <cmath>
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))
);
}
}
}

View File

@ -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 <cstddef>
namespace Aquila
{
/**
* Barlett (triangular) window.
*/
class AQUILA_EXPORT BarlettWindow : public SignalSource
{
public:
BarlettWindow(std::size_t size);
};
}
#endif // BARLETTWINDOW_H

View File

@ -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 <cmath>
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))
);
}
}
}

View File

@ -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 <cstddef>
namespace Aquila
{
/**
* Blackman window.
*/
class AQUILA_EXPORT BlackmanWindow : public SignalSource
{
public:
BlackmanWindow(std::size_t size);
};
}
#endif // BLACKMANWINDOW_H

View File

@ -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 <cmath>
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))
);
}
}
}

View File

@ -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 <cstddef>
namespace Aquila
{
/**
* Flat-top window.
*/
class AQUILA_EXPORT FlattopWindow : public SignalSource
{
public:
FlattopWindow(std::size_t size);
};
}
#endif // FLATTOPWINDOW_H

View File

@ -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 <cmath>
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))
);
}
}
}

View File

@ -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 <cstddef>
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

View File

@ -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 <cmath>
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))
);
}
}
}

View File

@ -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 <cstddef>
namespace Aquila
{
/**
* Hamming window.
*/
class AQUILA_EXPORT HammingWindow : public SignalSource
{
public:
HammingWindow(std::size_t size);
};
}
#endif // HAMMINGWINDOW_H

View File

@ -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 <cmath>
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)))
);
}
}
}

View File

@ -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 <cstddef>
namespace Aquila
{
/**
* Hann window.
*/
class AQUILA_EXPORT HannWindow : public SignalSource
{
public:
HannWindow(std::size_t size);
};
}
#endif // HANNWINDOW_H

View File

@ -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 <cstddef>
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

View File

@ -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

View File

@ -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 <SFML/Audio.hpp>
#include <SFML/System.hpp>
#include <algorithm>
#include <cstddef>
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<std::size_t>(m_sampleFrequency / frequency);
std::size_t totalSamples = static_cast<std::size_t>(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;
}
}

View File

@ -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

View File

@ -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 <SFML/Audio.hpp>
#include <SFML/System.hpp>
#include <cstddef>
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<std::size_t>(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));
}
}

View File

@ -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

View File

@ -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 <SFML/System/Sleep.hpp>
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);
}
}
}

View File

@ -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 <SFML/System.hpp>
#include <SFML/Audio.hpp>
#include <map>
#include <string>
namespace Aquila
{
/**
* Type of the mapping from note names to frequencies.
*/
typedef std::map<std::string, FrequencyType> 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

View File

@ -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

View File

@ -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<double> absSpectrum(halfLength);
for (std::size_t i = 0; i < halfLength; ++i)
{
absSpectrum[i] = std::abs(spectrum[i]);
}
plot(absSpectrum);
}
}

View File

@ -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 <algorithm>
#include <cmath>
#include <cstddef>
#include <iostream>
#include <string>
#include <vector>
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<typename Numeric>
void plot(const std::vector<Numeric>& 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<typename Numeric>
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<std::vector<char>> 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 <typename Iterator>
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::size_t>(
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

View File

@ -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

View File

@ -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 <algorithm>
#include <cmath>
#include <iostream>
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<unsigned int>(
std::log(static_cast<double>(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<unsigned int>(
std::log(static_cast<double>(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<double>(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<double>(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;
}
}
}

View File

@ -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 <map>
#include <utility>
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<fftWiCacheKeyType, ComplexType**> fftWiCacheType;
/**
* Twiddle factor cache - implemented as a map.
*/
fftWiCacheType fftWiCache;
ComplexType** getCachedFftWi(unsigned int numStages);
void clearFftWiCache();
};
}
#endif // AQUILAFFT_H

View File

@ -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 <algorithm>
#include <cmath>
#include <iterator>
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<double> Dct::dct(const std::vector<double>& data, std::size_t outputLength)
{
// zero-initialize output vector
std::vector<double> 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;
}
}
}

View File

@ -0,0 +1,2 @@
Dct.o: src/aquila/transform/Dct.cpp src/aquila/transform/Dct.h \
src/aquila/transform/../global.h

View File

@ -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 <cstddef>
#include <map>
#include <utility>
#include <vector>
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<double> dct(const std::vector<double>& data, std::size_t outputLength);
private:
/**
* Key type for the cache, using input and output length.
*/
typedef std::pair<std::size_t, std::size_t> cosineCacheKeyType;
/**
* Cache type.
*/
typedef std::map<cosineCacheKeyType, double**> cosineCacheType;
/**
* Cache object, implemented as a map.
*/
cosineCacheType cosineCache;
double** getCachedCosines(std::size_t inputLength, std::size_t outputLength);
void clearCosineCache();
};
}
#endif // DCT_H

Binary file not shown.

View File

@ -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 <algorithm>
#include <cmath>
#include <complex>
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<double>(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<double>(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<int>(n * k));
}
x[k] = std::abs(sum) / static_cast<double>(N);
}
}
}

View File

@ -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

View File

@ -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

Binary file not shown.

View File

@ -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 <cstddef>
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

View File

@ -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<Fft> FftFactory::getFft(std::size_t length)
{
return std::shared_ptr<Fft>(new OouraFft(length));
}
}

View File

@ -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

View File

@ -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 <cstddef>
#include <memory>
namespace Aquila
{
/**
* A factory class to manage the creation of FFT calculation objects.
*/
class AQUILA_EXPORT FftFactory
{
public:
static std::shared_ptr<Fft> getFft(std::size_t length);
};
}
#endif // FFTFACTORY_H

Binary file not shown.

View File

@ -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<double> 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);
}
}

View File

@ -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 <cstddef>
#include <vector>
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<double> 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<Fft> m_fft;
};
}
#endif // MFCC_H

View File

@ -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 <algorithm>
#include <cmath>
#include <cstddef>
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<std::size_t>(2 + std::sqrt(static_cast<double>(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<double> 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<ComplexType*>(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<double> 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<double*>(&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<double>(N);
}
delete [] a;
}
}

View File

@ -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

View File

@ -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

Binary file not shown.

View File

@ -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());
}
}
}

View File

@ -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 <cstddef>
#include <memory>
#include <vector>
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<SpectrumType> 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<Fft> m_fft;
/**
* A shared pointer to spectrogram data.
*/
std::shared_ptr<SpectrogramDataType> m_data;
};
}
#endif // SPECTROGRAM_H

View File

@ -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 <SFML/System.hpp>
#include <algorithm>
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<unsigned int>(source.getSampleFrequency()));
delete [] samples;
return result;
}
}

View File

@ -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 <SFML/Audio.hpp>
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

View File

@ -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<s.get_length()) {
cerr<<pos/(float)s.get_length()*100<<endl;
cerr<<'\r';
cerr<<"adding: "<<pos/(float)s.get_length()*100;
sample region;
s.get_region(region,pos,pos+block_size-1);
m_blocks.push_back(brain_block("",region,44100,ditchpcm));
@ -62,6 +63,9 @@ const sample &brain::get_block_pcm(u32 index) const {
return m_blocks[index].get_pcm();
}
const brain_block &brain::get_block(u32 index) const {
return m_blocks[index];
}
// returns index to block
u32 brain::search(const brain_block &target, float ratio) const {
@ -86,11 +90,14 @@ void brain::resynth(const string &filename, const brain &other, float ratio){
out.zero();
u32 pos = 0;
u32 count = 0;
cerr<<other.m_blocks.size()<<" brain blocks..."<<endl;
for (vector<brain_block>::iterator i=m_blocks.begin(); i!=m_blocks.end(); ++i) {
cerr<<count/float(m_blocks.size())*100<<endl;
cerr<<'\r';
cerr<<"searching: "<<count/float(m_blocks.size())*100;
u32 index = other.search(*i,ratio);
cerr<<index<<endl;
out.mix(other.get_block_pcm(index),pos);
//cerr<<index<<endl;
out.mul_mix(other.get_block_pcm(index),pos,0.2);
if (count%1000==0) {
save_sample(filename,out);
@ -100,9 +107,6 @@ void brain::resynth(const string &filename, const brain &other, float ratio){
pos += (m_block_size-m_overlap);
}
save_sample(filename,out);
cerr<<m_blocks.size()<<" brain blocks..."<<endl;
}

View File

@ -24,11 +24,16 @@ public:
void resynth(const std::string &filename, const brain &other, float ratio);
const sample &get_block_pcm(u32 index) const;
const brain_block &get_block(u32 index) const;
const u32 &get_block_size() const { return m_block_size; }
const u32 &get_overlap() const { return m_overlap; }
u32 search(const brain_block &target, float ratio) const;
static bool unit_test();
private:
u32 search(const brain_block &target, float ratio) const;
void chop_and_add(const sample &s, u32 block_size, u32 overlap, bool ditchpcm=false);
vector<brain_block> m_blocks;
@ -36,6 +41,7 @@ private:
u32 m_block_size;
u32 m_overlap;
};
#endif

View File

@ -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; i<start; ++i) {
@ -28,18 +29,32 @@ brain_block::brain_block(const string &filename, const sample &pcm, u32 rate, bo
m_orig_filename(filename)
{
init_fft(m_pcm.get_length());
assert(m_mfcc_proc!=NULL);
assert(m_fftw!=NULL);
enveloper(m_pcm,50,50);
m_fftw->impulse2freq(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<std::complex<double>> mfspec;
for (u32 i=0; i<m_block_size; ++i) {
m_fft[i]=m_fftw->m_spectrum[i][0];
mfspec.push_back(std::complex<double>(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<MFCC_FILTERS; i++) {
// m_mfcc[i] = GetCoefficient(m_fft.get_non_const_buffer(), rate, MFCC_FILTERS, m_block_size, i);
// }
// calculate mfcc
std::vector<double> m = m_mfcc_proc->calculate(mfspec,MFCC_FILTERS);
for (u32 i=0; i<MFCC_FILTERS; ++i) {
m_mfcc[i] = m[i];
}
}
void brain_block::init_fft(u32 block_size)
@ -47,28 +62,44 @@ void brain_block::init_fft(u32 block_size)
if (m_fftw == NULL || m_fftw->m_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<m_fft.get_length(); ++i) {
acc+=(m_fft[i]-other.m_fft[i]) * (m_fft[i]-other.m_fft[i]);
}
double mfcc_acc=0;
double fft_acc=0;
//for (u32 i=0; i<MFCC_FILTERS; ++i) {
// acc+=(m_mfcc[i]-other.m_mfcc[i]) * (m_mfcc[i]-other.m_mfcc[i]);
//}
if (ratio==0) {
for (u32 i=0; i<m_fft.get_length(); ++i) {
fft_acc+=(m_fft[i]-other.m_fft[i]) * (m_fft[i]-other.m_fft[i]);
}
return fft_acc/(float)m_fft.get_length();
}
return acc;
if (ratio==1) {
for (u32 i=0; i<MFCC_FILTERS; ++i) {
mfcc_acc+=(m_mfcc[i]-other.m_mfcc[i]) * (m_mfcc[i]-other.m_mfcc[i]);
}
return mfcc_acc/(float)MFCC_FILTERS;
}
// calculate both
for (u32 i=0; i<m_fft.get_length(); ++i) {
fft_acc+=(m_fft[i]-other.m_fft[i]) * (m_fft[i]-other.m_fft[i]);
}
for (u32 i=0; i<MFCC_FILTERS; ++i) {
mfcc_acc+=(m_mfcc[i]-other.m_mfcc[i]) * (m_mfcc[i]-other.m_mfcc[i]);
}
return (fft_acc/(float)m_fft.get_length())*(1-ratio) +
(mfcc_acc/(float)MFCC_FILTERS)*ratio;
}
bool brain_block::unit_test() {
sample data(20);
sample data(200);
for (u32 i=0; i<data.get_length(); i++) {
data[i]=i/(float)data.get_length();
}
@ -76,7 +107,7 @@ bool brain_block::unit_test() {
brain_block bb("test",data,44100);
assert(bb.m_pcm.get_length()==data.get_length());
assert(bb.m_fft.get_length()==data.get_length());
//assert(bb.m_fft.get_length()==data.get_length());
assert(bb.m_mfcc.get_length()==MFCC_FILTERS);
assert(bb.m_orig_filename==string("test"));
assert(bb.m_rate==44100);
@ -84,21 +115,24 @@ bool brain_block::unit_test() {
brain_block bb2("test",data,44100);
assert(bb.compare(bb2,1)==0);
assert(bb.compare(bb2,0)==0);
assert(bb.compare(bb2,0.5)==0);
sample data2(20);
sample data2(200);
for (u32 i=0; i<data.get_length(); i++) {
data[i]=i;
data[i]=i%10;
}
brain_block cpy("test",data,100);
{
brain_block bb3("test",data2,44100);
cerr<<bb.compare(bb3,1)<<endl;
assert(bb.compare(bb3,1)!=0);
assert(bb.compare(bb3,0)!=0);
assert(bb.compare(bb3,0.5)!=0);
cpy=bb3;
}
assert(cpy.m_pcm.get_length()==20);
assert(cpy.m_pcm.get_length()==200);
return true;
}

View File

@ -2,6 +2,7 @@
#include "jellyfish/fluxa/sample.h"
#include "jellyfish/core/types.h"
#include "fft.h"
#include "mfcc.h"
#ifndef BRAIN_BLOCK
#define BRAIN_BLOCK
@ -30,6 +31,7 @@ private:
u32 m_rate;
std::string m_orig_filename;
static FFT *m_fftw;
static Aquila::Mfcc *m_mfcc_proc;
};

View File

@ -1,4 +1,5 @@
#include <fft.h>
#include <jellyfish/core/types.h>
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<m_length; i++)
{
out[i] = m_spectrum[i][0];
}
}
void FFT::raw_impulse2freq()
{
fftw_execute(m_plan);
}

Some files were not shown because too many files have changed in this diff Show More