From 2ff0dfe30c74be864e6331bd4e1b30fd2d32f06a Mon Sep 17 00:00:00 2001 From: "Long N. Le" Date: Fri, 3 Mar 2017 18:12:48 -0600 Subject: [PATCH 1/8] added event_detector/ --- .../apps/audio_module/event_detector/Makefile | 18 ++++ .../apps/audio_module/event_detector/main.c | 88 +++++++++++++++++++ 2 files changed, 106 insertions(+) create mode 100644 software/apps/audio_module/event_detector/Makefile create mode 100644 software/apps/audio_module/event_detector/main.c diff --git a/software/apps/audio_module/event_detector/Makefile b/software/apps/audio_module/event_detector/Makefile new file mode 100644 index 00000000..0e331e06 --- /dev/null +++ b/software/apps/audio_module/event_detector/Makefile @@ -0,0 +1,18 @@ +# makefile for user application + +# the current directory +APP_DIR := $(dir $(abspath $(lastword $(MAKEFILE_LIST)))) + +LIB_DIR := $(HOME)/gcwa/c + +# files needed for this code +C_SRCS := $(wildcard *.c) $(LIB_DIR)/ridgeTracker.c $(LIB_DIR)/dynArray.c +INCLUDE_PATHS += . $(LIB_DIR) $(LIB_DIR)/kiss_fft $(LIB_DIR)/kiss_fft/tools + +CFLAGS += -DFIXED_POINT=16 -std=c99 + +TOCK_BOARD = audio_module + +# include makefile settings that are shared between applications +include ../../AppMakefile.mk + diff --git a/software/apps/audio_module/event_detector/main.c b/software/apps/audio_module/event_detector/main.c new file mode 100644 index 00000000..f3e25d83 --- /dev/null +++ b/software/apps/audio_module/event_detector/main.c @@ -0,0 +1,88 @@ + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "app_watchdog.h" +#include "signpost_api.h" +#include "simple_post.h" +#include "_kiss_fft_guts.h" +#include "kiss_fftr.h" +#include "common.h" +#include "ridgeTracker.h" + +static int freqAnalyze(kiss_fft_scalar *frame, kiss_fft_scalar *spec); + +int main (void) { + int err = SUCCESS; + printf("[Audio Module] Event Detector\n"); + + // initialize ADC + err = adc_initialize(); + if (err < SUCCESS) { + printf("Initialize errored: %d\n", err); + } + + static kiss_fft_scalar frame[BUF_LEN]; + static kiss_fft_scalar spec[FRE_LEN]; + static kiss_fft_scalar snrOut[FRE_LEN]; + memset(frame, 0, BUF_LEN*sizeof(kiss_fft_scalar)); + size_t frameIdx = 0; + size_t globTI = 0; + ridgeTracker_init(); + while (true) { + // read data from ADC + err = adc_read_single_sample(3); + if (err < SUCCESS) { + printf("ADC read error: %d\n", err); + } + uint16_t sample = err & 0xFFFF; + + frameIdx = (frameIdx+INC_LEN) % BUF_LEN; + // perform FFT + if (freqAnalyze(frame, spec)){ + printf("freqAnalyze() failed!\n"); + return 1; + } + + // ridgetracker update + ridgeTracker_update(spec, snrOut); + } +} + +static int freqAnalyze(kiss_fft_scalar *frame, kiss_fft_scalar *spec){ + // FFT data structure + kiss_fftr_cfg cfg; + kiss_fft_scalar in[BUF_LEN]; + kiss_fft_cpx out[FRE_LEN + 1]; + + static kiss_fft_scalar win[BUF_LEN] = {0x0000,0x0001,0x0004,0x000b,0x0013,0x001e,0x002c,0x003c,0x004f,0x0064,0x007b,0x0095,0x00b2,0x00d0,0x00f2,0x0115,0x013c,0x0164,0x018f,0x01bd,0x01ec,0x021f,0x0253,0x028a,0x02c4,0x0300,0x033e,0x037e,0x03c1,0x0406,0x044e,0x0497,0x04e3,0x0532,0x0583,0x05d5,0x062b,0x0682,0x06dc,0x0737,0x0796,0x07f6,0x0858,0x08bd,0x0923,0x098c,0x09f7,0x0a64,0x0ad3,0x0b44,0x0bb8,0x0c2d,0x0ca4,0x0d1d,0x0d98,0x0e15,0x0e94,0x0f15,0x0f98,0x101d,0x10a4,0x112c,0x11b6,0x1242,0x12d0,0x1360,0x13f1,0x1484,0x1519,0x15af,0x1647,0x16e0,0x177c,0x1818,0x18b7,0x1956,0x19f8,0x1a9a,0x1b3f,0x1be4,0x1c8b,0x1d34,0x1ddd,0x1e88,0x1f35,0x1fe2,0x2091,0x2141,0x21f3,0x22a5,0x2359,0x240d,0x24c3,0x257a,0x2632,0x26eb,0x27a5,0x285f,0x291b,0x29d8,0x2a95,0x2b53,0x2c12,0x2cd2,0x2d93,0x2e54,0x2f16,0x2fd9,0x309c,0x3160,0x3224,0x32e9,0x33ae,0x3474,0x353b,0x3601,0x36c9,0x3790,0x3858,0x3920,0x39e9,0x3ab1,0x3b7a,0x3c43,0x3d0c,0x3dd6,0x3e9f,0x3f68,0x4032,0x40fb,0x41c5,0x428e,0x4357,0x4420,0x44e9,0x45b2,0x467b,0x4743,0x480b,0x48d3,0x499a,0x4a61,0x4b28,0x4bee,0x4cb3,0x4d79,0x4e3d,0x4f01,0x4fc5,0x5088,0x514a,0x520c,0x52cd,0x538d,0x544c,0x550b,0x55c9,0x5686,0x5742,0x57fd,0x58b8,0x5971,0x5a29,0x5ae1,0x5b97,0x5c4c,0x5d00,0x5db3,0x5e65,0x5f16,0x5fc5,0x6074,0x6121,0x61cc,0x6277,0x6320,0x63c7,0x646e,0x6513,0x65b6,0x6658,0x66f9,0x6798,0x6835,0x68d1,0x696c,0x6a04,0x6a9c,0x6b31,0x6bc5,0x6c57,0x6ce7,0x6d76,0x6e03,0x6e8e,0x6f17,0x6f9f,0x7025,0x70a8,0x712a,0x71aa,0x7228,0x72a5,0x731f,0x7397,0x740d,0x7481,0x74f3,0x7564,0x75d2,0x763e,0x76a7,0x770f,0x7775,0x77d8,0x783a,0x7899,0x78f6,0x7950,0x79a9,0x79ff,0x7a53,0x7aa5,0x7af5,0x7b42,0x7b8d,0x7bd5,0x7c1c,0x7c60,0x7ca1,0x7ce1,0x7d1e,0x7d58,0x7d90,0x7dc6,0x7dfa,0x7e2b,0x7e59,0x7e86,0x7eaf,0x7ed7,0x7efc,0x7f1e,0x7f3e,0x7f5c,0x7f77,0x7f90,0x7fa6,0x7fba,0x7fcb,0x7fda,0x7fe6,0x7ff0,0x7ff8,0x7ffd,0x7fff,0x7fff,0x7ffd,0x7ff8,0x7ff0,0x7fe6,0x7fda,0x7fcb,0x7fba,0x7fa6,0x7f90,0x7f77,0x7f5c,0x7f3e,0x7f1e,0x7efc,0x7ed7,0x7eaf,0x7e86,0x7e59,0x7e2b,0x7dfa,0x7dc6,0x7d90,0x7d58,0x7d1e,0x7ce1,0x7ca1,0x7c60,0x7c1c,0x7bd5,0x7b8d,0x7b42,0x7af5,0x7aa5,0x7a53,0x79ff,0x79a9,0x7950,0x78f6,0x7899,0x783a,0x77d8,0x7775,0x770f,0x76a7,0x763e,0x75d2,0x7564,0x74f3,0x7481,0x740d,0x7397,0x731f,0x72a5,0x7228,0x71aa,0x712a,0x70a8,0x7025,0x6f9f,0x6f17,0x6e8e,0x6e03,0x6d76,0x6ce7,0x6c57,0x6bc5,0x6b31,0x6a9c,0x6a04,0x696c,0x68d1,0x6835,0x6798,0x66f9,0x6658,0x65b6,0x6513,0x646e,0x63c7,0x6320,0x6277,0x61cc,0x6121,0x6074,0x5fc5,0x5f16,0x5e65,0x5db3,0x5d00,0x5c4c,0x5b97,0x5ae1,0x5a29,0x5971,0x58b8,0x57fd,0x5742,0x5686,0x55c9,0x550b,0x544c,0x538d,0x52cd,0x520c,0x514a,0x5088,0x4fc5,0x4f01,0x4e3d,0x4d79,0x4cb3,0x4bee,0x4b28,0x4a61,0x499a,0x48d3,0x480b,0x4743,0x467b,0x45b2,0x44e9,0x4420,0x4357,0x428e,0x41c5,0x40fb, + 0x4032,0x3f68,0x3e9f,0x3dd6,0x3d0c,0x3c43,0x3b7a,0x3ab1,0x39e9,0x3920,0x3858,0x3790,0x36c9,0x3601,0x353b,0x3474,0x33ae,0x32e9,0x3224,0x3160,0x309c,0x2fd9,0x2f16,0x2e54,0x2d93,0x2cd2,0x2c12,0x2b53,0x2a95,0x29d8,0x291b,0x285f,0x27a5,0x26eb,0x2632,0x257a,0x24c3,0x240d,0x2359,0x22a5,0x21f3,0x2141,0x2091,0x1fe2,0x1f35,0x1e88,0x1ddd,0x1d34,0x1c8b,0x1be4,0x1b3f,0x1a9a,0x19f8,0x1956,0x18b7,0x1818,0x177c,0x16e0,0x1647,0x15af,0x1519,0x1484,0x13f1,0x1360,0x12d0,0x1242,0x11b6,0x112c,0x10a4,0x101d,0x0f98,0x0f15,0x0e94,0x0e15,0x0d98,0x0d1d,0x0ca4,0x0c2d,0x0bb8,0x0b44,0x0ad3,0x0a64,0x09f7,0x098c,0x0923,0x08bd,0x0858,0x07f6,0x0796,0x0737,0x06dc,0x0682,0x062b,0x05d5,0x0583,0x0532,0x04e3,0x0497,0x044e,0x0406,0x03c1,0x037e,0x033e,0x0300,0x02c4,0x028a,0x0253,0x021f,0x01ec,0x01bd,0x018f,0x0164,0x013c,0x0115,0x00f2,0x00d0,0x00b2,0x0095,0x007b,0x0064,0x004f,0x003c,0x002c,0x001e,0x0013,0x000b,0x0004,0x0001,0x0000}; + // windowing and prescaling + for (size_t i = 0; i < BUF_LEN; i++){ + in[i] = S_MUL(frame[i],win[i]) << 4; + } + + if ((cfg = kiss_fftr_alloc(BUF_LEN, 0/*is_inverse_fft*/, NULL, NULL)) == NULL){ + printf("Not enough memory?\n"); + return 1; + }else{ + kiss_fftr(cfg, in, out); + free(cfg); + } + + for (size_t k = 0; k < FRE_LEN; k++){ + spec[k] = MAG(out[k].r,out[k].i); + } + + return 0; +} From bbac57a54a3a7b928040bd0b26dbcb6627c15d55 Mon Sep 17 00:00:00 2001 From: "Long N. Le" Date: Tue, 7 Mar 2017 10:59:03 -0600 Subject: [PATCH 2/8] compiled event detector --- .../apps/audio_module/event_detector/Makefile | 7 +- .../audio_module/event_detector/README.md | 14 + .../apps/audio_module/event_detector/common.h | 21 + .../audio_module/event_detector/dynArray.c | 64 +++ .../audio_module/event_detector/dynArray.h | 21 + .../audio_module/event_detector/kiss_fft.c | 408 ++++++++++++++++++ .../audio_module/event_detector/kiss_fftr.c | 159 +++++++ .../audio_module/event_detector/log2fix.c | 69 +++ .../event_detector/ridgeTracker.c | 184 ++++++++ .../event_detector/ridgeTracker.h | 21 + 10 files changed, 964 insertions(+), 4 deletions(-) create mode 100644 software/apps/audio_module/event_detector/README.md create mode 100644 software/apps/audio_module/event_detector/common.h create mode 100644 software/apps/audio_module/event_detector/dynArray.c create mode 100644 software/apps/audio_module/event_detector/dynArray.h create mode 100644 software/apps/audio_module/event_detector/kiss_fft.c create mode 100644 software/apps/audio_module/event_detector/kiss_fftr.c create mode 100644 software/apps/audio_module/event_detector/log2fix.c create mode 100644 software/apps/audio_module/event_detector/ridgeTracker.c create mode 100644 software/apps/audio_module/event_detector/ridgeTracker.h diff --git a/software/apps/audio_module/event_detector/Makefile b/software/apps/audio_module/event_detector/Makefile index 0e331e06..e32011b8 100644 --- a/software/apps/audio_module/event_detector/Makefile +++ b/software/apps/audio_module/event_detector/Makefile @@ -3,13 +3,12 @@ # the current directory APP_DIR := $(dir $(abspath $(lastword $(MAKEFILE_LIST)))) -LIB_DIR := $(HOME)/gcwa/c - # files needed for this code -C_SRCS := $(wildcard *.c) $(LIB_DIR)/ridgeTracker.c $(LIB_DIR)/dynArray.c -INCLUDE_PATHS += . $(LIB_DIR) $(LIB_DIR)/kiss_fft $(LIB_DIR)/kiss_fft/tools +C_SRCS := $(wildcard *.c) +INCLUDE_PATHS += . ./kiss_fft ./kiss_fft/tools ./log2fix/ CFLAGS += -DFIXED_POINT=16 -std=c99 +CFLAGS += -Wno-type-limits -Wno-sign-compare TOCK_BOARD = audio_module diff --git a/software/apps/audio_module/event_detector/README.md b/software/apps/audio_module/event_detector/README.md new file mode 100644 index 00000000..1041d620 --- /dev/null +++ b/software/apps/audio_module/event_detector/README.md @@ -0,0 +1,14 @@ +# Fixed-point ridge tracker in C + +## Quick start +``` +make clean +make +./main input.wav +``` +Use test.ipynb to visualize the resulting \*.mat files + +## Prerequisites +* Fixed-point FFT: git clone https://github.com/longle2718/kiss_fft +* Fixed-point log: git clone https://github.com/dmoulding/log2fix +* libsndfile: sudo apt-get install libsndfile1-dev diff --git a/software/apps/audio_module/event_detector/common.h b/software/apps/audio_module/event_detector/common.h new file mode 100644 index 00000000..8a15780d --- /dev/null +++ b/software/apps/audio_module/event_detector/common.h @@ -0,0 +1,21 @@ +#ifndef __COMMON_H__ +#define __COMMON_H__ + +#include "_kiss_fft_guts.h" + +#define BUF_LEN (512) +#define INC_LEN (BUF_LEN/2) +#define FRE_LEN (BUF_LEN/2) + +#define S_SQ(a) S_MUL((a),(a)) // scalar square +#define S_DIV(a,b) S_MUL((a),SAMP_MAX/(b)) // a < b +#define MAX(a,b) (((a) > (b)) ? (a) : (b)) +#define MIN(a,b) (((a) < (b)) ? (a) : (b)) +#define ABS(a) (((a) < 0) ? -(a): (a)) +// Magnitude Estimator +// http://dspguru.com/book/export/html/62 +#define ALPHA ((kiss_fft_scalar)(0.947543636291*SAMP_MAX)) +#define BETA ((kiss_fft_scalar)(0.392485425092*SAMP_MAX)) +#define MAG(r,i) ( S_MUL(ALPHA,MAX(ABS((r)),ABS((i)))) + S_MUL(BETA,MIN(ABS((r)),ABS((i)))) ) + +#endif diff --git a/software/apps/audio_module/event_detector/dynArray.c b/software/apps/audio_module/event_detector/dynArray.c new file mode 100644 index 00000000..07d82b65 --- /dev/null +++ b/software/apps/audio_module/event_detector/dynArray.c @@ -0,0 +1,64 @@ +#include +#include +#include "dynArray.h" +#include "_kiss_fft_guts.h" + +void initArray(Array *a, size_t initialSize) { + a->maxSize = initialSize; + a->used = 0; + + a->size = a->maxSize; + printf("Init array size = %zu\n",a->size); + a->SNR = malloc(a->size * sizeof(kiss_fft_scalar)); + a->FI = malloc(a->size * sizeof(size_t)); + a->TI = malloc(a->size * sizeof(size_t)); +} + +void insertArray(Array *a, kiss_fft_scalar snr, size_t fi, size_t ti) { + // a->used is the number of used entries, because a->array[a->used++] updates a->used only *after* the array has been accessed. + // Therefore a->used can go up to a->size + if (a->size > 0 && a->used == a->size) { + a->maxSize *= 2; + printf("Doubled array max size to %zu\n",a->maxSize); + } + if (a->size < a->maxSize) { + a->size = a->maxSize; + // with NULL ptr, realloc == malloc + a->SNR = realloc(a->SNR, a->size * sizeof(kiss_fft_scalar)); + a->FI = realloc(a->FI, a->size * sizeof(size_t)); + a->TI = realloc(a->TI, a->size * sizeof(size_t)); + } + a->SNR[a->used] = snr; + a->FI[a->used] = fi; + a->TI[a->used] = ti; + a->used++; +} + +void freeArray(Array *a) { + free(a->SNR); + free(a->FI); + free(a->TI); + a->SNR = NULL; + a->FI = NULL; + a->TI = NULL; + a->used = a->size = 0; +} +size_t getMaxDurArray(Array a){ + if (a.used == 0) + return 0; + + size_t minTI = a.TI[0]; + size_t maxTI = a.TI[a.used-1]; + return maxTI-minTI; +} + +kiss_fft_scalar getAvgSNRArray(Array a){ + if (a.used == 0) + return 0; + + SAMPPROD sum = 0; + for (size_t k=0; k < a.used; k++){ + sum += a.SNR[k]; + } + return (kiss_fft_scalar)(sum/a.used); +} diff --git a/software/apps/audio_module/event_detector/dynArray.h b/software/apps/audio_module/event_detector/dynArray.h new file mode 100644 index 00000000..47918f37 --- /dev/null +++ b/software/apps/audio_module/event_detector/dynArray.h @@ -0,0 +1,21 @@ +#ifndef __DYNARRAY_H__ +#define __DYNARRAY_H__ + +#include "kiss_fft.h" + +typedef struct { + kiss_fft_scalar *SNR; // log SNR + size_t *FI; // frequency index + size_t *TI; // time index + size_t size; + size_t used; + size_t maxSize; +} Array; + +void initArray(Array *a, size_t initialSize); +void insertArray(Array *a, kiss_fft_scalar snr, size_t fi, size_t ti); +void freeArray(Array *a); + +size_t getMaxDurArray(Array a); +kiss_fft_scalar getAvgSNRArray(Array a); +#endif diff --git a/software/apps/audio_module/event_detector/kiss_fft.c b/software/apps/audio_module/event_detector/kiss_fft.c new file mode 100644 index 00000000..465d6c97 --- /dev/null +++ b/software/apps/audio_module/event_detector/kiss_fft.c @@ -0,0 +1,408 @@ +/* +Copyright (c) 2003-2010, Mark Borgerding + +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. + * Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + + +#include "_kiss_fft_guts.h" +/* The guts header contains all the multiplication and addition macros that are defined for + fixed or floating point complex numbers. It also delares the kf_ internal functions. + */ + +static void kf_bfly2( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + int m + ) +{ + kiss_fft_cpx * Fout2; + kiss_fft_cpx * tw1 = st->twiddles; + kiss_fft_cpx t; + Fout2 = Fout + m; + do{ + C_FIXDIV(*Fout,2); C_FIXDIV(*Fout2,2); + + C_MUL (t, *Fout2 , *tw1); + tw1 += fstride; + C_SUB( *Fout2 , *Fout , t ); + C_ADDTO( *Fout , t ); + ++Fout2; + ++Fout; + }while (--m); +} + +static void kf_bfly4( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + const size_t m + ) +{ + kiss_fft_cpx *tw1,*tw2,*tw3; + kiss_fft_cpx scratch[6]; + size_t k=m; + const size_t m2=2*m; + const size_t m3=3*m; + + + tw3 = tw2 = tw1 = st->twiddles; + + do { + C_FIXDIV(*Fout,4); C_FIXDIV(Fout[m],4); C_FIXDIV(Fout[m2],4); C_FIXDIV(Fout[m3],4); + + C_MUL(scratch[0],Fout[m] , *tw1 ); + C_MUL(scratch[1],Fout[m2] , *tw2 ); + C_MUL(scratch[2],Fout[m3] , *tw3 ); + + C_SUB( scratch[5] , *Fout, scratch[1] ); + C_ADDTO(*Fout, scratch[1]); + C_ADD( scratch[3] , scratch[0] , scratch[2] ); + C_SUB( scratch[4] , scratch[0] , scratch[2] ); + C_SUB( Fout[m2], *Fout, scratch[3] ); + tw1 += fstride; + tw2 += fstride*2; + tw3 += fstride*3; + C_ADDTO( *Fout , scratch[3] ); + + if(st->inverse) { + Fout[m].r = scratch[5].r - scratch[4].i; + Fout[m].i = scratch[5].i + scratch[4].r; + Fout[m3].r = scratch[5].r + scratch[4].i; + Fout[m3].i = scratch[5].i - scratch[4].r; + }else{ + Fout[m].r = scratch[5].r + scratch[4].i; + Fout[m].i = scratch[5].i - scratch[4].r; + Fout[m3].r = scratch[5].r - scratch[4].i; + Fout[m3].i = scratch[5].i + scratch[4].r; + } + ++Fout; + }while(--k); +} + +static void kf_bfly3( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + size_t m + ) +{ + size_t k=m; + const size_t m2 = 2*m; + kiss_fft_cpx *tw1,*tw2; + kiss_fft_cpx scratch[5]; + kiss_fft_cpx epi3; + epi3 = st->twiddles[fstride*m]; + + tw1=tw2=st->twiddles; + + do{ + C_FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3); + + C_MUL(scratch[1],Fout[m] , *tw1); + C_MUL(scratch[2],Fout[m2] , *tw2); + + C_ADD(scratch[3],scratch[1],scratch[2]); + C_SUB(scratch[0],scratch[1],scratch[2]); + tw1 += fstride; + tw2 += fstride*2; + + Fout[m].r = Fout->r - HALF_OF(scratch[3].r); + Fout[m].i = Fout->i - HALF_OF(scratch[3].i); + + C_MULBYSCALAR( scratch[0] , epi3.i ); + + C_ADDTO(*Fout,scratch[3]); + + Fout[m2].r = Fout[m].r + scratch[0].i; + Fout[m2].i = Fout[m].i - scratch[0].r; + + Fout[m].r -= scratch[0].i; + Fout[m].i += scratch[0].r; + + ++Fout; + }while(--k); +} + +static void kf_bfly5( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + int m + ) +{ + kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4; + int u; + kiss_fft_cpx scratch[13]; + kiss_fft_cpx * twiddles = st->twiddles; + kiss_fft_cpx *tw; + kiss_fft_cpx ya,yb; + ya = twiddles[fstride*m]; + yb = twiddles[fstride*2*m]; + + Fout0=Fout; + Fout1=Fout0+m; + Fout2=Fout0+2*m; + Fout3=Fout0+3*m; + Fout4=Fout0+4*m; + + tw=st->twiddles; + for ( u=0; ur += scratch[7].r + scratch[8].r; + Fout0->i += scratch[7].i + scratch[8].i; + + scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r); + scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r); + + scratch[6].r = S_MUL(scratch[10].i,ya.i) + S_MUL(scratch[9].i,yb.i); + scratch[6].i = -S_MUL(scratch[10].r,ya.i) - S_MUL(scratch[9].r,yb.i); + + C_SUB(*Fout1,scratch[5],scratch[6]); + C_ADD(*Fout4,scratch[5],scratch[6]); + + scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r); + scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r); + scratch[12].r = - S_MUL(scratch[10].i,yb.i) + S_MUL(scratch[9].i,ya.i); + scratch[12].i = S_MUL(scratch[10].r,yb.i) - S_MUL(scratch[9].r,ya.i); + + C_ADD(*Fout2,scratch[11],scratch[12]); + C_SUB(*Fout3,scratch[11],scratch[12]); + + ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4; + } +} + +/* perform the butterfly for one stage of a mixed radix FFT */ +static void kf_bfly_generic( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + int m, + int p + ) +{ + int u,k,q1,q; + kiss_fft_cpx * twiddles = st->twiddles; + kiss_fft_cpx t; + int Norig = st->nfft; + + kiss_fft_cpx * scratch = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC(sizeof(kiss_fft_cpx)*p); + + for ( u=0; u=Norig) twidx-=Norig; + C_MUL(t,scratch[q] , twiddles[twidx] ); + C_ADDTO( Fout[ k ] ,t); + } + k += m; + } + } + KISS_FFT_TMP_FREE(scratch); +} + +static +void kf_work( + kiss_fft_cpx * Fout, + const kiss_fft_cpx * f, + const size_t fstride, + int in_stride, + int * factors, + const kiss_fft_cfg st + ) +{ + kiss_fft_cpx * Fout_beg=Fout; + const int p=*factors++; /* the radix */ + const int m=*factors++; /* stage's fft length/p */ + const kiss_fft_cpx * Fout_end = Fout + p*m; + +#ifdef _OPENMP + // use openmp extensions at the + // top-level (not recursive) + if (fstride==1 && p<=5) + { + int k; + + // execute the p different work units in different threads +# pragma omp parallel for + for (k=0;k floor_sqrt) + p = n; /* no more factors, skip to end */ + } + n /= p; + *facbuf++ = p; + *facbuf++ = n; + } while (n > 1); +} + +/* + * + * User-callable function to allocate all necessary storage space for the fft. + * + * The return value is a contiguous block of memory, allocated with malloc. As such, + * It can be freed with free(), rather than a kiss_fft-specific function. + * */ +kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem ) +{ + kiss_fft_cfg st=NULL; + size_t memneeded = sizeof(struct kiss_fft_state) + + sizeof(kiss_fft_cpx)*(nfft-1); /* twiddle factors*/ + + if ( lenmem==NULL ) { + st = ( kiss_fft_cfg)KISS_FFT_MALLOC( memneeded ); + }else{ + if (mem != NULL && *lenmem >= memneeded) + st = (kiss_fft_cfg)mem; + *lenmem = memneeded; + } + if (st) { + int i; + st->nfft=nfft; + st->inverse = inverse_fft; + + for (i=0;iinverse) + phase *= -1; + kf_cexp(st->twiddles+i, phase ); + } + + kf_factor(nfft,st->factors); + } + return st; +} + + +void kiss_fft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride) +{ + if (fin == fout) { + //NOTE: this is not really an in-place FFT algorithm. + //It just performs an out-of-place FFT into a temp buffer + kiss_fft_cpx * tmpbuf = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC( sizeof(kiss_fft_cpx)*st->nfft); + kf_work(tmpbuf,fin,1,in_stride, st->factors,st); + memcpy(fout,tmpbuf,sizeof(kiss_fft_cpx)*st->nfft); + KISS_FFT_TMP_FREE(tmpbuf); + }else{ + kf_work( fout, fin, 1,in_stride, st->factors,st ); + } +} + +void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout) +{ + kiss_fft_stride(cfg,fin,fout,1); +} + + +void kiss_fft_cleanup(void) +{ + // nothing needed any more +} + +int kiss_fft_next_fast_size(int n) +{ + while(1) { + int m=n; + while ( (m%2) == 0 ) m/=2; + while ( (m%3) == 0 ) m/=3; + while ( (m%5) == 0 ) m/=5; + if (m<=1) + break; /* n is completely factorable by twos, threes, and fives */ + n++; + } + return n; +} diff --git a/software/apps/audio_module/event_detector/kiss_fftr.c b/software/apps/audio_module/event_detector/kiss_fftr.c new file mode 100644 index 00000000..b8e238b1 --- /dev/null +++ b/software/apps/audio_module/event_detector/kiss_fftr.c @@ -0,0 +1,159 @@ +/* +Copyright (c) 2003-2004, Mark Borgerding + +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. + * Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + +#include "kiss_fftr.h" +#include "_kiss_fft_guts.h" + +struct kiss_fftr_state{ + kiss_fft_cfg substate; + kiss_fft_cpx * tmpbuf; + kiss_fft_cpx * super_twiddles; +#ifdef USE_SIMD + void * pad; +#endif +}; + +kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem) +{ + int i; + kiss_fftr_cfg st = NULL; + size_t subsize, memneeded; + + if (nfft & 1) { + fprintf(stderr,"Real FFT optimization must be even.\n"); + return NULL; + } + nfft >>= 1; + + kiss_fft_alloc (nfft, inverse_fft, NULL, &subsize); + memneeded = sizeof(struct kiss_fftr_state) + subsize + sizeof(kiss_fft_cpx) * ( nfft * 3 / 2); + + if (lenmem == NULL) { + st = (kiss_fftr_cfg) KISS_FFT_MALLOC (memneeded); + } else { + if (*lenmem >= memneeded) + st = (kiss_fftr_cfg) mem; + *lenmem = memneeded; + } + if (!st) + return NULL; + + st->substate = (kiss_fft_cfg) (st + 1); /*just beyond kiss_fftr_state struct */ + st->tmpbuf = (kiss_fft_cpx *) (((char *) st->substate) + subsize); + st->super_twiddles = st->tmpbuf + nfft; + kiss_fft_alloc(nfft, inverse_fft, st->substate, &subsize); + + for (i = 0; i < nfft/2; ++i) { + double phase = + -3.14159265358979323846264338327 * ((double) (i+1) / nfft + .5); + if (inverse_fft) + phase *= -1; + kf_cexp (st->super_twiddles+i,phase); + } + return st; +} + +void kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata) +{ + /* input buffer timedata is stored row-wise */ + int k,ncfft; + kiss_fft_cpx fpnk,fpk,f1k,f2k,tw,tdc; + + if ( st->substate->inverse) { + fprintf(stderr,"kiss fft usage error: improper alloc\n"); + exit(1); + } + + ncfft = st->substate->nfft; + + /*perform the parallel fft of two real signals packed in real,imag*/ + kiss_fft( st->substate , (const kiss_fft_cpx*)timedata, st->tmpbuf ); + /* The real part of the DC element of the frequency spectrum in st->tmpbuf + * contains the sum of the even-numbered elements of the input time sequence + * The imag part is the sum of the odd-numbered elements + * + * The sum of tdc.r and tdc.i is the sum of the input time sequence. + * yielding DC of input time sequence + * The difference of tdc.r - tdc.i is the sum of the input (dot product) [1,-1,1,-1... + * yielding Nyquist bin of input time sequence + */ + + tdc.r = st->tmpbuf[0].r; + tdc.i = st->tmpbuf[0].i; + C_FIXDIV(tdc,2); + CHECK_OVERFLOW_OP(tdc.r ,+, tdc.i); + CHECK_OVERFLOW_OP(tdc.r ,-, tdc.i); + freqdata[0].r = tdc.r + tdc.i; + freqdata[ncfft].r = tdc.r - tdc.i; +#ifdef USE_SIMD + freqdata[ncfft].i = freqdata[0].i = _mm_set1_ps(0); +#else + freqdata[ncfft].i = freqdata[0].i = 0; +#endif + + for ( k=1;k <= ncfft/2 ; ++k ) { + fpk = st->tmpbuf[k]; + fpnk.r = st->tmpbuf[ncfft-k].r; + fpnk.i = - st->tmpbuf[ncfft-k].i; + C_FIXDIV(fpk,2); + C_FIXDIV(fpnk,2); + + C_ADD( f1k, fpk , fpnk ); + C_SUB( f2k, fpk , fpnk ); + C_MUL( tw , f2k , st->super_twiddles[k-1]); + + freqdata[k].r = HALF_OF(f1k.r + tw.r); + freqdata[k].i = HALF_OF(f1k.i + tw.i); + freqdata[ncfft-k].r = HALF_OF(f1k.r - tw.r); + freqdata[ncfft-k].i = HALF_OF(tw.i - f1k.i); + } +} + +void kiss_fftri(kiss_fftr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata) +{ + /* input buffer timedata is stored row-wise */ + int k, ncfft; + + if (st->substate->inverse == 0) { + fprintf (stderr, "kiss fft usage error: improper alloc\n"); + exit (1); + } + + ncfft = st->substate->nfft; + + st->tmpbuf[0].r = freqdata[0].r + freqdata[ncfft].r; + st->tmpbuf[0].i = freqdata[0].r - freqdata[ncfft].r; + C_FIXDIV(st->tmpbuf[0],2); + + for (k = 1; k <= ncfft / 2; ++k) { + kiss_fft_cpx fk, fnkc, fek, fok, tmp; + fk = freqdata[k]; + fnkc.r = freqdata[ncfft - k].r; + fnkc.i = -freqdata[ncfft - k].i; + C_FIXDIV( fk , 2 ); + C_FIXDIV( fnkc , 2 ); + + C_ADD (fek, fk, fnkc); + C_SUB (tmp, fk, fnkc); + C_MUL (fok, tmp, st->super_twiddles[k-1]); + C_ADD (st->tmpbuf[k], fek, fok); + C_SUB (st->tmpbuf[ncfft - k], fek, fok); +#ifdef USE_SIMD + st->tmpbuf[ncfft - k].i *= _mm_set1_ps(-1.0); +#else + st->tmpbuf[ncfft - k].i *= -1; +#endif + } + kiss_fft (st->substate, st->tmpbuf, (kiss_fft_cpx *) timedata); +} diff --git a/software/apps/audio_module/event_detector/log2fix.c b/software/apps/audio_module/event_detector/log2fix.c new file mode 100644 index 00000000..b39d70e1 --- /dev/null +++ b/software/apps/audio_module/event_detector/log2fix.c @@ -0,0 +1,69 @@ +#include +#include + +#include "log2fix.h" + +#define INV_LOG2_E_Q1DOT31 UINT64_C(0x58b90bfc) // Inverse log base 2 of e +#define INV_LOG2_10_Q1DOT31 UINT64_C(0x268826a1) // Inverse log base 2 of 10 + +int32_t log2fix (uint32_t x, size_t precision) +{ + // This implementation is based on Clay. S. Turner's fast binary logarithm + // algorithm[1]. + + int32_t b = 1U << (precision - 1); + int32_t y = 0; + + if (precision < 1 || precision > 31) { + errno = EINVAL; + return INT32_MAX; // indicates an error + } + + if (x == 0) { + return INT32_MIN; // represents negative infinity + } + + while (x < 1U << precision) { + x <<= 1; + y -= 1U << precision; + } + + while (x >= 2U << precision) { + x >>= 1; + y += 1U << precision; + } + + uint64_t z = x; + + for (size_t i = 0; i < precision; i++) { + z = z * z >> precision; + if (z >= 2U << precision) { + z >>= 1; + y += b; + } + b >>= 1; + } + + return y; +} + +int32_t logfix (uint32_t x, size_t precision) +{ + uint64_t t; + + t = log2fix(x, precision) * INV_LOG2_E_Q1DOT31; + + return t >> 31; +} + +int32_t log10fix (uint32_t x, size_t precision) +{ + uint64_t t; + + t = log2fix(x, precision) * INV_LOG2_10_Q1DOT31; + + return t >> 31; +} + +// [1] C. S. Turner, "A Fast Binary Logarithm Algorithm", IEEE Signal +// Processing Mag., pp. 124,140, Sep. 2010. diff --git a/software/apps/audio_module/event_detector/ridgeTracker.c b/software/apps/audio_module/event_detector/ridgeTracker.c new file mode 100644 index 00000000..8ad80403 --- /dev/null +++ b/software/apps/audio_module/event_detector/ridgeTracker.c @@ -0,0 +1,184 @@ +/* + * Ridge tracker + * + * Long Le + * University of Illinois + */ + +#include +#include + +#include "_kiss_fft_guts.h" +#include "common.h" +#include "log2fix/log2fix.h" +#include "dynArray.h" + +#define FOFF (2) +// SAMP_MAX is fixed-point 1 +#define BTLEN (10) // int(btTime/tInc) +#define ALP ((kiss_fft_scalar)(0.9*SAMP_MAX)) // exp(-tInc/btTime) +#define SUPTHRESH ((kiss_fft_scalar)(0.3*SAMP_MAX)) +#define NADA ((kiss_fft_scalar)(0.01*SAMP_MAX)) +#define EPSILON ((kiss_fft_scalar)(0.001*SAMP_MAX)) //noise adaptive param +kiss_fft_scalar *ind=NULL; +kiss_fft_scalar *noiseFloor=NULL; +kiss_fft_scalar *snrAcc=NULL; +kiss_fft_scalar *freAcc=NULL; + +#define PROBTARGET ((kiss_fft_scalar)(0.5*SAMP_MAX)) +#define MUP ((kiss_fft_scalar)(0.99*SAMP_MAX)) +#define MUT ((kiss_fft_scalar)(0.99*SAMP_MAX)) +int curTime; +kiss_fft_scalar adaThresh; +kiss_fft_scalar probEst; + +bool ridgeTracker_isReady; +Array ridgeTracker_out; + +void ridgeTracker_init(void){ + ind = malloc(FRE_LEN*sizeof(kiss_fft_scalar)); + noiseFloor = malloc(FRE_LEN*sizeof(kiss_fft_scalar)); + snrAcc = malloc(FRE_LEN*sizeof(kiss_fft_scalar)); + freAcc = malloc(FRE_LEN*sizeof(kiss_fft_scalar)); + for (size_t f=0; f noiseFloor[f]){ + ind[f] -= 1; + if (ind[f] < 0){ + noiseFloor[f] += S_MUL(noiseFloor[f],NADA); + }else{ + noiseFloor[f] += S_MUL(noiseFloor[f],NADA/2); + } + }else{ + noiseFloor[f] = MAX(EPSILON, noiseFloor[f]-S_MUL(noiseFloor[f],NADA)); + ind[f] = BTLEN; + } + + // snr update + uint32_t r = spec[f]/noiseFloor[f]; + if (r >= 10){ + snr = SAMP_MAX; + }else if (r > 0){ + // see log2fix/main.c + snr = log10fix(r*(1< maxVal){ + maxVal = val; + maxIdx = l; + } + } + snrAcc[f] = maxVal; + freAcc[f] = maxIdx; + } + + // max-pooling + memcpy(snrAccLast, snrAcc, FRE_LEN*sizeof(kiss_fft_scalar)); + for (size_t f=0; f maxVal){ + maxVal = snrAccLast[l]; + maxIdx = l; + } + } + if (maxIdx == f){ + snrAcc[f] = maxVal; + }else{ + snrAcc[f] = 0; + } + } + + // per-bin suppressed stats + memset(snrOut, 0, FRE_LEN*sizeof(kiss_fft_scalar)); + bool binActive = false; + for (size_t f=0; f SUPTHRESH){ + binActive = true; + snrOut[f] = snrAcc[f] - SUPTHRESH; + + if (curTime == -1){ + printf("\nStart timer\n"); + curTime = BTLEN; + } + insertArray(&ridgeTracker_out,snrOut[f],f,(size_t)curTime); + } + } + + // output if applicable + if (curTime != -1){ + if (!binActive){ + printf("Potential output\n"); + // longer time-granuality approximation + if (getMaxDurArray(ridgeTracker_out) >= 50){ + printf("Sufficiently long temporal approximation\n"); + // adaptive thresholding + printf("avgSNR = %d\n", getAvgSNRArray(ridgeTracker_out)); + if (getAvgSNRArray(ridgeTracker_out) >= adaThresh){ + ridgeTracker_isReady = true; + probEst = S_MUL(MUP,probEst) + (SAMP_MAX-MUP); + } else{ + ridgeTracker_reset(); + probEst = S_MUL(MUP,probEst); + } + } else{ + ridgeTracker_reset(); + } + adaThresh = adaThresh + S_MUL(SAMP_MAX-MUT,probEst-PROBTARGET); + + } else{ + curTime += 1; + } + } + + // cleanup + free(snrAccLast); +} + +void ridgeTracker_destroy(void){ + free(ind); + free(noiseFloor); + free(snrAcc); + free(freAcc); + + printf("ridgeTracker_destroy() completed\n"); +} diff --git a/software/apps/audio_module/event_detector/ridgeTracker.h b/software/apps/audio_module/event_detector/ridgeTracker.h new file mode 100644 index 00000000..114eab1a --- /dev/null +++ b/software/apps/audio_module/event_detector/ridgeTracker.h @@ -0,0 +1,21 @@ +// FIXED POINT MATHS +// Long Le +// University of Illinois +// + +#ifndef __RIDGETRACKER_H__ +#define __RIDGETRACKER_H__ + +#include +#include "dynArray.h" + +// The extern vars +extern bool ridgeTracker_isReady; +extern Array ridgeTracker_out; + +void ridgeTracker_init(void); +void ridgeTracker_update(kiss_fft_scalar* spec, kiss_fft_scalar* snrOut); +void ridgeTracker_destroy(void); +void ridgeTracker_reset(void); + +#endif From f16db9fedf6e45cf09d891d1a6b84d7f9bbf3dd4 Mon Sep 17 00:00:00 2001 From: "Long N. Le" Date: Tue, 7 Mar 2017 11:00:21 -0600 Subject: [PATCH 3/8] updated README.md --- software/apps/audio_module/event_detector/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/software/apps/audio_module/event_detector/README.md b/software/apps/audio_module/event_detector/README.md index 1041d620..c35aafac 100644 --- a/software/apps/audio_module/event_detector/README.md +++ b/software/apps/audio_module/event_detector/README.md @@ -1,4 +1,4 @@ -# Fixed-point ridge tracker in C +# Fixed-point ridge tracker for acoustic event detection in C ## Quick start ``` From ed2315ed7a80149e1df38b29892bdc3b0d80b7a5 Mon Sep 17 00:00:00 2001 From: Long Le Date: Tue, 7 Mar 2017 16:35:03 -0600 Subject: [PATCH 4/8] completed main.c --- .../apps/audio_module/event_detector/main.c | 83 +++++++++++++++++-- 1 file changed, 78 insertions(+), 5 deletions(-) diff --git a/software/apps/audio_module/event_detector/main.c b/software/apps/audio_module/event_detector/main.c index f3e25d83..8570acf2 100644 --- a/software/apps/audio_module/event_detector/main.c +++ b/software/apps/audio_module/event_detector/main.c @@ -21,6 +21,8 @@ #include "ridgeTracker.h" static int freqAnalyze(kiss_fft_scalar *frame, kiss_fft_scalar *spec); +static int printMat(char *filename, kiss_fft_scalar *in, size_t N); +static int printMat2(char *filename, size_t *in, size_t N); int main (void) { int err = SUCCESS; @@ -41,13 +43,16 @@ int main (void) { ridgeTracker_init(); while (true) { // read data from ADC - err = adc_read_single_sample(3); - if (err < SUCCESS) { - printf("ADC read error: %d\n", err); + for (size_t k = frameIdx; k < frameIdx+INC_LEN; k++){ + err = adc_read_single_sample(3); + if (err < SUCCESS) { + printf("ADC read error: %d\n", err); + } + uint16_t sample = err & 0xFFFF; + frame[k] = (kiss_fft_scalar)sample; } - uint16_t sample = err & 0xFFFF; - frameIdx = (frameIdx+INC_LEN) % BUF_LEN; + // perform FFT if (freqAnalyze(frame, spec)){ printf("freqAnalyze() failed!\n"); @@ -56,7 +61,39 @@ int main (void) { // ridgetracker update ridgeTracker_update(spec, snrOut); + + // check for available event + if (ridgeTracker_isReady){ + printf("Event detected!\n"); + printf("ridgeTracker_out.used = %zu\n",ridgeTracker_out.used); + + if (printMat("SNR.mat",ridgeTracker_out.SNR,ridgeTracker_out.used)){ + printf("printMat() failed!\n"); + return 1; + } + if (printMat2("FI.mat",ridgeTracker_out.FI,ridgeTracker_out.used)){ + printf("printMat2() failed!\n"); + return 1; + } + size_t TIE = ridgeTracker_out.TI[ridgeTracker_out.used-1]; + for (size_t k=0; k Date: Sun, 12 Mar 2017 10:43:08 -0500 Subject: [PATCH 5/8] updated README.md according to the Amit's review --- .../audio_module/event_detector/README.md | 27 +++++++++++++------ 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/software/apps/audio_module/event_detector/README.md b/software/apps/audio_module/event_detector/README.md index c35aafac..6c69be95 100644 --- a/software/apps/audio_module/event_detector/README.md +++ b/software/apps/audio_module/event_detector/README.md @@ -1,14 +1,25 @@ -# Fixed-point ridge tracker for acoustic event detection in C +# Universal acoustic event detection app +The first step into any audio processing chain: + +Detect a large class of "acoustic events" (bird calls, gunshots, broken glasses, etc.) from noisy background. +Specific classification are reserved for downstream. ## Quick start -``` -make clean -make -./main input.wav -``` -Use test.ipynb to visualize the resulting \*.mat files +1. Plug the Audio Module into the `Module 1` slot + +2. **Important** Make sure the programming knob is turned to `MOD1`. + +3. Get all the prerequisites (See below) + +4. Flash the app + + ```bash + cd signpost/software/apps/audio_module/event_detector + make flash + ``` ## Prerequisites * Fixed-point FFT: git clone https://github.com/longle2718/kiss_fft * Fixed-point log: git clone https://github.com/dmoulding/log2fix -* libsndfile: sudo apt-get install libsndfile1-dev +* Ridge tracker: git clone https://bitbucket.org/longle1/gcwa +..* Fixed point implementation under gcwa/c/ From 972b5c2876c16715acf195d158d6b10f328d04a2 Mon Sep 17 00:00:00 2001 From: Long Le Date: Sun, 12 Mar 2017 11:24:20 -0500 Subject: [PATCH 6/8] Update README.md --- software/apps/audio_module/event_detector/README.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/software/apps/audio_module/event_detector/README.md b/software/apps/audio_module/event_detector/README.md index 6c69be95..b478e9e1 100644 --- a/software/apps/audio_module/event_detector/README.md +++ b/software/apps/audio_module/event_detector/README.md @@ -1,8 +1,8 @@ # Universal acoustic event detection app -The first step into any audio processing chain: +The first step into any audio processing chain. Detect a large class of "acoustic events" (bird calls, gunshots, broken glasses, etc.) from noisy background. -Specific classification are reserved for downstream. +Further classifications are reserved for downstream processing. ## Quick start 1. Plug the Audio Module into the `Module 1` slot @@ -22,4 +22,4 @@ Specific classification are reserved for downstream. * Fixed-point FFT: git clone https://github.com/longle2718/kiss_fft * Fixed-point log: git clone https://github.com/dmoulding/log2fix * Ridge tracker: git clone https://bitbucket.org/longle1/gcwa -..* Fixed point implementation under gcwa/c/ + * Fixed point implementation under gcwa/c/ From 8bfec1e25ff6abe14107f241e765f85c7adda464 Mon Sep 17 00:00:00 2001 From: Long Le Date: Sun, 12 Mar 2017 11:24:55 -0500 Subject: [PATCH 7/8] Update README.md --- software/apps/audio_module/event_detector/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/software/apps/audio_module/event_detector/README.md b/software/apps/audio_module/event_detector/README.md index b478e9e1..59887354 100644 --- a/software/apps/audio_module/event_detector/README.md +++ b/software/apps/audio_module/event_detector/README.md @@ -22,4 +22,4 @@ Further classifications are reserved for downstream processing. * Fixed-point FFT: git clone https://github.com/longle2718/kiss_fft * Fixed-point log: git clone https://github.com/dmoulding/log2fix * Ridge tracker: git clone https://bitbucket.org/longle1/gcwa - * Fixed point implementation under gcwa/c/ + * Fixed-point implementation in C is under gcwa/c/ From a6bba7819adbe3d2a0c3b7dad1d765a7fc2ec5aa Mon Sep 17 00:00:00 2001 From: Long Le Date: Sun, 12 Mar 2017 11:30:27 -0500 Subject: [PATCH 8/8] working on making a hierarchical dir struct --- .../apps/audio_module/event_detector/Makefile | 9 +- .../apps/audio_module/event_detector/common.h | 21 - .../audio_module/event_detector/dynArray.c | 64 --- .../audio_module/event_detector/dynArray.h | 21 - .../audio_module/event_detector/kiss_fft.c | 408 ------------------ .../audio_module/event_detector/kiss_fftr.c | 159 ------- .../audio_module/event_detector/log2fix.c | 69 --- .../event_detector/ridgeTracker.c | 184 -------- .../event_detector/ridgeTracker.h | 21 - 9 files changed, 8 insertions(+), 948 deletions(-) delete mode 100644 software/apps/audio_module/event_detector/common.h delete mode 100644 software/apps/audio_module/event_detector/dynArray.c delete mode 100644 software/apps/audio_module/event_detector/dynArray.h delete mode 100644 software/apps/audio_module/event_detector/kiss_fft.c delete mode 100644 software/apps/audio_module/event_detector/kiss_fftr.c delete mode 100644 software/apps/audio_module/event_detector/log2fix.c delete mode 100644 software/apps/audio_module/event_detector/ridgeTracker.c delete mode 100644 software/apps/audio_module/event_detector/ridgeTracker.h diff --git a/software/apps/audio_module/event_detector/Makefile b/software/apps/audio_module/event_detector/Makefile index e32011b8..79be72e1 100644 --- a/software/apps/audio_module/event_detector/Makefile +++ b/software/apps/audio_module/event_detector/Makefile @@ -5,7 +5,14 @@ APP_DIR := $(dir $(abspath $(lastword $(MAKEFILE_LIST)))) # files needed for this code C_SRCS := $(wildcard *.c) -INCLUDE_PATHS += . ./kiss_fft ./kiss_fft/tools ./log2fix/ +C_SRCS := $(wildcard **/log2fix.c) +#C_SRCS += ./log2fix/log2fix.c +#C_SRCS += ./kiss_fft/kiss_fft.c +#C_SRCS += ./kiss_fft/tools/kiss_fftr.c +#C_SRCS += ./gcwa/c/ridgeTracker.c +#C_SRCS += ./gcwa/c/dynArray.c + +INCLUDE_PATHS += . ./kiss_fft ./kiss_fft/tools ./log2fix ./gcwa/c CFLAGS += -DFIXED_POINT=16 -std=c99 CFLAGS += -Wno-type-limits -Wno-sign-compare diff --git a/software/apps/audio_module/event_detector/common.h b/software/apps/audio_module/event_detector/common.h deleted file mode 100644 index 8a15780d..00000000 --- a/software/apps/audio_module/event_detector/common.h +++ /dev/null @@ -1,21 +0,0 @@ -#ifndef __COMMON_H__ -#define __COMMON_H__ - -#include "_kiss_fft_guts.h" - -#define BUF_LEN (512) -#define INC_LEN (BUF_LEN/2) -#define FRE_LEN (BUF_LEN/2) - -#define S_SQ(a) S_MUL((a),(a)) // scalar square -#define S_DIV(a,b) S_MUL((a),SAMP_MAX/(b)) // a < b -#define MAX(a,b) (((a) > (b)) ? (a) : (b)) -#define MIN(a,b) (((a) < (b)) ? (a) : (b)) -#define ABS(a) (((a) < 0) ? -(a): (a)) -// Magnitude Estimator -// http://dspguru.com/book/export/html/62 -#define ALPHA ((kiss_fft_scalar)(0.947543636291*SAMP_MAX)) -#define BETA ((kiss_fft_scalar)(0.392485425092*SAMP_MAX)) -#define MAG(r,i) ( S_MUL(ALPHA,MAX(ABS((r)),ABS((i)))) + S_MUL(BETA,MIN(ABS((r)),ABS((i)))) ) - -#endif diff --git a/software/apps/audio_module/event_detector/dynArray.c b/software/apps/audio_module/event_detector/dynArray.c deleted file mode 100644 index 07d82b65..00000000 --- a/software/apps/audio_module/event_detector/dynArray.c +++ /dev/null @@ -1,64 +0,0 @@ -#include -#include -#include "dynArray.h" -#include "_kiss_fft_guts.h" - -void initArray(Array *a, size_t initialSize) { - a->maxSize = initialSize; - a->used = 0; - - a->size = a->maxSize; - printf("Init array size = %zu\n",a->size); - a->SNR = malloc(a->size * sizeof(kiss_fft_scalar)); - a->FI = malloc(a->size * sizeof(size_t)); - a->TI = malloc(a->size * sizeof(size_t)); -} - -void insertArray(Array *a, kiss_fft_scalar snr, size_t fi, size_t ti) { - // a->used is the number of used entries, because a->array[a->used++] updates a->used only *after* the array has been accessed. - // Therefore a->used can go up to a->size - if (a->size > 0 && a->used == a->size) { - a->maxSize *= 2; - printf("Doubled array max size to %zu\n",a->maxSize); - } - if (a->size < a->maxSize) { - a->size = a->maxSize; - // with NULL ptr, realloc == malloc - a->SNR = realloc(a->SNR, a->size * sizeof(kiss_fft_scalar)); - a->FI = realloc(a->FI, a->size * sizeof(size_t)); - a->TI = realloc(a->TI, a->size * sizeof(size_t)); - } - a->SNR[a->used] = snr; - a->FI[a->used] = fi; - a->TI[a->used] = ti; - a->used++; -} - -void freeArray(Array *a) { - free(a->SNR); - free(a->FI); - free(a->TI); - a->SNR = NULL; - a->FI = NULL; - a->TI = NULL; - a->used = a->size = 0; -} -size_t getMaxDurArray(Array a){ - if (a.used == 0) - return 0; - - size_t minTI = a.TI[0]; - size_t maxTI = a.TI[a.used-1]; - return maxTI-minTI; -} - -kiss_fft_scalar getAvgSNRArray(Array a){ - if (a.used == 0) - return 0; - - SAMPPROD sum = 0; - for (size_t k=0; k < a.used; k++){ - sum += a.SNR[k]; - } - return (kiss_fft_scalar)(sum/a.used); -} diff --git a/software/apps/audio_module/event_detector/dynArray.h b/software/apps/audio_module/event_detector/dynArray.h deleted file mode 100644 index 47918f37..00000000 --- a/software/apps/audio_module/event_detector/dynArray.h +++ /dev/null @@ -1,21 +0,0 @@ -#ifndef __DYNARRAY_H__ -#define __DYNARRAY_H__ - -#include "kiss_fft.h" - -typedef struct { - kiss_fft_scalar *SNR; // log SNR - size_t *FI; // frequency index - size_t *TI; // time index - size_t size; - size_t used; - size_t maxSize; -} Array; - -void initArray(Array *a, size_t initialSize); -void insertArray(Array *a, kiss_fft_scalar snr, size_t fi, size_t ti); -void freeArray(Array *a); - -size_t getMaxDurArray(Array a); -kiss_fft_scalar getAvgSNRArray(Array a); -#endif diff --git a/software/apps/audio_module/event_detector/kiss_fft.c b/software/apps/audio_module/event_detector/kiss_fft.c deleted file mode 100644 index 465d6c97..00000000 --- a/software/apps/audio_module/event_detector/kiss_fft.c +++ /dev/null @@ -1,408 +0,0 @@ -/* -Copyright (c) 2003-2010, Mark Borgerding - -All rights reserved. - -Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: - - * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. - * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. - * Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission. - -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -*/ - - -#include "_kiss_fft_guts.h" -/* The guts header contains all the multiplication and addition macros that are defined for - fixed or floating point complex numbers. It also delares the kf_ internal functions. - */ - -static void kf_bfly2( - kiss_fft_cpx * Fout, - const size_t fstride, - const kiss_fft_cfg st, - int m - ) -{ - kiss_fft_cpx * Fout2; - kiss_fft_cpx * tw1 = st->twiddles; - kiss_fft_cpx t; - Fout2 = Fout + m; - do{ - C_FIXDIV(*Fout,2); C_FIXDIV(*Fout2,2); - - C_MUL (t, *Fout2 , *tw1); - tw1 += fstride; - C_SUB( *Fout2 , *Fout , t ); - C_ADDTO( *Fout , t ); - ++Fout2; - ++Fout; - }while (--m); -} - -static void kf_bfly4( - kiss_fft_cpx * Fout, - const size_t fstride, - const kiss_fft_cfg st, - const size_t m - ) -{ - kiss_fft_cpx *tw1,*tw2,*tw3; - kiss_fft_cpx scratch[6]; - size_t k=m; - const size_t m2=2*m; - const size_t m3=3*m; - - - tw3 = tw2 = tw1 = st->twiddles; - - do { - C_FIXDIV(*Fout,4); C_FIXDIV(Fout[m],4); C_FIXDIV(Fout[m2],4); C_FIXDIV(Fout[m3],4); - - C_MUL(scratch[0],Fout[m] , *tw1 ); - C_MUL(scratch[1],Fout[m2] , *tw2 ); - C_MUL(scratch[2],Fout[m3] , *tw3 ); - - C_SUB( scratch[5] , *Fout, scratch[1] ); - C_ADDTO(*Fout, scratch[1]); - C_ADD( scratch[3] , scratch[0] , scratch[2] ); - C_SUB( scratch[4] , scratch[0] , scratch[2] ); - C_SUB( Fout[m2], *Fout, scratch[3] ); - tw1 += fstride; - tw2 += fstride*2; - tw3 += fstride*3; - C_ADDTO( *Fout , scratch[3] ); - - if(st->inverse) { - Fout[m].r = scratch[5].r - scratch[4].i; - Fout[m].i = scratch[5].i + scratch[4].r; - Fout[m3].r = scratch[5].r + scratch[4].i; - Fout[m3].i = scratch[5].i - scratch[4].r; - }else{ - Fout[m].r = scratch[5].r + scratch[4].i; - Fout[m].i = scratch[5].i - scratch[4].r; - Fout[m3].r = scratch[5].r - scratch[4].i; - Fout[m3].i = scratch[5].i + scratch[4].r; - } - ++Fout; - }while(--k); -} - -static void kf_bfly3( - kiss_fft_cpx * Fout, - const size_t fstride, - const kiss_fft_cfg st, - size_t m - ) -{ - size_t k=m; - const size_t m2 = 2*m; - kiss_fft_cpx *tw1,*tw2; - kiss_fft_cpx scratch[5]; - kiss_fft_cpx epi3; - epi3 = st->twiddles[fstride*m]; - - tw1=tw2=st->twiddles; - - do{ - C_FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3); - - C_MUL(scratch[1],Fout[m] , *tw1); - C_MUL(scratch[2],Fout[m2] , *tw2); - - C_ADD(scratch[3],scratch[1],scratch[2]); - C_SUB(scratch[0],scratch[1],scratch[2]); - tw1 += fstride; - tw2 += fstride*2; - - Fout[m].r = Fout->r - HALF_OF(scratch[3].r); - Fout[m].i = Fout->i - HALF_OF(scratch[3].i); - - C_MULBYSCALAR( scratch[0] , epi3.i ); - - C_ADDTO(*Fout,scratch[3]); - - Fout[m2].r = Fout[m].r + scratch[0].i; - Fout[m2].i = Fout[m].i - scratch[0].r; - - Fout[m].r -= scratch[0].i; - Fout[m].i += scratch[0].r; - - ++Fout; - }while(--k); -} - -static void kf_bfly5( - kiss_fft_cpx * Fout, - const size_t fstride, - const kiss_fft_cfg st, - int m - ) -{ - kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4; - int u; - kiss_fft_cpx scratch[13]; - kiss_fft_cpx * twiddles = st->twiddles; - kiss_fft_cpx *tw; - kiss_fft_cpx ya,yb; - ya = twiddles[fstride*m]; - yb = twiddles[fstride*2*m]; - - Fout0=Fout; - Fout1=Fout0+m; - Fout2=Fout0+2*m; - Fout3=Fout0+3*m; - Fout4=Fout0+4*m; - - tw=st->twiddles; - for ( u=0; ur += scratch[7].r + scratch[8].r; - Fout0->i += scratch[7].i + scratch[8].i; - - scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r); - scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r); - - scratch[6].r = S_MUL(scratch[10].i,ya.i) + S_MUL(scratch[9].i,yb.i); - scratch[6].i = -S_MUL(scratch[10].r,ya.i) - S_MUL(scratch[9].r,yb.i); - - C_SUB(*Fout1,scratch[5],scratch[6]); - C_ADD(*Fout4,scratch[5],scratch[6]); - - scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r); - scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r); - scratch[12].r = - S_MUL(scratch[10].i,yb.i) + S_MUL(scratch[9].i,ya.i); - scratch[12].i = S_MUL(scratch[10].r,yb.i) - S_MUL(scratch[9].r,ya.i); - - C_ADD(*Fout2,scratch[11],scratch[12]); - C_SUB(*Fout3,scratch[11],scratch[12]); - - ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4; - } -} - -/* perform the butterfly for one stage of a mixed radix FFT */ -static void kf_bfly_generic( - kiss_fft_cpx * Fout, - const size_t fstride, - const kiss_fft_cfg st, - int m, - int p - ) -{ - int u,k,q1,q; - kiss_fft_cpx * twiddles = st->twiddles; - kiss_fft_cpx t; - int Norig = st->nfft; - - kiss_fft_cpx * scratch = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC(sizeof(kiss_fft_cpx)*p); - - for ( u=0; u=Norig) twidx-=Norig; - C_MUL(t,scratch[q] , twiddles[twidx] ); - C_ADDTO( Fout[ k ] ,t); - } - k += m; - } - } - KISS_FFT_TMP_FREE(scratch); -} - -static -void kf_work( - kiss_fft_cpx * Fout, - const kiss_fft_cpx * f, - const size_t fstride, - int in_stride, - int * factors, - const kiss_fft_cfg st - ) -{ - kiss_fft_cpx * Fout_beg=Fout; - const int p=*factors++; /* the radix */ - const int m=*factors++; /* stage's fft length/p */ - const kiss_fft_cpx * Fout_end = Fout + p*m; - -#ifdef _OPENMP - // use openmp extensions at the - // top-level (not recursive) - if (fstride==1 && p<=5) - { - int k; - - // execute the p different work units in different threads -# pragma omp parallel for - for (k=0;k floor_sqrt) - p = n; /* no more factors, skip to end */ - } - n /= p; - *facbuf++ = p; - *facbuf++ = n; - } while (n > 1); -} - -/* - * - * User-callable function to allocate all necessary storage space for the fft. - * - * The return value is a contiguous block of memory, allocated with malloc. As such, - * It can be freed with free(), rather than a kiss_fft-specific function. - * */ -kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem ) -{ - kiss_fft_cfg st=NULL; - size_t memneeded = sizeof(struct kiss_fft_state) - + sizeof(kiss_fft_cpx)*(nfft-1); /* twiddle factors*/ - - if ( lenmem==NULL ) { - st = ( kiss_fft_cfg)KISS_FFT_MALLOC( memneeded ); - }else{ - if (mem != NULL && *lenmem >= memneeded) - st = (kiss_fft_cfg)mem; - *lenmem = memneeded; - } - if (st) { - int i; - st->nfft=nfft; - st->inverse = inverse_fft; - - for (i=0;iinverse) - phase *= -1; - kf_cexp(st->twiddles+i, phase ); - } - - kf_factor(nfft,st->factors); - } - return st; -} - - -void kiss_fft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride) -{ - if (fin == fout) { - //NOTE: this is not really an in-place FFT algorithm. - //It just performs an out-of-place FFT into a temp buffer - kiss_fft_cpx * tmpbuf = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC( sizeof(kiss_fft_cpx)*st->nfft); - kf_work(tmpbuf,fin,1,in_stride, st->factors,st); - memcpy(fout,tmpbuf,sizeof(kiss_fft_cpx)*st->nfft); - KISS_FFT_TMP_FREE(tmpbuf); - }else{ - kf_work( fout, fin, 1,in_stride, st->factors,st ); - } -} - -void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout) -{ - kiss_fft_stride(cfg,fin,fout,1); -} - - -void kiss_fft_cleanup(void) -{ - // nothing needed any more -} - -int kiss_fft_next_fast_size(int n) -{ - while(1) { - int m=n; - while ( (m%2) == 0 ) m/=2; - while ( (m%3) == 0 ) m/=3; - while ( (m%5) == 0 ) m/=5; - if (m<=1) - break; /* n is completely factorable by twos, threes, and fives */ - n++; - } - return n; -} diff --git a/software/apps/audio_module/event_detector/kiss_fftr.c b/software/apps/audio_module/event_detector/kiss_fftr.c deleted file mode 100644 index b8e238b1..00000000 --- a/software/apps/audio_module/event_detector/kiss_fftr.c +++ /dev/null @@ -1,159 +0,0 @@ -/* -Copyright (c) 2003-2004, Mark Borgerding - -All rights reserved. - -Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: - - * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. - * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. - * Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission. - -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -*/ - -#include "kiss_fftr.h" -#include "_kiss_fft_guts.h" - -struct kiss_fftr_state{ - kiss_fft_cfg substate; - kiss_fft_cpx * tmpbuf; - kiss_fft_cpx * super_twiddles; -#ifdef USE_SIMD - void * pad; -#endif -}; - -kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem) -{ - int i; - kiss_fftr_cfg st = NULL; - size_t subsize, memneeded; - - if (nfft & 1) { - fprintf(stderr,"Real FFT optimization must be even.\n"); - return NULL; - } - nfft >>= 1; - - kiss_fft_alloc (nfft, inverse_fft, NULL, &subsize); - memneeded = sizeof(struct kiss_fftr_state) + subsize + sizeof(kiss_fft_cpx) * ( nfft * 3 / 2); - - if (lenmem == NULL) { - st = (kiss_fftr_cfg) KISS_FFT_MALLOC (memneeded); - } else { - if (*lenmem >= memneeded) - st = (kiss_fftr_cfg) mem; - *lenmem = memneeded; - } - if (!st) - return NULL; - - st->substate = (kiss_fft_cfg) (st + 1); /*just beyond kiss_fftr_state struct */ - st->tmpbuf = (kiss_fft_cpx *) (((char *) st->substate) + subsize); - st->super_twiddles = st->tmpbuf + nfft; - kiss_fft_alloc(nfft, inverse_fft, st->substate, &subsize); - - for (i = 0; i < nfft/2; ++i) { - double phase = - -3.14159265358979323846264338327 * ((double) (i+1) / nfft + .5); - if (inverse_fft) - phase *= -1; - kf_cexp (st->super_twiddles+i,phase); - } - return st; -} - -void kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata) -{ - /* input buffer timedata is stored row-wise */ - int k,ncfft; - kiss_fft_cpx fpnk,fpk,f1k,f2k,tw,tdc; - - if ( st->substate->inverse) { - fprintf(stderr,"kiss fft usage error: improper alloc\n"); - exit(1); - } - - ncfft = st->substate->nfft; - - /*perform the parallel fft of two real signals packed in real,imag*/ - kiss_fft( st->substate , (const kiss_fft_cpx*)timedata, st->tmpbuf ); - /* The real part of the DC element of the frequency spectrum in st->tmpbuf - * contains the sum of the even-numbered elements of the input time sequence - * The imag part is the sum of the odd-numbered elements - * - * The sum of tdc.r and tdc.i is the sum of the input time sequence. - * yielding DC of input time sequence - * The difference of tdc.r - tdc.i is the sum of the input (dot product) [1,-1,1,-1... - * yielding Nyquist bin of input time sequence - */ - - tdc.r = st->tmpbuf[0].r; - tdc.i = st->tmpbuf[0].i; - C_FIXDIV(tdc,2); - CHECK_OVERFLOW_OP(tdc.r ,+, tdc.i); - CHECK_OVERFLOW_OP(tdc.r ,-, tdc.i); - freqdata[0].r = tdc.r + tdc.i; - freqdata[ncfft].r = tdc.r - tdc.i; -#ifdef USE_SIMD - freqdata[ncfft].i = freqdata[0].i = _mm_set1_ps(0); -#else - freqdata[ncfft].i = freqdata[0].i = 0; -#endif - - for ( k=1;k <= ncfft/2 ; ++k ) { - fpk = st->tmpbuf[k]; - fpnk.r = st->tmpbuf[ncfft-k].r; - fpnk.i = - st->tmpbuf[ncfft-k].i; - C_FIXDIV(fpk,2); - C_FIXDIV(fpnk,2); - - C_ADD( f1k, fpk , fpnk ); - C_SUB( f2k, fpk , fpnk ); - C_MUL( tw , f2k , st->super_twiddles[k-1]); - - freqdata[k].r = HALF_OF(f1k.r + tw.r); - freqdata[k].i = HALF_OF(f1k.i + tw.i); - freqdata[ncfft-k].r = HALF_OF(f1k.r - tw.r); - freqdata[ncfft-k].i = HALF_OF(tw.i - f1k.i); - } -} - -void kiss_fftri(kiss_fftr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata) -{ - /* input buffer timedata is stored row-wise */ - int k, ncfft; - - if (st->substate->inverse == 0) { - fprintf (stderr, "kiss fft usage error: improper alloc\n"); - exit (1); - } - - ncfft = st->substate->nfft; - - st->tmpbuf[0].r = freqdata[0].r + freqdata[ncfft].r; - st->tmpbuf[0].i = freqdata[0].r - freqdata[ncfft].r; - C_FIXDIV(st->tmpbuf[0],2); - - for (k = 1; k <= ncfft / 2; ++k) { - kiss_fft_cpx fk, fnkc, fek, fok, tmp; - fk = freqdata[k]; - fnkc.r = freqdata[ncfft - k].r; - fnkc.i = -freqdata[ncfft - k].i; - C_FIXDIV( fk , 2 ); - C_FIXDIV( fnkc , 2 ); - - C_ADD (fek, fk, fnkc); - C_SUB (tmp, fk, fnkc); - C_MUL (fok, tmp, st->super_twiddles[k-1]); - C_ADD (st->tmpbuf[k], fek, fok); - C_SUB (st->tmpbuf[ncfft - k], fek, fok); -#ifdef USE_SIMD - st->tmpbuf[ncfft - k].i *= _mm_set1_ps(-1.0); -#else - st->tmpbuf[ncfft - k].i *= -1; -#endif - } - kiss_fft (st->substate, st->tmpbuf, (kiss_fft_cpx *) timedata); -} diff --git a/software/apps/audio_module/event_detector/log2fix.c b/software/apps/audio_module/event_detector/log2fix.c deleted file mode 100644 index b39d70e1..00000000 --- a/software/apps/audio_module/event_detector/log2fix.c +++ /dev/null @@ -1,69 +0,0 @@ -#include -#include - -#include "log2fix.h" - -#define INV_LOG2_E_Q1DOT31 UINT64_C(0x58b90bfc) // Inverse log base 2 of e -#define INV_LOG2_10_Q1DOT31 UINT64_C(0x268826a1) // Inverse log base 2 of 10 - -int32_t log2fix (uint32_t x, size_t precision) -{ - // This implementation is based on Clay. S. Turner's fast binary logarithm - // algorithm[1]. - - int32_t b = 1U << (precision - 1); - int32_t y = 0; - - if (precision < 1 || precision > 31) { - errno = EINVAL; - return INT32_MAX; // indicates an error - } - - if (x == 0) { - return INT32_MIN; // represents negative infinity - } - - while (x < 1U << precision) { - x <<= 1; - y -= 1U << precision; - } - - while (x >= 2U << precision) { - x >>= 1; - y += 1U << precision; - } - - uint64_t z = x; - - for (size_t i = 0; i < precision; i++) { - z = z * z >> precision; - if (z >= 2U << precision) { - z >>= 1; - y += b; - } - b >>= 1; - } - - return y; -} - -int32_t logfix (uint32_t x, size_t precision) -{ - uint64_t t; - - t = log2fix(x, precision) * INV_LOG2_E_Q1DOT31; - - return t >> 31; -} - -int32_t log10fix (uint32_t x, size_t precision) -{ - uint64_t t; - - t = log2fix(x, precision) * INV_LOG2_10_Q1DOT31; - - return t >> 31; -} - -// [1] C. S. Turner, "A Fast Binary Logarithm Algorithm", IEEE Signal -// Processing Mag., pp. 124,140, Sep. 2010. diff --git a/software/apps/audio_module/event_detector/ridgeTracker.c b/software/apps/audio_module/event_detector/ridgeTracker.c deleted file mode 100644 index 8ad80403..00000000 --- a/software/apps/audio_module/event_detector/ridgeTracker.c +++ /dev/null @@ -1,184 +0,0 @@ -/* - * Ridge tracker - * - * Long Le - * University of Illinois - */ - -#include -#include - -#include "_kiss_fft_guts.h" -#include "common.h" -#include "log2fix/log2fix.h" -#include "dynArray.h" - -#define FOFF (2) -// SAMP_MAX is fixed-point 1 -#define BTLEN (10) // int(btTime/tInc) -#define ALP ((kiss_fft_scalar)(0.9*SAMP_MAX)) // exp(-tInc/btTime) -#define SUPTHRESH ((kiss_fft_scalar)(0.3*SAMP_MAX)) -#define NADA ((kiss_fft_scalar)(0.01*SAMP_MAX)) -#define EPSILON ((kiss_fft_scalar)(0.001*SAMP_MAX)) //noise adaptive param -kiss_fft_scalar *ind=NULL; -kiss_fft_scalar *noiseFloor=NULL; -kiss_fft_scalar *snrAcc=NULL; -kiss_fft_scalar *freAcc=NULL; - -#define PROBTARGET ((kiss_fft_scalar)(0.5*SAMP_MAX)) -#define MUP ((kiss_fft_scalar)(0.99*SAMP_MAX)) -#define MUT ((kiss_fft_scalar)(0.99*SAMP_MAX)) -int curTime; -kiss_fft_scalar adaThresh; -kiss_fft_scalar probEst; - -bool ridgeTracker_isReady; -Array ridgeTracker_out; - -void ridgeTracker_init(void){ - ind = malloc(FRE_LEN*sizeof(kiss_fft_scalar)); - noiseFloor = malloc(FRE_LEN*sizeof(kiss_fft_scalar)); - snrAcc = malloc(FRE_LEN*sizeof(kiss_fft_scalar)); - freAcc = malloc(FRE_LEN*sizeof(kiss_fft_scalar)); - for (size_t f=0; f noiseFloor[f]){ - ind[f] -= 1; - if (ind[f] < 0){ - noiseFloor[f] += S_MUL(noiseFloor[f],NADA); - }else{ - noiseFloor[f] += S_MUL(noiseFloor[f],NADA/2); - } - }else{ - noiseFloor[f] = MAX(EPSILON, noiseFloor[f]-S_MUL(noiseFloor[f],NADA)); - ind[f] = BTLEN; - } - - // snr update - uint32_t r = spec[f]/noiseFloor[f]; - if (r >= 10){ - snr = SAMP_MAX; - }else if (r > 0){ - // see log2fix/main.c - snr = log10fix(r*(1< maxVal){ - maxVal = val; - maxIdx = l; - } - } - snrAcc[f] = maxVal; - freAcc[f] = maxIdx; - } - - // max-pooling - memcpy(snrAccLast, snrAcc, FRE_LEN*sizeof(kiss_fft_scalar)); - for (size_t f=0; f maxVal){ - maxVal = snrAccLast[l]; - maxIdx = l; - } - } - if (maxIdx == f){ - snrAcc[f] = maxVal; - }else{ - snrAcc[f] = 0; - } - } - - // per-bin suppressed stats - memset(snrOut, 0, FRE_LEN*sizeof(kiss_fft_scalar)); - bool binActive = false; - for (size_t f=0; f SUPTHRESH){ - binActive = true; - snrOut[f] = snrAcc[f] - SUPTHRESH; - - if (curTime == -1){ - printf("\nStart timer\n"); - curTime = BTLEN; - } - insertArray(&ridgeTracker_out,snrOut[f],f,(size_t)curTime); - } - } - - // output if applicable - if (curTime != -1){ - if (!binActive){ - printf("Potential output\n"); - // longer time-granuality approximation - if (getMaxDurArray(ridgeTracker_out) >= 50){ - printf("Sufficiently long temporal approximation\n"); - // adaptive thresholding - printf("avgSNR = %d\n", getAvgSNRArray(ridgeTracker_out)); - if (getAvgSNRArray(ridgeTracker_out) >= adaThresh){ - ridgeTracker_isReady = true; - probEst = S_MUL(MUP,probEst) + (SAMP_MAX-MUP); - } else{ - ridgeTracker_reset(); - probEst = S_MUL(MUP,probEst); - } - } else{ - ridgeTracker_reset(); - } - adaThresh = adaThresh + S_MUL(SAMP_MAX-MUT,probEst-PROBTARGET); - - } else{ - curTime += 1; - } - } - - // cleanup - free(snrAccLast); -} - -void ridgeTracker_destroy(void){ - free(ind); - free(noiseFloor); - free(snrAcc); - free(freAcc); - - printf("ridgeTracker_destroy() completed\n"); -} diff --git a/software/apps/audio_module/event_detector/ridgeTracker.h b/software/apps/audio_module/event_detector/ridgeTracker.h deleted file mode 100644 index 114eab1a..00000000 --- a/software/apps/audio_module/event_detector/ridgeTracker.h +++ /dev/null @@ -1,21 +0,0 @@ -// FIXED POINT MATHS -// Long Le -// University of Illinois -// - -#ifndef __RIDGETRACKER_H__ -#define __RIDGETRACKER_H__ - -#include -#include "dynArray.h" - -// The extern vars -extern bool ridgeTracker_isReady; -extern Array ridgeTracker_out; - -void ridgeTracker_init(void); -void ridgeTracker_update(kiss_fft_scalar* spec, kiss_fft_scalar* snrOut); -void ridgeTracker_destroy(void); -void ridgeTracker_reset(void); - -#endif