Commit 32ca64e1 authored by Votre Nom's avatar Votre Nom
Browse files

modenise C++

parent 3baf71c9
cmake_minimum_required(VERSION 2.6)
cmake_minimum_required(VERSION 3.13)
project(f_fftw_spectrum)
find_path(FFTW_INCLUDE_DIR fftw3.h)
find_library(FFTW_LIBRARIES NAMES fftw3f)
find_package(blc_channel REQUIRED)
find_package(blc_program REQUIRED)
find_package(blc REQUIRED)
add_definitions(${BL_DEFINITIONS})
include_directories(${BL_INCLUDE_DIRS} ${FFTW_INCLUDE_DIR})
add_executable(f_fftw_spectrum f_fftw_spectrum.cpp)
target_link_libraries(f_fftw_spectrum ${BL_LIBRARIES} ${FFTW_LIBRARIES})
target_link_options(f_fftw_spectrum PRIVATE -Wno-multichar)
target_link_libraries(f_fftw_spectrum PRIVATE blc rt ${FFTW_LIBRARIES})
#include "blc_text.h"
#include "blc_channel.h"
#include "blc_program.h"
#include <fftw3.h>
#include <unistd.h>
#include <math.h> //sqrtf
#include <pthread.h>
#include <memory>
using namespace std;
#define DEFAULT_OUTPUT_NAME "spectrum<pid>"
......@@ -13,31 +17,22 @@ struct timeval;
char const *display, *option_record;
int columns_nb=64, rows_nb=16;
static void init_output_channel(blc_channel *output_channel, char const *output_name, int output_length){
//Create a channel or reopen a previous one.
output_channel->create_or_open(output_name, BLC_CHANNEL_WRITE, 'FL32', 'NDEF', 1, output_length);
output_channel->publish();
blc_loop_try_add_waiting_semaphore(output_channel->sem_ack_data);
blc_loop_try_add_posting_semaphore(output_channel->sem_new_data);
}
int main(int argc, char**argv){
blc_channel input, output;
blc_channel input;
unique_ptr<blc_channel> output;
char const *input_name, *output_name, *spectrum_option, *period_str;
FILE *file=NULL;
FILE *file=nullptr;
int ret;
int output_length, period, pos=0;
size_t i;
fftwf_complex *intermediate;
fftwf_plan fftw_plan;
blc_program_add_option(&display, 'd', "display", NULL, "Display the result as text graph", NULL);
blc_program_add_option(&display, 'd', "display", nullptr, "Display the result as text graph", nullptr);
blc_program_add_option(&output_name, 'o', "output", "blc_channel|-", "channel containing the fast fourier transformation, or '-' for terminal output", DEFAULT_OUTPUT_NAME);
blc_program_add_option(&period_str, 'p', "period", "integer", "Period in ms (0 as fast as possible).", "0");
blc_program_add_option(&spectrum_option, 's', "spectrum", NULL, "Get the square value of the fftw signal", NULL);
blc_program_add_parameter(&input_name, "input blc_channel", 1, "Channel containing the sound data", NULL);
blc_program_add_option(&spectrum_option, 's', "spectrum", nullptr, "Get the square value of the fftw signal", nullptr);
blc_program_add_parameter(&input_name, "input blc_channel", 1, "Channel containing the sound data", nullptr);
blc_program_init(&argc, &argv, blc_quit);
SSCANF(1, period_str, "%d", &period);
......@@ -45,26 +40,34 @@ int main(int argc, char**argv){
input.open(input_name, BLC_CHANNEL_READ);
if (input.type!='FL32') EXIT_ON_CHANNEL_ERROR(&input,"Input must be of type 'FL32'");
output_length=input.total_length/2+1;
if (strcmp("-", output_name)==0){
output.init('FL32', 'NDEF', 1, output_length);
file=stdout;
}
else if (strcmp(output_name, DEFAULT_OUTPUT_NAME)==0) {
asprintf((char**)&output_name, ":spectrum%d", getpid());
init_output_channel(&output, output_name, output_length);
vector<size_t> dims = { static_cast<size_t>(output_length) };
if (strcmp(output_name, DEFAULT_OUTPUT_NAME)==0) {
asprintf((char**)&output_name, ":spectrum%d", getpid());
output=make_unique<blc_channel>(output_name, BLC_CHANNEL_WRITE, 'FL32', 'NDEF', dims);
output->publish();
blc_loop_try_add_waiting_semaphore(output->sem_ack_data);
blc_loop_try_add_posting_semaphore(output->sem_new_data);
}
else{
ret=sscanf(output_name, "%*[:/.^]%*[^/]%n", &pos);
if (pos==strlen(output_name))
ret=sscanf(output_name, "%*[:/.^]%*[^/]%n", &pos);
if (pos==strlen(output_name)){
output = make_unique<blc_channel> (output_name, BLC_CHANNEL_WRITE, 'FL32', 'NDEF', dims);
output->publish();
init_output_channel(&output, output_name, output_length); //C'est un blc_channel
else //C'est un format de fichier
{
if (strcmp(blc_get_filename_extension(output_name), "tsv")!=0) EXIT_ON_ERROR("Filename extension must be .tsv but your file is: '%s'.\nThis is to be sure you are not overwriting important file by mistake", output_name);
output.init('FL32', 'NDEF', 1, output_length);
SYSTEM_ERROR_CHECK(file=fopen(output_name,"w"), NULL, "opening '%s'", output_name);
}
blc_loop_try_add_waiting_semaphore(output->sem_ack_data);
blc_loop_try_add_posting_semaphore(output->sem_new_data);
}
else //C'est un format de fichier
{
EXIT_ON_ERROR("File not managed");
// if (strcmp(blc_get_filename_extension(output_name), "tsv")!=0) EXIT_ON_ERROR("Filename extension must be .tsv but your file is: '%s'.\nThis is to be sure you are not overwriting important file by mistake", output_name);
// output=make_unique<blc_channel>('FL32', 'NDEF', dims);
// SYSTEM_ERROR_CHECK(file=fopen(output_name,"w"), nullptr, "opening '%s'", output_name);
}
}
blc_loop_try_add_waiting_semaphore(input.sem_new_data);
......@@ -76,23 +79,23 @@ init_output_channel(&output, output_name, output_length); //C'est un blc_chann
BLC_COMMAND_LOOP(period*1000){
fftwf_execute(fftw_plan); /* repeat as needed */
FOR(i, output_length){
output.floats[i]=sqrtf(intermediate[i][0]*intermediate[i][0]+intermediate[i][1]*intermediate[i][1])/input.total_length;
output->floats[i]=sqrtf(intermediate[i][0]*intermediate[i][0]+intermediate[i][1]*intermediate[i][1])/input.total_length;
}
//if output is not a channel we export to standard output
if (file) {
// if (blc_loop_iteration && file==stdout && blc_output_terminal) blc_eprint_cursor_up(100);
output.fprint_tsv(file);
// if (blc_loop_iteration && file==stdout && blc_output_terminal) blc_eprint_cursor_up(100);
output->fprint_tsv(file);
}
if(display){
blc_fprint_float_graph(stderr, output.floats, output.total_length, "spectrum", columns_nb, rows_nb, 1, 0, "Frequency", "Intensity");
blc_fprint_float_graph(stderr, output->floats, output->total_length, "spectrum", columns_nb, rows_nb, 1, 0, "Frequency", "Intensity");
blc_eprint_cursor_up(rows_nb);
}
}
fftwf_destroy_plan(fftw_plan);
blc_eprint_cursor_down(rows_nb);
blc_eprint_cursor_down(rows_nb);
if (file) SYSTEM_ERROR_CHECK(fclose(file), -1, NULL);
if (file) SYSTEM_ERROR_CHECK(fclose(file), -1, nullptr);
return EXIT_SUCCESS;
}
cmake_minimum_required(VERSION 2.6)
project(f_mfcc)
find_package(blc_channel REQUIRED)
find_package(blc_program REQUIRED)
find_path(FFTW_INCLUDE_DIR fftw3.h REQUIRED)
find_library(FFTW_LIBRARIES NAMES fftw3f REQUIRED)
add_definitions(${BL_DEFINITIONS})
include_directories(${BL_INCLUDE_DIRS} libmfcc ${FFTW_INCLUDE_DIR})
add_executable(f_mfcc f_mfcc.cpp libmfcc/libmfcc.cpp)
target_link_libraries(f_mfcc ${BL_LIBRARIES} ${FFTW_LIBRARIES})
#include "blc_channel.h"
#include "blc_program.h"
#include <fftw3.h>
#include <unistd.h>
#include <math.h> //sqrtf
#include <pthread.h>
#include "libmfcc.h"
#define DEFAULT_OUTPUT_NAME "mfcc<pid>"
static const int columns_nb=64, rows_nb=16;
static const int mfcc_coeffs_nb=26;
static const int sample_rate=44100;
static const int filters_nb=48;
int main(int argc, char**argv){
char const *display;
blc_channel input, output;
char const *input_name, *output_name, *period_str;
int fftw_length, period;
size_t i;
fftwf_complex *intermediate;
fftwf_plan fftw_plan;
double *fftw_array;
blc_program_add_option(&display, 'd', "display", NULL, "Display the result as text graph", NULL);
blc_program_add_option(&output_name, 'o', "output", "blc_channel", "channel containing the mfcc", DEFAULT_OUTPUT_NAME);
blc_program_add_option(&period_str, 'p', "period", "integer", "Period in ms (0 as fast as possible).", "0");
blc_program_add_parameter(&input_name, "input blc_channel", 1, "Channel containing the sound data", NULL);
blc_program_init(&argc, &argv, blc_quit);
SSCANF(1, period_str, "%d", &period);
//Open input channel and check the type of values
input.open(input_name, BLC_CHANNEL_READ);
if (input.type!='FL32') EXIT_ON_CHANNEL_ERROR(&input,"Input must be of type 'FL32'");
fftw_length=input.total_length/2+1;
if (strcmp(output_name, DEFAULT_OUTPUT_NAME)==0) asprintf((char**)&output_name, ":mfcc%d", getpid()); //This will not be freed until the end (never mind)
//Create a channel or reopen a previous one.
output.create_or_open(output_name, BLC_CHANNEL_WRITE, 'FL32', 'NDEF', 1, mfcc_coeffs_nb);
output.publish();
blc_loop_try_add_waiting_semaphore(output.sem_ack_data);
blc_loop_try_add_posting_semaphore(output.sem_new_data);
blc_loop_try_add_waiting_semaphore(input.sem_new_data);
blc_loop_try_add_posting_semaphore(input.sem_ack_data);
fftw_array=MANY_ALLOCATIONS(fftw_length, double);
intermediate=fftwf_alloc_complex(fftw_length);
fftw_plan = fftwf_plan_dft_r2c_1d(input.total_length, input.floats, intermediate , FFTW_ESTIMATE);
BLC_COMMAND_LOOP(period*1000){
fftwf_execute(fftw_plan); /* repeat as needed */
FOR(i, fftw_length){
fftw_array[i]=sqrt(intermediate[i][0]*intermediate[i][0]+intermediate[i][1]*intermediate[i][1]);
}
FOR(i, mfcc_coeffs_nb) output.floats[i] = GetCoefficient(fftw_array, sample_rate, filters_nb, fftw_length, i);
if(display){
blc_fprint_float_graph(stderr, output.floats, output.total_length, "mfcc", columns_nb, rows_nb, 20, 0, "Coeff", "Intensity");
blc_eprint_cursor_up(rows_nb);
}
}
fftwf_destroy_plan(fftw_plan);
blc_eprint_cursor_down(rows_nb);
return EXIT_SUCCESS;
}
The MIT License
Copyright (c) 2010 Jeremy Sawruk
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
\ No newline at end of file
libMFCC - C library for computing Mel-Frequency Cepstral Coefficients
Version 0.1 - 5 October 2010
Implementation based on:
"Musical Instrument Identification using Multiscale Mel-Frequency Cepstral Coefficients" by Bob L. Sturm et al., 2010
Currently available at: http://imi.aau.dk/~bst/publications/Sturm2010b.pdf
Files in this distribution:
LICENSE => License to use and redistribute libMFCC (MIT License)
README => This file
TODO => Things not implemented in the current library
libmfcc.c => libMFCC algorithm implementation
libmfcc.h => libMFCC header file
Sturm2010b.pdf => Sturm et al. reference document
/libmfcc_example => Microsoft Visual Studio 2010 project with example code demonstrating use of libMFCC
\ No newline at end of file
TODO
- Add windowing functions
- Create internal FFT implementation?
- Add other MFCC implementations? (Is it possible to make an open-source implementation of the algorithm that Slaney uses?)
- Make into library (e.g. DLL) for use in other languages
- Performance optimizations
\ No newline at end of file
/*
* libmfcc.c - Code implementation for libMFCC
* Copyright (c) 2010 Jeremy Sawruk
*
* This code is released under the MIT License.
* For conditions of distribution and use, see the license in LICENSE
*/
#include <math.h>
#include "libmfcc.h"
/*
* Computes the specified (mth) MFCC
*
* spectralData - array of doubles containing the results of FFT computation. This data is already assumed to be purely real
* samplingRate - the rate that the original time-series data was sampled at (i.e 44100)
* NumFilters - the number of filters to use in the computation. Recommended value = 48
* binSize - the size of the spectralData array, usually a power of 2
* m - The mth MFCC coefficient to compute
*
*/
double GetCoefficient(double* spectralData, unsigned int samplingRate, unsigned int NumFilters, unsigned int binSize, unsigned int m)
{
double result = 0.0f;
double outerSum = 0.0f;
double innerSum = 0.0f;
unsigned int k, l;
// 0 <= m < L
if(m >= NumFilters)
{
// This represents an error condition - the specified coefficient is greater than or equal to the number of filters. The behavior in this case is undefined.
return 0.0f;
}
result = NormalizationFactor(NumFilters, m);
for(l = 1; l <= NumFilters; l++)
{
// Compute inner sum
innerSum = 0.0f;
for(k = 0; k < binSize - 1; k++)
{
innerSum += fabs(spectralData[k] * GetFilterParameter(samplingRate, binSize, k, l));
}
if(innerSum > 0.0f)
{
innerSum = log(innerSum); // The log of 0 is undefined, so don't use it
}
innerSum = innerSum * cos(((m * PI) / NumFilters) * (l - 0.5f));
outerSum += innerSum;
}
result *= outerSum;
return result;
}
/*
* Computes the Normalization Factor (Equation 6)
* Used for internal computation only - not to be called directly
*/
double NormalizationFactor(int NumFilters, int m)
{
double normalizationFactor = 0.0f;
if(m == 0)
{
normalizationFactor = sqrt(1.0f / NumFilters);
}
else
{
normalizationFactor = sqrt(2.0f / NumFilters);
}
return normalizationFactor;
}
/*
* Compute the filter parameter for the specified frequency and filter bands (Eq. 2)
* Used for internal computation only - not the be called directly
*/
double GetFilterParameter(unsigned int samplingRate, unsigned int binSize, unsigned int frequencyBand, unsigned int filterBand)
{
double filterParameter = 0.0f;
double boundary = (frequencyBand * samplingRate) / binSize; // k * Fs / N
double prevCenterFrequency = GetCenterFrequency(filterBand - 1); // fc(l - 1) etc.
double thisCenterFrequency = GetCenterFrequency(filterBand);
double nextCenterFrequency = GetCenterFrequency(filterBand + 1);
if(boundary >= 0 && boundary < prevCenterFrequency)
{
filterParameter = 0.0f;
}
else if(boundary >= prevCenterFrequency && boundary < thisCenterFrequency)
{
filterParameter = (boundary - prevCenterFrequency) / (thisCenterFrequency - prevCenterFrequency);
filterParameter *= GetMagnitudeFactor(filterBand);
}
else if(boundary >= thisCenterFrequency && boundary < nextCenterFrequency)
{
filterParameter = (boundary - nextCenterFrequency) / (thisCenterFrequency - nextCenterFrequency);
filterParameter *= GetMagnitudeFactor(filterBand);
}
else if(boundary >= nextCenterFrequency && boundary < samplingRate)
{
filterParameter = 0.0f;
}
return filterParameter;
}
/*
* Compute the band-dependent magnitude factor for the given filter band (Eq. 3)
* Used for internal computation only - not the be called directly
*/
double GetMagnitudeFactor(unsigned int filterBand)
{
double magnitudeFactor = 0.0f;
if(filterBand >= 1 && filterBand <= 14)
{
magnitudeFactor = 0.015;
}
else if(filterBand >= 15 && filterBand <= 48)
{
magnitudeFactor = 2.0f / (GetCenterFrequency(filterBand + 1) - GetCenterFrequency(filterBand -1));
}
return magnitudeFactor;
}
/*
* Compute the center frequency (fc) of the specified filter band (l) (Eq. 4)
* This where the mel-frequency scaling occurs. Filters are specified so that their
* center frequencies are equally spaced on the mel scale
* Used for internal computation only - not the be called directly
*/
double GetCenterFrequency(unsigned int filterBand)
{
double centerFrequency = 0.0f;
double exponent;
if(filterBand == 0)
{
centerFrequency = 0;
}
else if(filterBand >= 1 && filterBand <= 14)
{
centerFrequency = (200.0f * filterBand) / 3.0f;
}
else
{
exponent = filterBand - 14.0f;
centerFrequency = pow(1.0711703, exponent);
centerFrequency *= 1073.4;
}
return centerFrequency;
}
\ No newline at end of file
/*
* libmfcc.h - Header for libMFCC
* Copyright (c) 2010 Jeremy Sawruk
*
* This code is released under the MIT License.
* For conditions of distribution and use, see the license in LICENSE
*/
#pragma once
#define PI 3.14159265358979323846264338327
// Returns the specified (mth) MFCC
double GetCoefficient(double* spectralData, unsigned int samplingRate, unsigned int NumFilters, unsigned int binSize, unsigned int m);
// Compute the normalization factor (For internal computation only - not to be called directly)
double NormalizationFactor(int NumFilters, int m);
// Compute the filter parameter for the specified frequency and filter bands (For internal computation only - not the be called directly)
double GetFilterParameter(unsigned int samplingRate, unsigned int binSize, unsigned int frequencyBand, unsigned int filterBand);
// Compute the band-dependent magnitude factor for the given filter band (For internal computation only - not the be called directly)
double GetMagnitudeFactor(unsigned int filterBand);
// Compute the center frequency (fc) of the specified filter band (l) (For internal computation only - not the be called directly)
double GetCenterFrequency(unsigned int filterBand);
\ No newline at end of file
/*
* example.c - example of using libMFCC
* Written 2010 - J. Sawruk
*/
#include <stdio.h>
#include <string.h> // For memset
#include "libmfcc.h"
int main(void)
{
// Read in sample data from sample.dat
// sample.dat contains an 8192-point spectrum from a sine wave at 440Hz (A) in double precision
// Spectrum was computed using FFTW (http://www.fftw.org/)
// Data was not windowed (rectangular)
// Holds the spectrum data to be analyzed
double spectrum[8192];
// Pointer to the sample data file
FILE *sampleFile;
// Index counter - used to keep track of which data point is being read in
int i = 0;
// Determine which MFCC coefficient to compute
unsigned int coeff;
// Holds the value of the computed coefficient
double mfcc_result;
// Initialize the spectrum
memset(&spectrum, 0, sizeof(spectrum));
// Open the sample spectrum data
sampleFile = fopen("sample.dat","rb");
// Read in the contents of the sample file
while(fscanf(sampleFile, "%lf", &spectrum[i]) != EOF) // %lf tells fscanf to read a double
{
i++;
}
// Close the sample file
fclose(sampleFile);
// Compute the first 13 coefficients
for(coeff = 0; coeff < 13; coeff++)
{
mfcc_result = GetCoefficient(spectrum, 44100, 48, 128, coeff);
printf("%i %f\n", coeff, mfcc_result);
}
getchar();
return 0;
}
\ No newline at end of file
This diff is collapsed.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment