This paper focuses on the use of GPGPU (General-Purpose computing on Graphics Processing Units) for audio processing. This is a promising approach to problems where a high parallelization of tasks is desirable. Within the context of binaural spatialization we will develop a convolution engine having in mind both offline and real-time scenarios, and the support for multiple sound sources. Details on implementations and strategies used with both dominant technologies, namely CUDA and OpenCL, will be presented highlighting both advantages and issues. Comparisons between this approach and typical CPU implementations will be presented as well as between frequency (FFT) and time-domain approaches. Results will show that benefits exist in terms of execution time for a number of situations

Mauro, Davide Andrea

English

Marshall University

Marshall UniversityMarshall Digital ScholarWeisberg Division of Computer Science FacultyResearch Weisberg Division of Computer ScienceSpring 4-22-2012Audio convolution by the mean of GPU: CUDAand OpenCL implementationsDavide Andrea MauroMarshall University, maurod@marshall.eduFollow this and additional works at: http://mds.marshall.edu/wdcs_facultyPart of the Computer Sciences CommonsThis Article is brought to you for free and open access by the Weisberg Division of Computer Science at Marshall Digital Scholar. It has been acceptedfor inclusion in Weisberg Division of Computer Science Faculty Research by an authorized administrator of Marshall Digital Scholar. For moreinformation, please contact zhangj@marshall.edu, martj@marshall.edu.Recommended CitationMauro D.A. Audio convolution by the mean of GPU: CUDA and OpenCL implementations. Proceedings of the Acoustics 2012Nantes Conference, 23-27 April 2012, Nantes, France: pp.2863-2868, 2012.Audio convolution by the mean of GPU: CUDA andOpenCL implementationsD.A. MauroLaboratorio di Informatica Musicale (LIM), Dipartimento di Informatica e Comunicazione(DICo), Universita` degli Studi di Milano, Via Comelico 39/41, 20135 Milano, Italymauro@dico.unimi.itProceedings of the Acoustics 2012 Nantes Conference 23-27 April 2012, Nantes, France2857This paper focuses on the use of GPGPU (General-Purpose computing on Graphics Processing Units) for audioprocessing. This is a promising approach to problems where a high parallelization of tasks is desirable. Withinthe context of binaural spatialization we will develop a convolution engine having in mind both offline and real-time scenarios, and the support for multiple sound sources. Details on implementations and strategies used withboth dominant technologies, namely CUDA and OpenCL, will be presented highlighting both advantages andissues. Comparisons between this approach and typical CPU implementations will be presented as well as betweenfrequency (FFT) and time-domain approaches. Results will show that benefits exist in terms of execution time fora number of situations.Figure 1: The workflow diagram of the system.1 IntroductionWe first introduce the core of the work in terms of con-ceptualization and development of a model. Even if the pro-cess is well known and understood in terms of mathematics,the realization of implementations that work in real-life sce-narios is not trivial. One of the greatest obstacle is the com-putational complexity that convolution requires both in thetime and frequency domain approaches. This means that theproblem could be theoretically solved but the computer ar-chitecture does not allow it to be solved in a reasonable timefor some practical cases of interest.2 Convolution EnginesAs shown in Figure 1 the system requires as input an ane-choic signal (monophonic) and a impulse response (stereo)and the overall output will be two channel spatialized soundthat can feed both headphones or loudspeakers (with crosstalkcancelation algorithms [2]).We will focus on implementations of this system thanksto modern GPGPU techniques.2.1 State of the ArtIn the literature there are other systems that aim at re-alizing systems that achieve real-time auralization, or aug-mented reality. We present a brief sketch of the opportunitiesand the techniques employed. It is worth to cite the work ofKapralos et al. presented in [4] and [5] where the authors ap-ply GPGPU techniques to solve the problem of convolution.The main differences are that the authors use a time domainimplementation that exploits the use of OpenGL in order toprocess audio data. This means basically that they need totweak the system to threat audio data as RGB bitmaps.• TConvolutionUB∼: A Max/MSP external patch fromThomas Resch that extends the possibilities given bythe buffir∼ object allowing convolution with a filterthat has more than 255 points.• SIR2: An easy to use native audio-plugin to use forhigh quality reverberation. It’s available for the pluginformats VST and AudioUnit. Its use can be stretchedfrom a convolution reverb to a convolution engine forauralization given the flexibility of the program itself.• djbfft: A library for floating-point convolution. Thecurrent version provides power-of-2 complex FFTs, realFFTs at twice the speed, and fast multiplication of com-plex arrays. Single precision and double precision areequally supported.• BruteFIR: An open-source convolution engine, a pro-gram for applying long FIR filters to multi-channeldigital audio, either offline or in realtime, by AndersTorger [8]. Its basic operation is specified througha configuration file, and filters, attenuation and delaycan be changed at runtime through a simple commandline interface. The author states that the FIR filter al-gorithm used is an optimized frequency domain algo-rithm, partly implemented in hand-coded assembler,thus throughput is extremely high. In real-time, a stan-dard computer can typically run more than 10 channelswith more than 60000 filter taps each. It makes use ofthe partitioned convolution and overlap-save methodsthat are introduced in the following subsection.• AlmusVCU: From the author of BruteFIR this is a com-plete system that aims at an integrated environment forsound spatialization. It has been designed primarilywith Ambiophonics in mind and contains all process-ing needed for a complete Ambiophonics system.• Aurora Plugin: From Angelo Farina, is a suite of plug-ins for Adobe Audition: room acoustical impulse re-sponses can be measured and manipulated, for the recre-ation of audible, three-dimensional simulations of theacoustical space.2.2 Convolution in the Time DomainThis approach can be mathematically described by theformula:y(k) =∑j=1x1( j)x2(k − j + 1) (1)Proceedings of the Acoustics 2012 Nantes Conference23-27 April 2012, Nantes, France2858Where x1 and x2 are the input sequences of length m and nand y is the output sequence of length k = m + n − 1.When m = n, which is the normal case for other imple-mentations, this gives:w(1) = u(1)v(1)w(2) = u(1)v(2) + u(2)v(1)w(3) = u(1)v(3) + u(2)v(2) + u(3)v(1)· · ·w(n) = u(1)v(n) + u(2)v(n − 1) + · · · + u(n)v(1)· · ·w(2n − 1) = u(n)v(n)(2)The computational complexity for the time domain approachis O(n2).This is the underlying approach to every other method.Implementing a FIR (Finite Impulse Response) filter is obvi-ously the easiest idea but as can be seen from the complexityas the input size increase it could become impossible to pro-cess data in real-time.2.3 Convolution in the Frequency DomainThanks to the convolution theorem we can express theconvolution of two sequences as the multiplication of theirFourier transforms. Here the general layout for the frequencydomain approach is introduced. The approach that can beschematized as follows:• Zero-Pad input vectors x1 and x2 of length m and n sothe length of the sequences becomes m + n − 1;• Perform FFT of the input vectors;• Perform the pointwise multiplication of the two se-quences;• Perform the IFFT of the obtained sequence.The computational complexity for the frequency domain ap-proach is O(n log(n)).2.3.1 Overlap-add algorithmSince the size of the input can become very high, it isnot convenient to use a single window to transform the en-tire signal so a number of methods can be implemented toovercome this. We choose to use a method called Overlap-add (OA, OLA). It is an efficient way to evaluate the discreteconvolution of a very long signal x[n] with a finite impulseresponse (FIR) filter h[n]. The concept is to divide the prob-lem into multiple convolutions of h[n] with short segmentsof x[n]:y[n] = x[n] ∗ h[n] :=∞∑m=−∞h[m]x[n − m] =M∑m=1h[m]x[n − m](3)where h[m] = 0 for m outside the region [1,M].xk[n] :=⎧⎪⎪⎨⎪⎪⎩x[n + kL] n = 1, 2, ···, L0 otherwise(4)Figure 3: Schematic view of the overlap-add convolutionmethod.where L is an arbitrary segment length.x[n] =∑kxk[n − kL] (5)So y[n] can be written as a sum of convolutions:y[n] =⎛⎜⎜⎜⎜⎜⎝∑kxk[n − kL]⎞⎟⎟⎟⎟⎟⎠ ∗ h[n] =∑k(xk[n − kL] ∗ h[n]) (6)The method is depicted in Figure 3It is particularly useful for our tasks since it works onindependent pieces of input and thus is well suited for a par-allelized approach such as one that employs a GPU.3 Reference CPU implementationsIn order to make comparisons with the GPU implemen-tations that we will present we need a reference implementa-tion that can serve as a basis in terms of execution time andbitwise precision. For this reason three different prototypeshave been developed that use different algorithms.The first two prototypes are Matlab scripts that use both aTime Domain and a Frequency Domain approach. Since thecomputational complexity for the Time Domain approach isO(n2) this can not be used when the filter kernels are big. Inour experiments, according to a Max/MSP implementationthat will be introduced in the following section, we choose tolimit the size to 256 samples.The frequency domain implementation (presented in [7])will be used to validate the results in terms of bitwise preci-sion. Since Matlab is mainly intended as a prototyping en-vironment there is no focus on performance and every otherimplementation can outperform our Matlab testbase by or-ders of magnitude. Moreover, this implementation worksonly in “direct mode”; this implies that a single FFT is per-formed for the entire signal and therefore the algorithm mayProceedings of the Acoustics 2012 Nantes Conference 23-27 April 2012, Nantes, France2859Figure 2: A scheme of convolution in frequency domain.not be applicable for long sequences due to memory con-straints or implementation limits. Source code for both theMatlab implementations are available from the author.The last CPU implementation is written in C++ and isbased on the FFTW3 library (see [6]). It is based on the ar-chitecture presented in Figure 2 and implements both modal-ities (Direct and OLA) previously discussed.The FFTW library itself is based on Cooley-Tukey al-gorithm [3]. As presented by the authors, the interactionof the user with FFTW occurs in two stages: planning, inwhich FFTW adapts to the hardware, and execution, in whichFFTW performs useful work for the user. To compute a DFT,the user first invokes the FFTW planner, specifying the prob-lem to be solved. The problem is a data structure that de-scribes the “shape” of the input data - array sizes and mem-ory layouts - but does not contain the data itself. In return,the planner yields a plan, an executable data structure thataccepts the input data and computes the desired DFT. After-wards, the user can execute the plan as many times as desired.3.1 A CUDA convolution engineFor the CPU implementation with CUDA we were ableto implement both Direct and OLA algorithm. We considerthe benefits of both approaches in the following section whilepresenting performance comparisons. For FFT we use a li-brary called CUFFT which is actually based on FFTW3 li-brary with some other optimizations specifically designed forGPUs. One of the current issue is the CUFFT limit of 64 mil-lions of points.3.2 An OpenCL convolution engineOne of the current limitations is that the factorization al-gorithms works only for powers of 2 (radix-2). So the pay-load should be adapted to make the sum with the length ofthe filter kernel to be the closest greater power of 2.4 The CGPUconv prototypeFrom a number of the previously cited prototypes wederived a single application that allows the user to choosebetween a CPU- or a GPU-based algorithm and between adirect mode (a single window for the entire signal) and anOverlap-add mode. It is structured as a “wrapper” aroundthe single module that has the capability of opening audiofiles and writing them back to disk thanks to libsndfile (see[1]). It is a command line tool that compiles and executesboth on Microsoft Windows, Apple Mac OS X, and Linux aslong as they have, or there exists a version of:• Libsndfile for I/O;• FFTW3 library for CPU implementation;• CUDA Framework;• OpenCL driver.The program can be adapted by removing functionalities pro-vided by any subset of the previous requirements by remov-ing the components that make use of that prerequisite. Thesource code is available from the author athttp://www.lim.dico.unimi.it/CGPUconv.4.1 Performance ComparisonsPerformances of these algorithms depends on the size ofinput. Therefore, to characterize the “trade-off”, we testedthem with different input sizes. To make a reliable com-parison we choose to use as input signals a logarithmic sinesweep and its TRM (time reversal mirror) so the output shouldbe the δ function (Dirac delta function) or, to be more pre-cise, the limited bandwidth approximation of the sinc (sinuscardinalis) function.δ(x) =⎧⎪⎪⎨⎪⎪⎩+∞, x = 00, x  0(7)∫ ∞−∞δ(x)dx = 1 (8)sinc(x) =sin(x)x(9)We then compute the time spent on the convolution proce-dure, excluding the load procedure that reads from audio filesand the write to disk procedure for the results, which are col-lateral to our primary goal. A special case is represented bythe first execution for both the CUDA and OpenCL imple-mentation where for the former there exists some extra timedevoted to the load of the environment while for the latter,apart from the aforementioned setup, we have to take intoaccount the time that the driver allocate to compile kernelfunctions.Proceedings of the Acoustics 2012 Nantes Conference23-27 April 2012, Nantes, France2860Direct OLACPU - 9699CUDA - 6181OpenCL 7486 6699Table 1: Performance comparisons. Time in ms.The algorithms were executed on an OSX 10.6.8 equippedApple Macbook Pro 13.3” (MacBookPro5,5), Intel Core 2Duo processor @2,53 GHz, 8 GB Ram, NVIDIA GeForce9400GM VRAM 256 MB shared memory. OpenCL driversare provided by the operating system (1.1 compatible), andthe CUDA framework is version 4.0.All the audio files are high quality PCM uncompressedfiles and have a sample rate of 96 kHz and a quantizationword of 24 bit. With this bit depth the theoretical dynamicrange is ∼ 144 dB.For each algorithm we measured the difference computedbetween the signal under test and the reference (coming fromthe Matlab implementation) with a phase inversion. So thedifference on a sample by sample basis gives us a new signalthat can be used as a degree of similarity between the twooriginal signals. For each and every proposed approach thissignal is below -122 dB FS (dB on the full scale) meaningthere is no practical difference, and the result is in the orderof magnitude of the noise floor.Coming to the execution time of the algorithms we pro-pose a summary of the results presented in Figures 4, 5.Results are depicted as a function of the number of inputsamples, averaged over 100 runs.We also present in Table 1 results for a “real-case sce-nario”. We have a violin sound that is three minutes long anda reverberant impulse response of 1 s (sample rate 96 kHz)• Input: 17703123 samples (∼3’10”)• Kernel: 96000 (∼1”)Please note that “-” occurs when there is not enough freevideo RAM to handle the data. The idea here is to have a sys-tem that can run on most home computer so the relatively oldand low powerful graphic card is a good example of what canbe achieved with standard equipment. There are differencebetween implementations and this can be explained by thedifferent way of encoding real and complex numbers. Alsonote that there does not exist a concept of “paging” for videoRAM so if a structure is too big to fit in memory there is noautomatic way to handle the situation.5 Summary and Discussion of the re-sultsIn this paper we presented a number of prototypes thatare suitable for spatialization of sounds exploiting the poten-tialities of GPUs. Some issues are still present but we wantto point out that the basic concepts here expressed are validand mark a profitable direction.Performance results suggest that for a number of real caseapplications there are benefits that can be at least of 1/3 ofthe execution time (compared to the reference CPU imple-mentation) and can be further improved with other GPU-specific, but not hardware specific, optimizations. Benefitsare increasingly evident as the size of the filter kernel growsand this is particularly useful for convolution with long rever-berant impulse responses (e.g. BRIRs) that can be employedin order to render real environments.AcknowledgmentsThe author gratefully wish to acknowledge their stimu-lating discussions with G. Haus director of the Laboratoryof Music Informatics (LIM) of the Universita` degli Studi diMilano and the many researchers and graduate students ofLIM. The author also would like to thank Lorenzo Picinaliand Brian F.G. Katz for extensive conversations concerningthese researches.This work has been partially funded by the Enhanced Mu-sic Interactive Platform for Internet User (EMIPIU) project.References[1] E. Castro Lopo. Libsndfile [computer software]. Re-trieved December, 28:2005, 2005.[2] E. Y. Choueiri. Optimal crosstalk cancellation for binau-ral audio with two loudspeakers. 2011.[3] J. Cooley and J. Tukey. An algorithm for the machinecalculation of complex Fourier series. Math. Comput,19(90):297–301, 1965.[4] B. Cowan and B. Kapralos. Spatial sound for videogames and virtual environments utilizing real-time GPU-based convolution. In Proceedings of the ACM Future-Play 2008 International Conference on the Future ofGame Design and Technology, pages 166–172, Toronto,Ontario, Canada, November 3-5 2008.[5] B. Cowan and B. Kapralos. Real-time GPU-based con-volution: a follow-up. In Proceedings of the ACM Fu-turePlay @ GDC Canada 2009 International Confer-ence on the Future of Game Design and Technology,pages 25–26, Vancouver, British Columbia, Canada,May 12–13 2009.[6] M. Frigo and S. Johnson. The design and implemen-tation of FFTW3. Proc. IEEE (Special Issue on Pro-gram Generation, Optimization, and Platform Adapta-tion), 93:216–231, 2005.[7] D. A. Mauro. Effetti della distanza nella spazializzazionee localizzazione binaurale. B.A. Thesis, Universita` degliStudi di Milano, July 2006.[8] A. Torger. BruteFIR - an open-source general-purpose audio convolver.http://www.ludd.luth.se/torger/brutefir.html.Proceedings of the Acoustics 2012 Nantes Conference 23-27 April 2012, Nantes, France2861 0 5000 10000 15000 20000 2500027 28 29 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224Execution time (ms)N. of samplesComparation of execution time for Direct ModeCPUCUDAOpenCLFigure 4: Execution time for Direct mode depending on input size. 0 1000 2000 3000 4000 5000 6000 7000 8000 900027 28 29 210 211 212 213 214 215 216 217 218 219 220 221 222 223Execution time (ms)N. of samplesComparation of execution time for Overlap-addCPUCUDAOpenCLFigure 5: Execution time for Overlap-add depending on input size.Proceedings of the Acoustics 2012 Nantes Conference23-27 April 2012, Nantes, France2862

Audio convolution by the mean of GPU: CUDA and OpenCL implementations

https://mds.marshall.edu/cgi/viewcontent.cgi?article=1021&context=wdcs_faculty

Audio convolution by the mean of GPU: CUDA and OpenCL implementations

Abstract

Similar works

Full text

Available Versions

Marshall University