From 24aed96f59ed4043aac008af81b49c3040168fa9 Mon Sep 17 00:00:00 2001 From: Andrea759 Date: Tue, 12 Jul 2022 00:05:08 +0200 Subject: [PATCH] release and testing --- Jack_AntiLarsen/include/Analyzer.h | 14 ++-- Jack_AntiLarsen/include/DSP.h | 2 +- Jack_AntiLarsen/include/const.h | 10 ++- Jack_AntiLarsen/include/iir.h | 2 +- Jack_AntiLarsen/main.cpp | 15 +++-- Jack_AntiLarsen/src/Analyzer.cpp | 101 +++++++++++++---------------- Jack_AntiLarsen/src/DSP.cpp | 28 ++++---- Jack_AntiLarsen/src/iir.cpp | 10 +-- Jack_AntiLarsen/test.cpp | 90 +++++++++++++++++++++++++ 9 files changed, 181 insertions(+), 91 deletions(-) diff --git a/Jack_AntiLarsen/include/Analyzer.h b/Jack_AntiLarsen/include/Analyzer.h index e39e492..c8ad747 100644 --- a/Jack_AntiLarsen/include/Analyzer.h +++ b/Jack_AntiLarsen/include/Analyzer.h @@ -11,7 +11,7 @@ class Analyzer { public: explicit Analyzer(const bool *settings); - void run(const float *jackBuffer); + void analyzeBuffer(const float *jackBuffer); void setInputBuffer(const float *jackBuffer); void setEnableAlgo(const bool *); float getPhprThreshold() const; @@ -23,17 +23,19 @@ class Analyzer { int found_howls[N_PEAKS]; virtual ~Analyzer(); + std::complex *getFtOut() const; + float *getOutBuffer() const; + private: bool phpr(const float *); bool pnpr(const float *); bool imsd(const float *); void fftWrapper(const float *); - void updateBuffer(const float *); static void inline minHead(const float*, int *); void blackman_win(int); fftwf_plan ft_plan; - std::complex *ft_in; + float *ft_in; std::complex *ft_out; const float *jack_buffer; float *buffers[3]; @@ -41,9 +43,9 @@ class Analyzer { /* * Thresholds for howling frequencies detection (in dB) */ - float phpr_threshold = 10.0f; - float pnpr_threshold = 30.0f; - float imsd_threshold = 1.0f; + static constexpr float phpr_threshold = 10.0f; + static constexpr float pnpr_threshold = 25.5f; + static constexpr float imsd_threshold = 1.0f; bool run_phpr = true; bool run_pnpr = true; diff --git a/Jack_AntiLarsen/include/DSP.h b/Jack_AntiLarsen/include/DSP.h index 8fac590..96de891 100644 --- a/Jack_AntiLarsen/include/DSP.h +++ b/Jack_AntiLarsen/include/DSP.h @@ -23,13 +23,13 @@ class DSP{ bool applyFilters(); void reset_bank(); void add_filter_to_bank(int index,t_filter*); - t_bank bank; void setInOutBuffers(float *in, float* out, int nframes); float *getBufOut() const; private: float *buf_in; float *buf_out; int nframes = 0; + t_bank bank; }; #endif diff --git a/Jack_AntiLarsen/include/const.h b/Jack_AntiLarsen/include/const.h index 89fed12..708734c 100644 --- a/Jack_AntiLarsen/include/const.h +++ b/Jack_AntiLarsen/include/const.h @@ -6,9 +6,15 @@ #define JACK_ANTILARSEN_CONST_H #define BUF_LENGTH 1024 -#define N_PRE_FILTERS 1024 +#define N_PRE_FILTERS 512 #define MAX_ACTIVE_FILTERS 10 #define PI 3.1415926 #define N_PEAKS 10 - +#define FSTEP 46.875 //Hz +#define INT_TIME 21.333 //ms +/* + * r2c fftw output format: + * https://www.fftw.org/fftw3_doc/The-1d-Real_002ddata-DFT.html + */ +#define FFTOUT_BUF_LENGTH BUF_LENGTH/2+1 #endif diff --git a/Jack_AntiLarsen/include/iir.h b/Jack_AntiLarsen/include/iir.h index 205ce88..d259c18 100644 --- a/Jack_AntiLarsen/include/iir.h +++ b/Jack_AntiLarsen/include/iir.h @@ -26,7 +26,7 @@ class preFiltersBank{ void setQFactor(float new_q_factor); private: // Default values - float f_sampling = 44100; + float f_sampling = 48000; float gb = -10; float q_factor = 30; float f_step; diff --git a/Jack_AntiLarsen/main.cpp b/Jack_AntiLarsen/main.cpp index aa832c9..d68c3ee 100644 --- a/Jack_AntiLarsen/main.cpp +++ b/Jack_AntiLarsen/main.cpp @@ -7,6 +7,7 @@ #include "include/DSP.h" #include #include +#include jack_port_t *input_port; jack_port_t *output_port; @@ -21,7 +22,7 @@ int process (jack_nframes_t nframes, void *arg){ bool larsen = false; in = (float *)jack_port_get_buffer (input_port, nframes); out = (float *)jack_port_get_buffer (output_port, nframes); - analyzer->run(in); + analyzer->analyzeBuffer(in); dsp->setInOutBuffers(in,out, (int) nframes); for (auto f_idx: analyzer->found_howls) { @@ -29,7 +30,7 @@ int process (jack_nframes_t nframes, void *arg){ dsp->add_filter_to_bank(f_idx, filters->filters); larsen = true; // debug only - std::cout << f_idx << std::endl; + std::cout << f_idx*FSTEP << std::endl; } } if (larsen) dsp->applyFilters(); @@ -60,13 +61,13 @@ int main (int argc, char *argv[]) bool analyzer_settings[3]; /* Default settings */ - analyzer_settings[0]= true; // enable pnpr - analyzer_settings[1]= true; // enable phpr + analyzer_settings[0]= true; // enable phpr + analyzer_settings[1]= true; // enable pnpr analyzer_settings[2]= false; // enable imsd - float f_sampling = 44100; + float gb = -10; float q_factor = 30; - float f_min = 20; + float f_min = 0; /* open a client connection to the JACK server */ client = jack_client_open (client_name, options, &status, server_name); @@ -186,7 +187,7 @@ int main (int argc, char *argv[]) free (ports); /* keep running until stopped by the user */ - while(true); + sleep(-1); /* this is never reached but if the program * had some other way to exit besides being killed, * they would be important to call. diff --git a/Jack_AntiLarsen/src/Analyzer.cpp b/Jack_AntiLarsen/src/Analyzer.cpp index 47a7467..401243d 100644 --- a/Jack_AntiLarsen/src/Analyzer.cpp +++ b/Jack_AntiLarsen/src/Analyzer.cpp @@ -12,50 +12,47 @@ Analyzer::Analyzer(const bool settings[3]) { // locate jack audio buffer this->jack_buffer = nullptr; - /* - * todo: test fftw_complex *fftw_malloc() for SSE/AVX speedup - */ - this->ft_in = (std::complex *) calloc(BUF_LENGTH, sizeof(std::complex)); - this->ft_out = (std::complex *) calloc(BUF_LENGTH, sizeof(std::complex)); - for (auto &buf: this->buffers) buf = (float *) calloc(BUF_LENGTH, sizeof(float)); + this->ft_in = (float *) calloc(BUF_LENGTH,sizeof(float)); + this->ft_out = (std::complex *) fftwf_alloc_complex(FFTOUT_BUF_LENGTH); + + for (auto &buf: this->buffers) + buf = (float *) calloc(FFTOUT_BUF_LENGTH, sizeof(float)); + for (auto &i: this->found_howls) i = 0; + // select which algorithm will be run setEnableAlgo(settings); blackman_win(BUF_LENGTH); // fftw3 plan creation - this->ft_plan = fftwf_plan_dft_1d(BUF_LENGTH, reinterpret_cast(ft_in),\ - reinterpret_cast(ft_out), FFTW_FORWARD, FFTW_MEASURE); + this->ft_plan = fftwf_plan_dft_r2c_1d(BUF_LENGTH, ft_in,\ + reinterpret_cast(ft_out), FFTW_MEASURE); } void Analyzer::blackman_win(int nsamples) { - int nsamples_periodic = nsamples++; - for (int i = 0; i < (nsamples_periodic+1)/2; ++i) { + int nsamples_periodic = nsamples + 1; + this->blackman[0]=0.0f; + for (int i = 1; i < (nsamples_periodic+1)/2; ++i) { this->blackman[i]=0.42f - 0.5f*cosf(2.0f*PI*(float)i/(nsamples_periodic-1)) + 0.08f*cosf(4.0f*PI*(float)i/(nsamples_periodic-1)); - if (i + (nsamples_periodic+1)/2 < nsamples){ - this->blackman[i+(nsamples_periodic+1/2)]= this->blackman[i]; - } + this->blackman[nsamples_periodic-i-1]= this->blackman[i]; } } /* - * buf_out = abs(fft(buf_in)) + * buf_out = abs(fft(jack_buf)) */ -void Analyzer::fftWrapper(const float *buf_in) { +void Analyzer::fftWrapper(const float *jack_buf) { for (int i = 0; i < BUF_LENGTH; ++i) { - this->ft_in[i] = buf_in[i]*this->blackman[i]; + this->ft_in[i] = jack_buf[i]* this->blackman[i]; } fftwf_execute(this->ft_plan); - for (int i = 0; i < BUF_LENGTH; ++i) { - this->buffers[0][i] = std::abs(this->ft_out[i]); + for (int i = 0; i < (BUF_LENGTH/2)+1; ++i) { + this->buffers[0][i] = std::abs(this->ft_out[i])/213.0f; } } -void Analyzer::updateBuffer(const float *jack_buffer_update) { - this->jack_buffer = jack_buffer_update; -} /* * Peak to Harmonics Power Ratio: * Test if the frequency peaks at index[i] has significant harmonics. @@ -68,14 +65,14 @@ bool Analyzer::phpr(const float *buf_in) { for (int i = 0; i < N_PEAKS; ++i) { if (this->found_howls[i] != 0){ ret = true; - if (i*2 < BUF_LENGTH) { - if (2 * log10f_fast(buf_in[this->found_howls[i]] / buf_in[this->found_howls[i] * 2])\ + if (i*2 < FFTOUT_BUF_LENGTH) { + if (20 * log10f_fast(buf_in[this->found_howls[i]] / buf_in[this->found_howls[i] * 2])\ < phpr_threshold) { this->found_howls[i] = 0; } } - if (i*3 < BUF_LENGTH){ - if (2 * log10f_fast(buf_in[this->found_howls[i]] / buf_in[this->found_howls[i] * 3])\ + if (i*3 < FFTOUT_BUF_LENGTH){ + if (20 * log10f_fast(buf_in[this->found_howls[i]] / buf_in[this->found_howls[i] * 3])\ < phpr_threshold) { this->found_howls[i] = 0; } @@ -96,26 +93,14 @@ bool Analyzer::pnpr(const float *buf_in) { for (int i = 0; i < N_PEAKS; ++i) { if(this->found_howls[i] != 0){ ret = true; - if (i+1 < BUF_LENGTH){ - if (2*log10f_fast(buf_in[this->found_howls[i]] / buf_in[this->found_howls[i]+1])\ - < pnpr_threshold) { - this->found_howls[i] = 0; - } - } - if (i+2 < BUF_LENGTH) { - if (2*log10f_fast(buf_in[this->found_howls[i]] / buf_in[this->found_howls[i] + 2])\ + if (i+1 < FFTOUT_BUF_LENGTH){ + if (20 *log10f_fast(buf_in[this->found_howls[i]] / buf_in[this->found_howls[i]+1])\ < pnpr_threshold) { this->found_howls[i] = 0; } } if (i-1 > 0){ - if (2*log10f_fast(buf_in[this->found_howls[i]] / buf_in[this->found_howls[i] - 1])\ - < pnpr_threshold) { - this->found_howls[i] = 0; - } - } - if (i-2 > 0) { - if (2*log10f_fast(buf_in[this->found_howls[i]] / buf_in[this->found_howls[i] - 2])\ + if (20 *log10f_fast(buf_in[this->found_howls[i]] / buf_in[this->found_howls[i] - 1])\ < pnpr_threshold) { this->found_howls[i] = 0; } @@ -138,28 +123,29 @@ bool Analyzer::imsd(const float *buf_in) { * find 10 highest peaks, run phpr, pnpr, imsd if enabled. * write howling frequencies indexes in found_howls, if detected. */ -void Analyzer::run(const float *jackBuffer) { - updateBuffer(jackBuffer); +void Analyzer::analyzeBuffer(const float *jackBuffer) { + // swap buffers order + float *t = this->buffers[2]; + this->buffers[2]= this->buffers[1]; + this->buffers[1]= this->buffers[0]; + this->buffers[0]=t; + + setInputBuffer(jackBuffer); fftWrapper(this->jack_buffer); float *buf_in = this->buffers[0]; for (int & found_howl : this->found_howls) { found_howl = 0; } - for (int i = 0; i < BUF_LENGTH; ++i) { + for (int i = 0; i < FFTOUT_BUF_LENGTH; ++i) { if (buf_in[i] > buf_in[this->found_howls[0]]){ this->found_howls[0] = i; minHead(buf_in,this->found_howls); } } // recognize larsen effetcs - if(this->run_pnpr) { if (!pnpr(buf_in)) return; } if(this->run_phpr) { if (!phpr(buf_in)) return; } + if(this->run_pnpr) { if (!pnpr(buf_in)) return; } if(this->run_imsd) imsd(buf_in); - // swap buffers order - float *t = this->buffers[2]; - this->buffers[2]= this->buffers[1]; - this->buffers[1]= this->buffers[0]; - this->buffers[0]=t; } void inline Analyzer::minHead(const float *buf_in, int *peaks) { @@ -183,23 +169,21 @@ Analyzer::~Analyzer() { fftwf_destroy_plan(this->ft_plan); free(ft_in); free(ft_out); - free(found_howls); for (auto &buf: this->buffers) { free(buf); } - free(buffers); } float Analyzer::getPhprThreshold() const { - return phpr_threshold; + return this->phpr_threshold; } float Analyzer::getPnprThreshold() const { - return pnpr_threshold; + return this->pnpr_threshold; } float Analyzer::getImsdThreshold() const { - return imsd_threshold; + return this->imsd_threshold; } bool Analyzer::isRunPhpr() const { @@ -215,8 +199,13 @@ bool Analyzer::isRunImsd() const { } void Analyzer::setInputBuffer(const float *jackBuffer) { - jack_buffer = jackBuffer; + this->jack_buffer = jackBuffer; } +std::complex *Analyzer::getFtOut() const { + return ft_out; +} - +float *Analyzer::getOutBuffer() const { + return buffers[0]; +} diff --git a/Jack_AntiLarsen/src/DSP.cpp b/Jack_AntiLarsen/src/DSP.cpp index 42bc9d2..d692037 100644 --- a/Jack_AntiLarsen/src/DSP.cpp +++ b/Jack_AntiLarsen/src/DSP.cpp @@ -18,15 +18,17 @@ DSP::DSP() { } void DSP::add_filter_to_bank(int index, t_filter filters[]) { - if (this->bank.active_filters < MAX_ACTIVE_FILTERS){ - this->bank.p_filter[this->bank.active_filters] = &filters[index]; - this->bank.active_filters++; - } - if (this->bank.active_filters == MAX_ACTIVE_FILTERS){ - this->bank.p_filter[this->bank.next_insert] = &filters[index]; - this->bank.next_insert++; - if (this->bank.next_insert == 9){ - this->bank.next_insert = 0; + if (index < N_PRE_FILTERS){ + if (this->bank.active_filters < MAX_ACTIVE_FILTERS){ + this->bank.p_filter[this->bank.active_filters] = &filters[index]; + this->bank.active_filters++; + } + if (this->bank.active_filters == MAX_ACTIVE_FILTERS){ + this->bank.p_filter[this->bank.next_insert] = &filters[index]; + this->bank.next_insert++; + if (this->bank.next_insert == 9){ + this->bank.next_insert = 0; + } } } } @@ -41,11 +43,11 @@ bool DSP::applyFilters() { if (this->bank.p_filter[i] != nullptr){ for (int j = 0; j < BUF_LENGTH; ++j) { // IIR Filters Equation - buf_out[j] = this->bank.p_filter[i]->e * buf_in[j] + + buf_out[j] = 0.5f*(this->bank.p_filter[i]->e * buf_in[j] + this->bank.p_filter[i]->p * x_1 + - this->bank.p_filter[i]->d[0] * x_2 + - this->bank.p_filter[i]->d[1] * y_1 + - this->bank.p_filter[i]->d[2] * y_2; + this->bank.p_filter[i]->d[0] * x_2 - + this->bank.p_filter[i]->d[1] * y_1 - + this->bank.p_filter[i]->d[2] * y_2); // shift delayed x, y samples x_2 = x_1; x_1 = buf_in[j]; diff --git a/Jack_AntiLarsen/src/iir.cpp b/Jack_AntiLarsen/src/iir.cpp index d98536a..39f2901 100644 --- a/Jack_AntiLarsen/src/iir.cpp +++ b/Jack_AntiLarsen/src/iir.cpp @@ -24,19 +24,19 @@ preFiltersBank::preFiltersBank(float f_sampling_in, float gb_in, float q_factor_ this->f_sampling = f_sampling_in; this->gb = gb_in; this->q_factor = q_factor_in; - this->f_step = (f_sampling_in - f_min_in) / (float)N_PRE_FILTERS; + this->f_step = (f_sampling_in/2 - f_min_in) / (float)N_PRE_FILTERS; this->f_min = f_min_in; f= f_min_in; - damp = sqrtf(std::abs(1.0f - powf(gb_in, 2.0f))) / gb_in; + damp = sqrtf(std::abs(1.0f - powf(gb_in, 2.0f)))/ gb_in; for (auto &filter: this->filters) { - wo = 2.0f * PI * (f) / f_sampling_in; + wo = 2.0f * PI * (f) / (f_sampling_in/2); filter.e = 1.0f/(1.0f + damp*tanf(wo/(q_factor_in * 2.0f))); filter.p = cosf(wo); filter.d[0] = filter.e; - filter.d[1] = 2*filter.e*filter.p; - filter.d[2] = (2*filter.e-1); + filter.d[1] = 2.0f*filter.e*filter.p; + filter.d[2] = (2.0f*filter.e-1.0f); f += f_step; } } diff --git a/Jack_AntiLarsen/test.cpp b/Jack_AntiLarsen/test.cpp index b16ef37..224b92a 100644 --- a/Jack_AntiLarsen/test.cpp +++ b/Jack_AntiLarsen/test.cpp @@ -7,14 +7,104 @@ #include "include/iir.h" #include "include/Analyzer.h" #include "include/DSP.h" +#include "include/const.h" #include #include +#include +#include preFiltersBank *filters; DSP *dsp; Analyzer *analyzer; int main(){ + float out_buf_test[BUF_LENGTH]; + int ntests = 3; + bool settings[3]; + int test_peaks[10]; + settings[0] = true; + settings[1] = true; + settings[2] = false; + bool larsen; + const static float fsampling = 48000.0f; + filters = new preFiltersBank(fsampling,10,30,0); + dsp = new DSP; + analyzer = new Analyzer(settings); + float test[ntests][BUF_LENGTH] ; + for (int i = 0; i < BUF_LENGTH; ++i) { + test[0][i] = 0.9f * cosf((float)i*2.0f*PI*42.0f*FSTEP/fsampling); + } + for (int i = 0; i < BUF_LENGTH; ++i) { + test[1][i] = 0.9f * cosf((float)i*2.0f*PI*2.0f*FSTEP/fsampling); + } + for (int i = 0; i < BUF_LENGTH; ++i) { + test[2][i] = 0.9f * cosf((float)i*2.0f*PI*8000.f/fsampling); + } + if (analyzer->isRunPhpr()) std::cout << "phpr "; + if (analyzer->isRunPnpr()) std::cout << "pnpr "; + if (analyzer->isRunImsd()) std::cout << "imsd "; + std::cout << std::endl; + for (int i = 0; i < ntests; ++i) { + analyzer->analyzeBuffer(test[i]); + std::cout << "test " << i << ": "; + for (auto &n: analyzer->found_howls) { + std::cout << n << " "; + } + memcpy(test_peaks,analyzer->found_howls, 10*sizeof(int)); + std::cout << std::endl; + std::cout << "peaks: "; + for (auto &n: analyzer->found_howls) { + if (n != 0)std::cout << (analyzer->getOutBuffer()[n]) << " "; + } + std::cout << std::endl; + std::cout << "filtered: "; + analyzer->analyzeBuffer(test[i]); + dsp->setInOutBuffers(test[i],out_buf_test,BUF_LENGTH); + for (auto f_idx: analyzer->found_howls) { + if (f_idx != 0){ + dsp->add_filter_to_bank(f_idx, filters->filters); + larsen = true; + } + } + if (larsen) dsp->applyFilters(); + dsp->reset_bank(); + analyzer->analyzeBuffer(out_buf_test); + larsen = false; + for (auto &n: test_peaks) { + if (n != 0)std::cout << analyzer->getOutBuffer()[n] << " "; + } + std::cout << std::endl; + } + + + + + + + + // Timing test + for (int i = 0; i < ntests; ++i) { + auto begin = std::chrono::high_resolution_clock::now(); + analyzer->analyzeBuffer(test[i]); + dsp->setInOutBuffers(test[i],out_buf_test,BUF_LENGTH); + for (auto f_idx: analyzer->found_howls) { + if (f_idx != 0){ + dsp->add_filter_to_bank(f_idx*2, filters->filters); + larsen = true; + } + } + if (larsen) dsp->applyFilters(); + larsen = false; + dsp->reset_bank(); + auto end = std::chrono::high_resolution_clock::now(); + std::cout << "test " << i << " time: "<< \ + std::chrono::duration_cast(end-begin).count() \ + << "ns" << std::endl; + } + + delete analyzer; + delete dsp; + delete filters; return 0; } \ No newline at end of file