1 /*
2  * Copyright (C) 2020 Linux Studio Plugins Project <https://lsp-plug.in/>
3  *           (C) 2020 Stefano Tronci <stefano.tronci@protonmail.com>
4  *
5  * This file is part of lsp-plugins
6  * Created on: 14 Apr 2018
7  *
8  * lsp-plugins is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU Lesser General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * any later version.
12  *
13  * lsp-plugins is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  * GNU Lesser General Public License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public License
19  * along with lsp-plugins. If not, see <https://www.gnu.org/licenses/>.
20  */
21 
22 #include <plugins/nonlinear_convolver.h>
23 
24 #define LOAD_EXT        ".lspc"         // Loaded file extension
25 #define TMP_BUF_SIZE    1024
26 #define PRC_BUF_SIZE    (12 * 1024)     // Multiple of 3, 4 and 8
27 #define FFT_MAX_RANK    16
28 
29 namespace lsp
30 {
Loader(nonlinear_convolver_mono * base)31     nonlinear_convolver_mono::Loader::Loader(nonlinear_convolver_mono *base)
32     {
33         pCore = base;
34     }
35 
~Loader()36     nonlinear_convolver_mono::Loader::~Loader()
37     {
38         pCore = NULL;
39     }
40 
run()41     status_t nonlinear_convolver_mono::Loader::run()
42     {
43         pCore->bDataLoaded = false;
44 
45         path_t *path = pCore->pFile->getBuffer<path_t>();
46 
47         if ((path == NULL) || (!path->accepted()))
48         {
49             pCore->nStatus = STATUS_BAD_ARGUMENTS;
50             pCore->pStatus->setValue(pCore->nStatus);
51             return pCore->nStatus;
52         }
53         else
54         {
55             path->commit();
56         }
57 
58         status_t status = pCore->sSyncChirpProcessor.load_from_lspc(path->get_path());
59 
60 //        if (status != STATUS_OK)
61 //            return status;
62 //
63 //        status = pCore->sSyncChirpProcessor.postprocess_nonlinear_convolution(
64 //                    pCore->nModelOrder,
65 //                    pCore->calculate_rank(pCore->nWindowSize),
66 //                    pCore->nWindowSize / 2
67 //                    );
68 //
69         if (status == STATUS_OK)
70             pCore->bDataLoaded = true;
71 
72         return status;
73     }
74 
Preparator(nonlinear_convolver_mono * base)75     nonlinear_convolver_mono::Preparator::Preparator(nonlinear_convolver_mono *base)
76     {
77         pCore = base;
78     }
79 
~Preparator()80     nonlinear_convolver_mono::Preparator::~Preparator()
81     {
82         pCore = NULL;
83     }
84 
run()85     status_t nonlinear_convolver_mono::Preparator::run()
86     {
87         pCore->bDSP_Valid = false;
88 
89         if (!pCore->bDataLoaded)
90             return STATUS_NO_DATA;
91 
92         status_t status = STATUS_OK;
93 
94 //        status = pCore->sSyncChirpProcessor.postprocess_nonlinear_convolution(
95 //                            pCore->nModelOrder,
96 //                            false,
97 //                            10,
98 //                            10,
99 //                            windows::HANN,
100 //                            pCore->calculate_rank(pCore->nWindowSize)
101 //                            );
102 
103         if (status != STATUS_OK)
104             return status;
105 
106         // Handle allocation.
107         if (pCore->bReAllocate || pCore->bIsFirstAllocation)
108         {
109             // DSP Data
110             free_aligned(pCore->pDSPData);
111             pCore->mDSP = NULL;
112 
113             // Convolvers
114             for (size_t n = 0; n < pCore->misoFIRs.nBranches; ++n)
115             {
116                 if (pCore->misoFIRs.FIRConvolvers[n] != NULL)
117                 {
118                     pCore->misoFIRs.FIRConvolvers[n]->destroy();
119                     delete pCore->misoFIRs.FIRConvolvers[n];
120                     pCore->misoFIRs.FIRConvolvers[n] = NULL;
121                 }
122             }
123             delete [] pCore->misoFIRs.FIRConvolvers;
124             pCore->misoFIRs.nBranches   = 0;
125             pCore->misoFIRs.nTaps       = 0;
126 
127             // Set up oversampler
128             switch (pCore->nModelOrder)
129             {
130                 case 2:
131                 {
132                     pCore->sOverPrepare.set_mode(OM_LANCZOS_2X2);
133                     pCore->sOverProcess.set_mode(OM_LANCZOS_2X2);
134                 }
135                 break;
136                 case 3:
137                 {
138                     pCore->sOverPrepare.set_mode(OM_LANCZOS_3X2);
139                     pCore->sOverProcess.set_mode(OM_LANCZOS_3X2);
140                 }
141                 break;
142                 case 4:
143                 {
144                     pCore->sOverPrepare.set_mode(OM_LANCZOS_4X2);
145                     pCore->sOverProcess.set_mode(OM_LANCZOS_4X2);
146                 }
147                 break;
148                 case 6:
149                 {
150                     pCore->sOverPrepare.set_mode(OM_LANCZOS_6X2);
151                     pCore->sOverProcess.set_mode(OM_LANCZOS_6X2);
152                 }
153                 break;
154                 case 8:
155                 {
156                     pCore->sOverPrepare.set_mode(OM_LANCZOS_8X2);
157                     pCore->sOverProcess.set_mode(OM_LANCZOS_8X2);
158                 }
159                 break;
160             }
161 
162             switch (pCore->nDSP)
163             {
164                 case HAMMERSTEIN_FIR:
165                 {
166                     // DSP Data as order by window size matrix of FIR taps, but
167                     // oversampled by order
168                     size_t samples  = pCore->nModelOrder * pCore->nModelOrder * pCore->nWindowSize;
169 
170                     float *ptr      = alloc_aligned<float>(pCore->pDSPData, samples);
171                     if (ptr == NULL)
172                         return STATUS_NO_MEM;
173 
174                     lsp_guard_assert(float *save = ptr);
175                     pCore->mDSP     = ptr;
176                     ptr            += samples;
177 
178                     lsp_assert(ptr <= &save[samples]);
179 
180                     // Create the required convolvers.
181                     pCore->misoFIRs.nBranches   = pCore->nModelOrder;
182                     pCore->misoFIRs.nTaps       = pCore->nModelOrder * pCore->nWindowSize;
183 
184                     pCore->misoFIRs.FIRConvolvers = new Convolver*[pCore->misoFIRs.nBranches]();
185                     Convolver *tmpPtr;
186 
187                     for (size_t n = 0; n < pCore->misoFIRs.nBranches; ++n)
188                     {
189                         tmpPtr = new Convolver();
190 
191                         if (tmpPtr == NULL)
192                         {
193                             return STATUS_NO_MEM;
194                         }
195                         else
196                         {
197                             pCore->misoFIRs.FIRConvolvers[n] = tmpPtr;
198                         }
199                     }
200                 }
201                 break;
202 
203                 default:
204                     return STATUS_NOT_IMPLEMENTED;
205             }
206 
207             pCore->bIsFirstAllocation = false;
208 
209             pCore->nModelOrder_Previous = pCore->nModelOrder;
210             pCore->nWindowSize_Previous = pCore->nWindowSize;
211         }
212 
213         if (status != STATUS_OK)
214             return status;
215 
216         // Assign all DSP structures.
217         switch (pCore->nDSP)
218         {
219             case HAMMERSTEIN_FIR:
220             {
221                 // Randomize phase of the convolver
222                 uint32_t phase  = seed_addr(this);
223                 phase           = ((phase << 16) | (phase >> 16)) & 0x7fffffff;
224                 uint32_t step   = 0x80000000 / (pCore->misoFIRs.nBranches + 1);
225 
226                 for (size_t n = 0; n < pCore->misoFIRs.nBranches; ++n)
227                 {
228                     status = pCore->sSyncChirpProcessor.get_kernel_fir(
229                             &pCore->mDSP[n * pCore->misoFIRs.nTaps], n + 1
230                             );
231 
232                     if (status != STATUS_OK)
233                         return status;
234 
235                     pCore->sOverPrepare.upsample(
236                             &pCore->mDSP[n * pCore->misoFIRs.nTaps],
237                             &pCore->mDSP[n * pCore->misoFIRs.nTaps],
238                             pCore->nWindowSize
239                             );
240 
241                     if(
242                             !pCore->misoFIRs.FIRConvolvers[n]->init(
243                                    &pCore->mDSP[n * pCore->misoFIRs.nTaps],
244                                     pCore->misoFIRs.nTaps,
245                                     FFT_MAX_RANK,
246                                     float((phase + n * step) & 0x7fffffff) / float(0x80000000))
247                             )
248                         return STATUS_NO_MEM;
249                 }
250             }
251             break;
252 
253             default:
254                 return STATUS_NOT_IMPLEMENTED;
255 
256         }
257 
258         if (status == STATUS_OK)
259             pCore->bDSP_Valid = true;
260 
261         return status;
262     }
263 
nonlinear_convolver_mono()264     nonlinear_convolver_mono::nonlinear_convolver_mono(): plugin_t(metadata)
265     {
266         nState                  = IDLE;
267         nDSP                    = HAMMERSTEIN_FIR;
268 
269         pExecutor               = NULL;
270         pLoader                 = NULL;
271         pPreparator             = NULL;
272 
273         nSampleRate             = 0;
274         nStatus                 = STATUS_UNSPECIFIED;
275         fOutGain                = 0;
276         nModelOrder             = 0;
277         nModelOrder_Previous    = 0;
278         nWindowSize             = 0;
279         nWindowSize_Previous    = 0;
280 
281         misoFIRs.FIRConvolvers  = NULL;
282         misoFIRs.nBranches      = 0;
283         misoFIRs.nTaps          = 0;
284 
285         bBypass                 = true;
286         bIsFirstAllocation      = true;
287         bReAllocate             = true;
288         bDataLoaded             = false;
289         bDSP_Valid              = false;
290         bDSP_Prepare_Trigger    = false;
291         bSwitch2Loading         = false;
292         bSwitch2Prepare         = false;
293 
294         mDSP                    = NULL;
295         pDSPData                = NULL;
296 
297         vBuffer                 = NULL;
298         vProcessBufIn           = NULL;
299         vProcessBufTmp          = NULL;
300         vProcessBufOut          = NULL;
301         pData                   = NULL;
302 
303         pIn                     = NULL;
304         pOut                    = NULL;
305 
306         pBypass                 = NULL;
307 
308         pFile                   = NULL;
309         pStatus                 = NULL;
310         pOutGain                = NULL;
311         pOrder                  = NULL;
312         pWindowSize             = NULL;
313         pDSP_Prepare_Trigger    = NULL;
314         pKernelsMesh            = NULL;
315     }
316 
~nonlinear_convolver_mono()317     nonlinear_convolver_mono::~nonlinear_convolver_mono()
318     {
319 
320     }
321 
destroy()322     void nonlinear_convolver_mono::destroy()
323     {
324         if (pLoader != NULL)
325         {
326             delete pLoader;
327             pLoader     = NULL;
328         }
329 
330         if (pPreparator !=NULL)
331         {
332             delete pPreparator;
333             pPreparator = NULL;
334         }
335 
336         free_aligned(pDSPData);
337         mDSP            = NULL;
338 
339         free_aligned(pData);
340         vBuffer         = NULL;
341         vProcessBufIn   = NULL;
342         vProcessBufTmp  = NULL;
343         vProcessBufOut  = NULL;
344     }
345 
get_model_order(size_t order)346     size_t nonlinear_convolver_mono::get_model_order(size_t order)
347     {
348         switch (order)
349         {
350             case nonlinear_convolver_mono_metadata::MODEL_ORDER_2:
351                 return 2;
352             case nonlinear_convolver_mono_metadata::MODEL_ORDER_3:
353                 return 3;
354             case nonlinear_convolver_mono_metadata::MODEL_ORDER_4:
355                 return 4;
356             case nonlinear_convolver_mono_metadata::MODEL_ORDER_6:
357                 return 6;
358             case nonlinear_convolver_mono_metadata::MODEL_ORDER_8:
359                 return 8;
360             default:
361                 return 0;
362         }
363     }
364 
get_window_size(size_t windowSize)365     size_t nonlinear_convolver_mono::get_window_size(size_t windowSize)
366     {
367         switch (windowSize)
368         {
369             case nonlinear_convolver_mono_metadata::WSIZE_ORDER_512:
370                 return 512;
371             case nonlinear_convolver_mono_metadata::WSIZE_ORDER_1024:
372                 return 1024;
373             case nonlinear_convolver_mono_metadata::WSIZE_ORDER_2048:
374                 return 2048;
375             case nonlinear_convolver_mono_metadata::WSIZE_ORDER_4096:
376                 return 4096;
377             case nonlinear_convolver_mono_metadata::WSIZE_ORDER_8192:
378                 return 8192;
379             case nonlinear_convolver_mono_metadata::WSIZE_ORDER_16384:
380                 return 16384;
381             case nonlinear_convolver_mono_metadata::WSIZE_ORDER_32768:
382                 return 32768;
383             case nonlinear_convolver_mono_metadata::WSIZE_ORDER_65536:
384                 return 65536;
385             default:
386                 return 0;
387         }
388     }
389 
init(IWrapper * wrapper)390     void nonlinear_convolver_mono::init(IWrapper *wrapper)
391     {
392         plugin_t::init(wrapper);
393 
394         pExecutor               = wrapper->get_executor();
395 
396         pLoader                 = new Loader(this);
397         pPreparator             = new Preparator(this);
398 
399         sSyncChirpProcessor.init();
400         sOverPrepare.init();
401         sOverProcess.init();
402 
403         // 1X temporary DSP buffer + 3X main process buffers
404         size_t samples          = TMP_BUF_SIZE + 3 * PRC_BUF_SIZE;
405 
406         float *ptr              = alloc_aligned<float>(pData, samples);
407         if (ptr == NULL)
408             return;
409 
410         lsp_guard_assert(float *save = ptr);
411         vBuffer                 = ptr;
412         ptr                    += TMP_BUF_SIZE;
413         vProcessBufIn           = ptr;
414         ptr                    += PRC_BUF_SIZE;
415         vProcessBufTmp          = ptr;
416         ptr                    += PRC_BUF_SIZE;
417         vProcessBufOut          = ptr;
418         ptr                    += PRC_BUF_SIZE;
419 
420         lsp_assert(ptr <= &save[samples]);
421 
422         size_t port_id          = 0;
423 
424         pIn                     = vPorts[port_id++];
425         pOut                    = vPorts[port_id++];
426 
427         pBypass                 = vPorts[port_id++];
428 
429         pFile                   = vPorts[port_id++];
430         pStatus                 = vPorts[port_id++];
431         pOutGain                = vPorts[port_id++];
432         pOrder                  = vPorts[port_id++];
433         pWindowSize             = vPorts[port_id++];
434         pDSP_Prepare_Trigger    = vPorts[port_id++];
435         pKernelsMesh            = vPorts[port_id++];
436     }
437 
update_sample_rate(long sr)438     void nonlinear_convolver_mono::update_sample_rate(long sr)
439     {
440         nSampleRate = sr;
441 
442         sBypass.init(sr);
443         sSyncChirpProcessor.set_sample_rate(sr);
444         sOverPrepare.set_sample_rate(sr);
445         sOverProcess.set_sample_rate(sr);
446     }
447 
fastIntPow(float base,size_t exponent)448     float nonlinear_convolver_mono::fastIntPow(float base, size_t exponent)
449     {
450         return powf(base, exponent);
451 //        if( exponent == 0)
452 //           return 1;
453 //
454 //        float temp = fastIntPow(base, exponent / 2);
455 //
456 //        if (exponent % 2 == 0)
457 //        {
458 //            return temp * temp;
459 //        }
460 //        else
461 //        {
462 //            if(exponent > 0)
463 //                return base * temp * temp;
464 //            else
465 //                return (temp * temp) / base;
466 //        }
467     }
468 
apply_fastIntPow(float * dst,float * src,size_t exponent,size_t count)469     void nonlinear_convolver_mono::apply_fastIntPow(float *dst, float *src, size_t exponent, size_t count)
470     {
471         for (size_t n = 0; n < count; ++n)
472         {
473             dst[n] = fastIntPow(src[n], exponent);
474         }
475     }
476 
calculate_rank(size_t taps)477     size_t nonlinear_convolver_mono::calculate_rank(size_t taps)
478     {
479         size_t rank         = 0;
480         size_t rankSamples  = 1;
481 
482         while (rankSamples < taps)
483         {
484             rankSamples <<= 1;
485             rank += 1;
486         }
487 
488         return rank;
489     }
490 
process_hammerstein_fir(float * dst,float * src,size_t count)491     void nonlinear_convolver_mono::process_hammerstein_fir(float *dst, float *src, size_t count)
492     {
493         size_t nOversampling    = sOverProcess.get_oversampling();
494         size_t max_to_do        = PRC_BUF_SIZE / sOverProcess.get_oversampling();
495 
496         while (count > 0)
497         {
498             size_t to_do = (count < max_to_do) ? count : max_to_do;
499 
500             sOverProcess.upsample(vProcessBufIn, src, to_do);
501 
502             dsp::fill_zero(vProcessBufOut, to_do * nOversampling);
503 
504             for (size_t b = 1; b <= misoFIRs.nBranches; ++b)
505             {
506                 apply_fastIntPow(vProcessBufTmp, vProcessBufIn, b, to_do * nOversampling);
507                 misoFIRs.FIRConvolvers[b - 1]->process(vProcessBufTmp, vProcessBufTmp, to_do * nOversampling);
508                 dsp::add2(vProcessBufOut, vProcessBufTmp, to_do * nOversampling);
509             }
510 
511             sOverProcess.downsample(dst, vProcessBufOut, to_do);
512 
513             count   -= to_do;
514             src     += to_do;
515             dst     += to_do;
516         }
517     }
518 
process(size_t samples)519     void nonlinear_convolver_mono::process(size_t samples)
520     {
521         float *in = pIn->getBuffer<float>();
522         if (in == NULL)
523             return;
524 
525         float *out = pOut->getBuffer<float>();
526         if (out == NULL)
527             return;
528 
529         if (bSwitch2Loading)
530         {
531             pLoader->reset();
532             pPreparator->reset();
533 
534             nState = LOADING;
535 
536             bSwitch2Loading = false;
537         }
538 
539         if (bSwitch2Prepare)
540         {
541             pLoader->reset();
542             pPreparator->reset();
543 
544             nState = PREPARE;
545 
546             bSwitch2Prepare = false;
547         }
548 
549         while (samples > 0)
550         {
551             size_t to_do = (samples > TMP_BUF_SIZE) ? TMP_BUF_SIZE : samples;
552 
553             switch (nState)
554             {
555                 case LOADING:
556                 {
557                     if (pLoader->idle())
558                         pExecutor->submit(pLoader);
559 
560                     if (pLoader->completed())
561                     {
562                         if (pLoader->successful())
563                             nState = PREPARE;
564                         else
565                             nState = IDLE;
566 
567                         pLoader->reset();
568                     }
569 
570                     dsp::fill_zero(vBuffer, to_do);
571                 }
572                 break;
573 
574                 case PREPARE:
575                 {
576                     if (pPreparator->idle())
577                         pExecutor->submit(pPreparator);
578 
579                     if (pPreparator->completed())
580                     {
581                         if (pPreparator->successful())
582                             nState = PROCESSING;
583                         else
584                             nState = IDLE;
585 
586                         pPreparator->reset();
587                     }
588 
589                     dsp::fill_zero(vBuffer, to_do);
590                 }
591                 break;
592 
593                 case PROCESSING:
594                 {
595                     switch (nDSP)
596                     {
597                         case HAMMERSTEIN_FIR:
598                             // Replace with cool function:
599                             process_hammerstein_fir(vBuffer, in, to_do);
600                             break;
601                         case WIENER_HAMMERSTEIN_FIR:
602                             // Replace with cool function:
603                             dsp::fill_zero(vBuffer, to_do);
604                             break;
605                         case HAMMERSTEIN_IIR_BIQUADS:
606                             // Replace with cool function:
607                             dsp::fill_zero(vBuffer, to_do);
608                             break;
609                         case WIENER_HAMMERSTEIN_IIR_BIQUADS:
610                             // Replace with cool function:
611                             dsp::fill_zero(vBuffer, to_do);
612                             break;
613                         default:
614                             dsp::fill_zero(vBuffer, to_do);
615                     }
616 
617                 }
618                 break;
619 
620                 case IDLE:
621                 default:
622                 {
623                     dsp::fill_zero(vBuffer, to_do);
624                 }
625             }
626 
627             dsp::mul_k2(vBuffer, fOutGain, to_do);
628 
629             sBypass.process(out, in, vBuffer, to_do);
630 
631             in         += to_do;
632             out        += to_do;
633             samples    -= to_do;
634         }
635     }
636 
update_settings()637     void nonlinear_convolver_mono::update_settings()
638     {
639         bBypass                 = pBypass->getValue() >= 0.5f;
640         sBypass.set_bypass(bBypass);
641 
642         // File extension check
643         path_t *path            = pFile->getBuffer<path_t>();
644 
645         if ((path != NULL) && (path->pending()))
646         {
647             const char *fname   = path->get_path();
648             size_t len          = strlen(fname);
649 
650             const char *ext     = LOAD_EXT;
651             size_t extlen       = strlen(ext);
652 
653             if (len < extlen)
654                 nStatus         = (len > 0) ? STATUS_BAD_ARGUMENTS : STATUS_UNSPECIFIED;
655             else
656             {
657                 nStatus         = STATUS_OK;
658 
659                 for (size_t n   = 0; n < extlen; ++n)
660                 {
661                     if (fname[len - extlen + n] != ext[n]) // This is case sensitive
662                     {
663                         nStatus = STATUS_BAD_ARGUMENTS;
664                         break;
665                     }
666                 }
667 
668             }
669 
670             path->accept();
671             bSwitch2Loading     = true;
672         }
673 
674         pStatus->setValue(nStatus);
675 
676         fOutGain                = pOutGain->getValue();
677 
678 //        size_t previousOrder    = nModelOrder;
679         nModelOrder             = get_model_order(pOrder->getValue());
680 
681 //        size_t previousWSize    = nWindowSize;
682         nWindowSize             = get_window_size(pWindowSize->getValue());
683 
684         bReAllocate             =
685                                   (nModelOrder_Previous != nModelOrder) ||
686                                   (nWindowSize_Previous != nWindowSize);
687 
688         bool previousTrigger    = bDSP_Prepare_Trigger;
689         bDSP_Prepare_Trigger    = pDSP_Prepare_Trigger->getValue() >= 0.5f;
690 
691         // Was the trigger pressed (are we passing from false to true?). Also,
692         // was a new file loaded? In the case, we cannot jump to prepare, but
693         // we will go to loading first.
694         if (!previousTrigger && bDSP_Prepare_Trigger && !bSwitch2Loading)
695             bSwitch2Prepare     = true;
696         else
697             bSwitch2Prepare     = false;
698     }
699 }
700