diff --git a/.gitmodules b/.gitmodules index f604849..72ddb71 100644 --- a/.gitmodules +++ b/.gitmodules @@ -4,3 +4,6 @@ [submodule "lib/MIPP"] path = lib/MIPP url = https://github.com/hayguen/MIPP.git +[submodule "lib/gaborator"] + path = lib/gaborator + url = https://git.gammaspectra.live/S.O.N.G/TheGaborator.git diff --git a/lib/gaborator b/lib/gaborator new file mode 160000 index 0000000..80ef6d4 --- /dev/null +++ b/lib/gaborator @@ -0,0 +1 @@ +Subproject commit 80ef6d4c8703d165660f7ef8d321fce9b526adfc diff --git a/lib/gaborator/CHANGES b/lib/gaborator/CHANGES deleted file mode 100644 index 3294dc5..0000000 --- a/lib/gaborator/CHANGES +++ /dev/null @@ -1,99 +0,0 @@ - -1.7 - -Miscellaneous bug fixes. - -Support lower numbers of bands per octave, down to 4. - -Further improve the performance of analyzing short signal blocks. - -The "Frequency-Domain Filtering" and "Streaming" examples now use -a white noise and impulse signal, respectively. - -1.6 - -Add "API Introduction" documentation section that was missing -from version 1.5, causing broken links. - -Improve analysis and resynthesis performance when using PFFFT or vDSP -by automatically enabling the use of real rather than complex FFTs -where applicable. - -1.5 - -Add navigation links to the HTML documentation. - -Add a code example demonstrating synthesis of musical notes. - -Add a function process() for iterating over coefficients sets with -greater flexibility than apply(). Also add a function fill() for -algorithmically creating new coefficients. - -Make the C++ declarations in the API reference documents more closely -resemble actual C++ code. - -Add a method gaborator::analyzer::band_ref() returning the band number -corresponding to the reference frequency. - -1.4 - -Support building the library as C++17, while retaining compatibility -with C++11. - -Further improve the performance of analyzing short signal blocks, and -of signal blocks not aligned to large powers of two. - -Add a code example mesasuring the resynthesis signal-to-noise -ratio (SNR). - -1.3 - -Eliminate some compiler warnings. - -Declare gaborator::analyzer::band_ff() const, making the code match -the documentation. - -Fix incorrect return type of gaborator::analyzer::band_ff() in the -documentation. - -Improve performance of analyzing short signal blocks. - -Remove special-case optimization of analyzing signal slices of all -zeros, as it caused incorrect results. - -Support up to 384 bands per octave. - -1.2 - -Add overview documentation. - -Add real-time FAQ. - -Actually include version.h in the release. - -Fix off-by-one error in defintion of analyzer constructor ff_min -argument. - -Fix incorrect return value of band_ff() for DC band. - -Add streaming example code. - -Add analyzer::analysis_support() and analyzer::synthesis_support(). - -Document analyzer::band_ff(). - -Improve signal to noise ratio at low numbers of bands per octave. - -Note the need for -mfpu=neon on ARM in render.html. - -1.1 - -Added CHANGES file. - -Added reference documentation. - -New include file gaborator/version.h. - -1.0 - -Initial release diff --git a/lib/gaborator/LICENSE b/lib/gaborator/LICENSE deleted file mode 100644 index abcfc9f..0000000 --- a/lib/gaborator/LICENSE +++ /dev/null @@ -1,11 +0,0 @@ - -The Gaborator library is Copyright (C) 1992-2019 Andreas Gustafsson. - -License to distribute and modify the code is hereby granted under the -terms of the GNU Affero General Public License, version 3 (henceforth, -the AGPLv3), but not under other versions of the AGPL. See the file -doc/agpl-3.0.txt for the full text of the AGPLv3. - -If the terms of the AGPLv3 are not acceptable to you, commercial -licensing under different terms is possible. Please contact -info@gaborator.com for more information. diff --git a/lib/gaborator/README b/lib/gaborator/README deleted file mode 100644 index 582593b..0000000 --- a/lib/gaborator/README +++ /dev/null @@ -1 +0,0 @@ -See doc/index.html for HTML documentation. diff --git a/lib/gaborator/doc/agpl-3.0.txt b/lib/gaborator/doc/agpl-3.0.txt deleted file mode 100644 index be3f7b2..0000000 --- a/lib/gaborator/doc/agpl-3.0.txt +++ /dev/null @@ -1,661 +0,0 @@ - GNU AFFERO GENERAL PUBLIC LICENSE - Version 3, 19 November 2007 - - Copyright (C) 2007 Free Software Foundation, Inc. - Everyone is permitted to copy and distribute verbatim copies - of this license document, but changing it is not allowed. - - Preamble - - The GNU Affero General Public License is a free, copyleft license for -software and other kinds of works, specifically designed to ensure -cooperation with the community in the case of network server software. - - The licenses for most software and other practical works are designed -to take away your freedom to share and change the works. By contrast, -our General Public Licenses are intended to guarantee your freedom to -share and change all versions of a program--to make sure it remains free -software for all its users. - - When we speak of free software, we are referring to freedom, not -price. Our General Public Licenses are designed to make sure that you -have the freedom to distribute copies of free software (and charge for -them if you wish), that you receive source code or can get it if you -want it, that you can change the software or use pieces of it in new -free programs, and that you know you can do these things. - - Developers that use our General Public Licenses protect your rights -with two steps: (1) assert copyright on the software, and (2) offer -you this License which gives you legal permission to copy, distribute -and/or modify the software. - - A secondary benefit of defending all users' freedom is that -improvements made in alternate versions of the program, if they -receive widespread use, become available for other developers to -incorporate. Many developers of free software are heartened and -encouraged by the resulting cooperation. However, in the case of -software used on network servers, this result may fail to come about. -The GNU General Public License permits making a modified version and -letting the public access it on a server without ever releasing its -source code to the public. - - The GNU Affero General Public License is designed specifically to -ensure that, in such cases, the modified source code becomes available -to the community. It requires the operator of a network server to -provide the source code of the modified version running there to the -users of that server. Therefore, public use of a modified version, on -a publicly accessible server, gives the public access to the source -code of the modified version. - - An older license, called the Affero General Public License and -published by Affero, was designed to accomplish similar goals. This is -a different license, not a version of the Affero GPL, but Affero has -released a new version of the Affero GPL which permits relicensing under -this license. - - The precise terms and conditions for copying, distribution and -modification follow. - - TERMS AND CONDITIONS - - 0. Definitions. - - "This License" refers to version 3 of the GNU Affero General Public License. - - "Copyright" also means copyright-like laws that apply to other kinds of -works, such as semiconductor masks. - - "The Program" refers to any copyrightable work licensed under this -License. Each licensee is addressed as "you". "Licensees" and -"recipients" may be individuals or organizations. - - To "modify" a work means to copy from or adapt all or part of the work -in a fashion requiring copyright permission, other than the making of an -exact copy. The resulting work is called a "modified version" of the -earlier work or a work "based on" the earlier work. - - A "covered work" means either the unmodified Program or a work based -on the Program. - - To "propagate" a work means to do anything with it that, without -permission, would make you directly or secondarily liable for -infringement under applicable copyright law, except executing it on a -computer or modifying a private copy. Propagation includes copying, -distribution (with or without modification), making available to the -public, and in some countries other activities as well. - - To "convey" a work means any kind of propagation that enables other -parties to make or receive copies. Mere interaction with a user through -a computer network, with no transfer of a copy, is not conveying. - - An interactive user interface displays "Appropriate Legal Notices" -to the extent that it includes a convenient and prominently visible -feature that (1) displays an appropriate copyright notice, and (2) -tells the user that there is no warranty for the work (except to the -extent that warranties are provided), that licensees may convey the -work under this License, and how to view a copy of this License. If -the interface presents a list of user commands or options, such as a -menu, a prominent item in the list meets this criterion. - - 1. Source Code. - - The "source code" for a work means the preferred form of the work -for making modifications to it. "Object code" means any non-source -form of a work. - - A "Standard Interface" means an interface that either is an official -standard defined by a recognized standards body, or, in the case of -interfaces specified for a particular programming language, one that -is widely used among developers working in that language. - - The "System Libraries" of an executable work include anything, other -than the work as a whole, that (a) is included in the normal form of -packaging a Major Component, but which is not part of that Major -Component, and (b) serves only to enable use of the work with that -Major Component, or to implement a Standard Interface for which an -implementation is available to the public in source code form. A -"Major Component", in this context, means a major essential component -(kernel, window system, and so on) of the specific operating system -(if any) on which the executable work runs, or a compiler used to -produce the work, or an object code interpreter used to run it. - - The "Corresponding Source" for a work in object code form means all -the source code needed to generate, install, and (for an executable -work) run the object code and to modify the work, including scripts to -control those activities. However, it does not include the work's -System Libraries, or general-purpose tools or generally available free -programs which are used unmodified in performing those activities but -which are not part of the work. For example, Corresponding Source -includes interface definition files associated with source files for -the work, and the source code for shared libraries and dynamically -linked subprograms that the work is specifically designed to require, -such as by intimate data communication or control flow between those -subprograms and other parts of the work. - - The Corresponding Source need not include anything that users -can regenerate automatically from other parts of the Corresponding -Source. - - The Corresponding Source for a work in source code form is that -same work. - - 2. Basic Permissions. - - All rights granted under this License are granted for the term of -copyright on the Program, and are irrevocable provided the stated -conditions are met. This License explicitly affirms your unlimited -permission to run the unmodified Program. The output from running a -covered work is covered by this License only if the output, given its -content, constitutes a covered work. This License acknowledges your -rights of fair use or other equivalent, as provided by copyright law. - - You may make, run and propagate covered works that you do not -convey, without conditions so long as your license otherwise remains -in force. You may convey covered works to others for the sole purpose -of having them make modifications exclusively for you, or provide you -with facilities for running those works, provided that you comply with -the terms of this License in conveying all material for which you do -not control copyright. Those thus making or running the covered works -for you must do so exclusively on your behalf, under your direction -and control, on terms that prohibit them from making any copies of -your copyrighted material outside their relationship with you. - - Conveying under any other circumstances is permitted solely under -the conditions stated below. Sublicensing is not allowed; section 10 -makes it unnecessary. - - 3. Protecting Users' Legal Rights From Anti-Circumvention Law. - - No covered work shall be deemed part of an effective technological -measure under any applicable law fulfilling obligations under article -11 of the WIPO copyright treaty adopted on 20 December 1996, or -similar laws prohibiting or restricting circumvention of such -measures. - - When you convey a covered work, you waive any legal power to forbid -circumvention of technological measures to the extent such circumvention -is effected by exercising rights under this License with respect to -the covered work, and you disclaim any intention to limit operation or -modification of the work as a means of enforcing, against the work's -users, your or third parties' legal rights to forbid circumvention of -technological measures. - - 4. Conveying Verbatim Copies. - - You may convey verbatim copies of the Program's source code as you -receive it, in any medium, provided that you conspicuously and -appropriately publish on each copy an appropriate copyright notice; -keep intact all notices stating that this License and any -non-permissive terms added in accord with section 7 apply to the code; -keep intact all notices of the absence of any warranty; and give all -recipients a copy of this License along with the Program. - - You may charge any price or no price for each copy that you convey, -and you may offer support or warranty protection for a fee. - - 5. Conveying Modified Source Versions. - - You may convey a work based on the Program, or the modifications to -produce it from the Program, in the form of source code under the -terms of section 4, provided that you also meet all of these conditions: - - a) The work must carry prominent notices stating that you modified - it, and giving a relevant date. - - b) The work must carry prominent notices stating that it is - released under this License and any conditions added under section - 7. This requirement modifies the requirement in section 4 to - "keep intact all notices". - - c) You must license the entire work, as a whole, under this - License to anyone who comes into possession of a copy. This - License will therefore apply, along with any applicable section 7 - additional terms, to the whole of the work, and all its parts, - regardless of how they are packaged. This License gives no - permission to license the work in any other way, but it does not - invalidate such permission if you have separately received it. - - d) If the work has interactive user interfaces, each must display - Appropriate Legal Notices; however, if the Program has interactive - interfaces that do not display Appropriate Legal Notices, your - work need not make them do so. - - A compilation of a covered work with other separate and independent -works, which are not by their nature extensions of the covered work, -and which are not combined with it such as to form a larger program, -in or on a volume of a storage or distribution medium, is called an -"aggregate" if the compilation and its resulting copyright are not -used to limit the access or legal rights of the compilation's users -beyond what the individual works permit. Inclusion of a covered work -in an aggregate does not cause this License to apply to the other -parts of the aggregate. - - 6. Conveying Non-Source Forms. - - You may convey a covered work in object code form under the terms -of sections 4 and 5, provided that you also convey the -machine-readable Corresponding Source under the terms of this License, -in one of these ways: - - a) Convey the object code in, or embodied in, a physical product - (including a physical distribution medium), accompanied by the - Corresponding Source fixed on a durable physical medium - customarily used for software interchange. - - b) Convey the object code in, or embodied in, a physical product - (including a physical distribution medium), accompanied by a - written offer, valid for at least three years and valid for as - long as you offer spare parts or customer support for that product - model, to give anyone who possesses the object code either (1) a - copy of the Corresponding Source for all the software in the - product that is covered by this License, on a durable physical - medium customarily used for software interchange, for a price no - more than your reasonable cost of physically performing this - conveying of source, or (2) access to copy the - Corresponding Source from a network server at no charge. - - c) Convey individual copies of the object code with a copy of the - written offer to provide the Corresponding Source. This - alternative is allowed only occasionally and noncommercially, and - only if you received the object code with such an offer, in accord - with subsection 6b. - - d) Convey the object code by offering access from a designated - place (gratis or for a charge), and offer equivalent access to the - Corresponding Source in the same way through the same place at no - further charge. You need not require recipients to copy the - Corresponding Source along with the object code. If the place to - copy the object code is a network server, the Corresponding Source - may be on a different server (operated by you or a third party) - that supports equivalent copying facilities, provided you maintain - clear directions next to the object code saying where to find the - Corresponding Source. Regardless of what server hosts the - Corresponding Source, you remain obligated to ensure that it is - available for as long as needed to satisfy these requirements. - - e) Convey the object code using peer-to-peer transmission, provided - you inform other peers where the object code and Corresponding - Source of the work are being offered to the general public at no - charge under subsection 6d. - - A separable portion of the object code, whose source code is excluded -from the Corresponding Source as a System Library, need not be -included in conveying the object code work. - - A "User Product" is either (1) a "consumer product", which means any -tangible personal property which is normally used for personal, family, -or household purposes, or (2) anything designed or sold for incorporation -into a dwelling. In determining whether a product is a consumer product, -doubtful cases shall be resolved in favor of coverage. For a particular -product received by a particular user, "normally used" refers to a -typical or common use of that class of product, regardless of the status -of the particular user or of the way in which the particular user -actually uses, or expects or is expected to use, the product. A product -is a consumer product regardless of whether the product has substantial -commercial, industrial or non-consumer uses, unless such uses represent -the only significant mode of use of the product. - - "Installation Information" for a User Product means any methods, -procedures, authorization keys, or other information required to install -and execute modified versions of a covered work in that User Product from -a modified version of its Corresponding Source. The information must -suffice to ensure that the continued functioning of the modified object -code is in no case prevented or interfered with solely because -modification has been made. - - If you convey an object code work under this section in, or with, or -specifically for use in, a User Product, and the conveying occurs as -part of a transaction in which the right of possession and use of the -User Product is transferred to the recipient in perpetuity or for a -fixed term (regardless of how the transaction is characterized), the -Corresponding Source conveyed under this section must be accompanied -by the Installation Information. But this requirement does not apply -if neither you nor any third party retains the ability to install -modified object code on the User Product (for example, the work has -been installed in ROM). - - The requirement to provide Installation Information does not include a -requirement to continue to provide support service, warranty, or updates -for a work that has been modified or installed by the recipient, or for -the User Product in which it has been modified or installed. Access to a -network may be denied when the modification itself materially and -adversely affects the operation of the network or violates the rules and -protocols for communication across the network. - - Corresponding Source conveyed, and Installation Information provided, -in accord with this section must be in a format that is publicly -documented (and with an implementation available to the public in -source code form), and must require no special password or key for -unpacking, reading or copying. - - 7. Additional Terms. - - "Additional permissions" are terms that supplement the terms of this -License by making exceptions from one or more of its conditions. -Additional permissions that are applicable to the entire Program shall -be treated as though they were included in this License, to the extent -that they are valid under applicable law. If additional permissions -apply only to part of the Program, that part may be used separately -under those permissions, but the entire Program remains governed by -this License without regard to the additional permissions. - - When you convey a copy of a covered work, you may at your option -remove any additional permissions from that copy, or from any part of -it. (Additional permissions may be written to require their own -removal in certain cases when you modify the work.) You may place -additional permissions on material, added by you to a covered work, -for which you have or can give appropriate copyright permission. - - Notwithstanding any other provision of this License, for material you -add to a covered work, you may (if authorized by the copyright holders of -that material) supplement the terms of this License with terms: - - a) Disclaiming warranty or limiting liability differently from the - terms of sections 15 and 16 of this License; or - - b) Requiring preservation of specified reasonable legal notices or - author attributions in that material or in the Appropriate Legal - Notices displayed by works containing it; or - - c) Prohibiting misrepresentation of the origin of that material, or - requiring that modified versions of such material be marked in - reasonable ways as different from the original version; or - - d) Limiting the use for publicity purposes of names of licensors or - authors of the material; or - - e) Declining to grant rights under trademark law for use of some - trade names, trademarks, or service marks; or - - f) Requiring indemnification of licensors and authors of that - material by anyone who conveys the material (or modified versions of - it) with contractual assumptions of liability to the recipient, for - any liability that these contractual assumptions directly impose on - those licensors and authors. - - All other non-permissive additional terms are considered "further -restrictions" within the meaning of section 10. If the Program as you -received it, or any part of it, contains a notice stating that it is -governed by this License along with a term that is a further -restriction, you may remove that term. If a license document contains -a further restriction but permits relicensing or conveying under this -License, you may add to a covered work material governed by the terms -of that license document, provided that the further restriction does -not survive such relicensing or conveying. - - If you add terms to a covered work in accord with this section, you -must place, in the relevant source files, a statement of the -additional terms that apply to those files, or a notice indicating -where to find the applicable terms. - - Additional terms, permissive or non-permissive, may be stated in the -form of a separately written license, or stated as exceptions; -the above requirements apply either way. - - 8. Termination. - - You may not propagate or modify a covered work except as expressly -provided under this License. Any attempt otherwise to propagate or -modify it is void, and will automatically terminate your rights under -this License (including any patent licenses granted under the third -paragraph of section 11). - - However, if you cease all violation of this License, then your -license from a particular copyright holder is reinstated (a) -provisionally, unless and until the copyright holder explicitly and -finally terminates your license, and (b) permanently, if the copyright -holder fails to notify you of the violation by some reasonable means -prior to 60 days after the cessation. - - Moreover, your license from a particular copyright holder is -reinstated permanently if the copyright holder notifies you of the -violation by some reasonable means, this is the first time you have -received notice of violation of this License (for any work) from that -copyright holder, and you cure the violation prior to 30 days after -your receipt of the notice. - - Termination of your rights under this section does not terminate the -licenses of parties who have received copies or rights from you under -this License. If your rights have been terminated and not permanently -reinstated, you do not qualify to receive new licenses for the same -material under section 10. - - 9. Acceptance Not Required for Having Copies. - - You are not required to accept this License in order to receive or -run a copy of the Program. Ancillary propagation of a covered work -occurring solely as a consequence of using peer-to-peer transmission -to receive a copy likewise does not require acceptance. However, -nothing other than this License grants you permission to propagate or -modify any covered work. These actions infringe copyright if you do -not accept this License. Therefore, by modifying or propagating a -covered work, you indicate your acceptance of this License to do so. - - 10. Automatic Licensing of Downstream Recipients. - - Each time you convey a covered work, the recipient automatically -receives a license from the original licensors, to run, modify and -propagate that work, subject to this License. You are not responsible -for enforcing compliance by third parties with this License. - - An "entity transaction" is a transaction transferring control of an -organization, or substantially all assets of one, or subdividing an -organization, or merging organizations. If propagation of a covered -work results from an entity transaction, each party to that -transaction who receives a copy of the work also receives whatever -licenses to the work the party's predecessor in interest had or could -give under the previous paragraph, plus a right to possession of the -Corresponding Source of the work from the predecessor in interest, if -the predecessor has it or can get it with reasonable efforts. - - You may not impose any further restrictions on the exercise of the -rights granted or affirmed under this License. For example, you may -not impose a license fee, royalty, or other charge for exercise of -rights granted under this License, and you may not initiate litigation -(including a cross-claim or counterclaim in a lawsuit) alleging that -any patent claim is infringed by making, using, selling, offering for -sale, or importing the Program or any portion of it. - - 11. Patents. - - A "contributor" is a copyright holder who authorizes use under this -License of the Program or a work on which the Program is based. The -work thus licensed is called the contributor's "contributor version". - - A contributor's "essential patent claims" are all patent claims -owned or controlled by the contributor, whether already acquired or -hereafter acquired, that would be infringed by some manner, permitted -by this License, of making, using, or selling its contributor version, -but do not include claims that would be infringed only as a -consequence of further modification of the contributor version. For -purposes of this definition, "control" includes the right to grant -patent sublicenses in a manner consistent with the requirements of -this License. - - Each contributor grants you a non-exclusive, worldwide, royalty-free -patent license under the contributor's essential patent claims, to -make, use, sell, offer for sale, import and otherwise run, modify and -propagate the contents of its contributor version. - - In the following three paragraphs, a "patent license" is any express -agreement or commitment, however denominated, not to enforce a patent -(such as an express permission to practice a patent or covenant not to -sue for patent infringement). To "grant" such a patent license to a -party means to make such an agreement or commitment not to enforce a -patent against the party. - - If you convey a covered work, knowingly relying on a patent license, -and the Corresponding Source of the work is not available for anyone -to copy, free of charge and under the terms of this License, through a -publicly available network server or other readily accessible means, -then you must either (1) cause the Corresponding Source to be so -available, or (2) arrange to deprive yourself of the benefit of the -patent license for this particular work, or (3) arrange, in a manner -consistent with the requirements of this License, to extend the patent -license to downstream recipients. "Knowingly relying" means you have -actual knowledge that, but for the patent license, your conveying the -covered work in a country, or your recipient's use of the covered work -in a country, would infringe one or more identifiable patents in that -country that you have reason to believe are valid. - - If, pursuant to or in connection with a single transaction or -arrangement, you convey, or propagate by procuring conveyance of, a -covered work, and grant a patent license to some of the parties -receiving the covered work authorizing them to use, propagate, modify -or convey a specific copy of the covered work, then the patent license -you grant is automatically extended to all recipients of the covered -work and works based on it. - - A patent license is "discriminatory" if it does not include within -the scope of its coverage, prohibits the exercise of, or is -conditioned on the non-exercise of one or more of the rights that are -specifically granted under this License. You may not convey a covered -work if you are a party to an arrangement with a third party that is -in the business of distributing software, under which you make payment -to the third party based on the extent of your activity of conveying -the work, and under which the third party grants, to any of the -parties who would receive the covered work from you, a discriminatory -patent license (a) in connection with copies of the covered work -conveyed by you (or copies made from those copies), or (b) primarily -for and in connection with specific products or compilations that -contain the covered work, unless you entered into that arrangement, -or that patent license was granted, prior to 28 March 2007. - - Nothing in this License shall be construed as excluding or limiting -any implied license or other defenses to infringement that may -otherwise be available to you under applicable patent law. - - 12. No Surrender of Others' Freedom. - - If conditions are imposed on you (whether by court order, agreement or -otherwise) that contradict the conditions of this License, they do not -excuse you from the conditions of this License. If you cannot convey a -covered work so as to satisfy simultaneously your obligations under this -License and any other pertinent obligations, then as a consequence you may -not convey it at all. For example, if you agree to terms that obligate you -to collect a royalty for further conveying from those to whom you convey -the Program, the only way you could satisfy both those terms and this -License would be to refrain entirely from conveying the Program. - - 13. Remote Network Interaction; Use with the GNU General Public License. - - Notwithstanding any other provision of this License, if you modify the -Program, your modified version must prominently offer all users -interacting with it remotely through a computer network (if your version -supports such interaction) an opportunity to receive the Corresponding -Source of your version by providing access to the Corresponding Source -from a network server at no charge, through some standard or customary -means of facilitating copying of software. This Corresponding Source -shall include the Corresponding Source for any work covered by version 3 -of the GNU General Public License that is incorporated pursuant to the -following paragraph. - - Notwithstanding any other provision of this License, you have -permission to link or combine any covered work with a work licensed -under version 3 of the GNU General Public License into a single -combined work, and to convey the resulting work. The terms of this -License will continue to apply to the part which is the covered work, -but the work with which it is combined will remain governed by version -3 of the GNU General Public License. - - 14. Revised Versions of this License. - - The Free Software Foundation may publish revised and/or new versions of -the GNU Affero General Public License from time to time. Such new versions -will be similar in spirit to the present version, but may differ in detail to -address new problems or concerns. - - Each version is given a distinguishing version number. If the -Program specifies that a certain numbered version of the GNU Affero General -Public License "or any later version" applies to it, you have the -option of following the terms and conditions either of that numbered -version or of any later version published by the Free Software -Foundation. If the Program does not specify a version number of the -GNU Affero General Public License, you may choose any version ever published -by the Free Software Foundation. - - If the Program specifies that a proxy can decide which future -versions of the GNU Affero General Public License can be used, that proxy's -public statement of acceptance of a version permanently authorizes you -to choose that version for the Program. - - Later license versions may give you additional or different -permissions. However, no additional obligations are imposed on any -author or copyright holder as a result of your choosing to follow a -later version. - - 15. Disclaimer of Warranty. - - THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY -APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT -HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY -OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, -THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM -IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF -ALL NECESSARY SERVICING, REPAIR OR CORRECTION. - - 16. Limitation of Liability. - - IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING -WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS -THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY -GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE -USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF -DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD -PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), -EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF -SUCH DAMAGES. - - 17. Interpretation of Sections 15 and 16. - - If the disclaimer of warranty and limitation of liability provided -above cannot be given local legal effect according to their terms, -reviewing courts shall apply local law that most closely approximates -an absolute waiver of all civil liability in connection with the -Program, unless a warranty or assumption of liability accompanies a -copy of the Program in return for a fee. - - END OF TERMS AND CONDITIONS - - How to Apply These Terms to Your New Programs - - If you develop a new program, and you want it to be of the greatest -possible use to the public, the best way to achieve this is to make it -free software which everyone can redistribute and change under these terms. - - To do so, attach the following notices to the program. It is safest -to attach them to the start of each source file to most effectively -state the exclusion of warranty; and each file should have at least -the "copyright" line and a pointer to where the full notice is found. - - - Copyright (C) - - This program is free software: you can redistribute it and/or modify - it under the terms of the GNU Affero General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Affero General Public License for more details. - - You should have received a copy of the GNU Affero General Public License - along with this program. If not, see . - -Also add information on how to contact you by electronic and paper mail. - - If your software can interact with users remotely through a computer -network, you should also make sure that it provides a way for users to -get its source. For example, if your program is a web application, its -interface could display a "Source" link that leads users to an archive -of the code. There are many ways you could offer source, and different -solutions will be better for different programs; see section 13 for the -specific requirements. - - You should also get your employer (if you work as a programmer) or school, -if any, to sign a "copyright disclaimer" for the program, if necessary. -For more information on this, and how to apply and follow the GNU AGPL, see -. diff --git a/lib/gaborator/doc/doc.css b/lib/gaborator/doc/doc.css deleted file mode 100644 index 86aab82..0000000 --- a/lib/gaborator/doc/doc.css +++ /dev/null @@ -1,53 +0,0 @@ -body { - background-color: #000; - font-family: sans-serif; - margin: 10%; - color: #eee; -} -a:link, a:visited, a:hover, a:active { - color: currentColor; -} -pre { - border: 1px solid #888; - margin: 20px; - margin-left: 0px; - padding: 10px; - background: #222; - /* To avoid text extending outside the border on narrow displays */ - white-space:pre-wrap; -} -img { - border: 1px solid #888; - margin: 20px; - margin-left: 0px; - padding: 10px; - background: #000; -} -h2 { - margin-top: 2em; -} -h3 { - margin-top: 1.5em; -} -/* http://code.stephenmorley.org/html-and-css/fixing-browsers-broken-monospace-font-handling/ */ -pre, code, kbd, samp, tt { - font-family:monospace,monospace; - font-size:1em; -} -pre.forward_decl { -/* Needed for syntax checking, but avoid clutter for human readers */ - display: none; -} -div.class_def { - margin-left: 2em; -} -div.nav { - margin-top: 30px; - font-style: oblique; -} -div.nav span.prev { - float: left; -} -div.nav span.next { - float: right; -} \ No newline at end of file diff --git a/lib/gaborator/doc/filter-response.png b/lib/gaborator/doc/filter-response.png deleted file mode 100644 index 7be505e..0000000 Binary files a/lib/gaborator/doc/filter-response.png and /dev/null differ diff --git a/lib/gaborator/doc/filter.html b/lib/gaborator/doc/filter.html deleted file mode 100644 index 25aa939..0000000 --- a/lib/gaborator/doc/filter.html +++ /dev/null @@ -1,266 +0,0 @@ - - - - - -Gaborator Example 2: Frequency-Domain Filtering - - -

Example 2: Frequency-Domain Filtering

- -

Introduction

- -

This example shows how to apply a filter to an audio file using -the Gaborator library, by turning the audio into spectrogram -coefficients, modifying the coefficients, and resynthesizing audio -from them.

- -

The specific filter implemented here is a 3 dB/octave lowpass -filter. This is sometimes called a pinking filter because it -can be used to produce pink noise from white noise. In practice, the -3 dB/octave slope is only applied above some minimum frequency, for -example 20 Hz, because otherwise the gain of the filter would approach -infinity as the frequency approaches 0, and the impulse response would -have to be infinitely wide. -

- -

Since the slope of this filter is not a multiple of 6 dB/octave, it -is difficult to implement as an analog filter, but by filtering -digitally in the frequency domain, arbitrary filter responses such as -this can easily be achieved. -

- -

Preamble

- -
-#include <memory.h>
-#include <iostream>
-#include <sndfile.h>
-#include <gaborator/gaborator.h>
-
-int main(int argc, char **argv) {
-    if (argc < 3) {
-        std::cerr << "usage: filter input.wav output.wav\n";
-        exit(1);
-    }
-
- -

Reading the Audio

- -

The code for reading the input audio file is identical to -that in Example 1:

- -
-    SF_INFO sfinfo;
-    memset(&sfinfo, 0, sizeof(sfinfo));
-    SNDFILE *sf_in = sf_open(argv[1], SFM_READ, &sfinfo);
-    if (! sf_in) {
-        std::cerr << "could not open input audio file: "
-            << sf_strerror(sf_in) << "\n";
-        exit(1);
-    }
-    double fs = sfinfo.samplerate;
-    sf_count_t n_frames = sfinfo.frames;
-    sf_count_t n_samples = sfinfo.frames * sfinfo.channels;
-    std::vector<float> audio(n_samples);
-    sf_count_t n_read = sf_readf_float(sf_in, audio.data(), n_frames);
-    if (n_read != n_frames) {
-        std::cerr << "read error\n";
-        exit(1);
-    }
-    sf_close(sf_in);
-
- -

Spectrum Analysis Parameters

- -

The spectrum analysis works much the same as in Example 1, -but uses slightly different parameters. -We use a larger number of frequency bands per octave (100) -to minimize ripple in the frequency response, and the -reference frequency argument is omitted as we don't care about the -exact alignment of the bands with respect to a musical scale.

-
-    gaborator::parameters params(100, 20.0 / fs);
-    gaborator::analyzer<float> analyzer(params);
-
- -

Precalculating Gains

- -

The filtering will be done by multiplying each spectrogram -coefficient with a frequency-dependent gain. To avoid having to -calculate the gain on the fly for each coefficient, which would -be slow, we will precalculate the gains into a vector band_gains -of one gain value per band, including one for the -special lowpass band that contains the frequencies from 0 to 20 Hz.

- -
-    std::vector<float> band_gains(analyzer.bands_end());
-
- -

First, we calculate the gains for the bandpass bands. -For a 3 dB/octave lowpass filter, the voltage gain needs to be -proportional to the square root of the inverse of the frequency. -To get the frequency of each band, we call the -analyzer method band_ff(), which -returns the center frequency of the band in units of the -sampling frequency. The gain is normalized to unity at 20 Hz. -

-
-    for (int band = analyzer.bandpass_bands_begin(); band < analyzer.bandpass_bands_end(); band++) {
-        double f_hz = analyzer.band_ff(band) * fs;
-        band_gains[band] = 1.0 / sqrt(f_hz / 20.0);
-    }
-
- -

The gain of the lowpass band is set to the the same value as the -lowest-frequency bandpass band, so that the overall filter gain -plateaus smoothly to a constant value below 20 Hz.

- -
-    band_gains[analyzer.band_lowpass()] = band_gains[analyzer.bandpass_bands_end() - 1];
-
- -

De-interleaving

- -

To handle stereo and other multi-channel audio files, -we will loop over the channels and filter each channel separately. -Since libsndfile produces interleaved samples, we first -de-interleave the current channel into a temporary vector called -channel:

-
-    for (sf_count_t ch = 0; ch < sfinfo.channels; ch++) {
-        std::vector<float> channel(n_frames);
-        for (sf_count_t i = 0; i < n_frames; i++)
-            channel[i] = audio[i * sfinfo.channels + ch];
-
-

Spectrum Analysis

-

Now we can spectrum analyze the current channel, producing -a set of coefficients:

-
-        gaborator::coefs<float> coefs(analyzer);
-        analyzer.analyze(channel.data(), 0, channel.size(), coefs);
-
- -

Filtering

-

-The filtering is done using the function -process(), which applies a user-defined function -to each spectrogram coefficient. Here, that user-defined function is a -lambda expression that multiplies the coefficient by the appropriate -precalculated frequency-dependent gain, modifying the coefficient in -place. The unused int64_t argument is the time in units -of samples; this could be use to implement a time-varying filter if -desired.

-

-The second and third argument to process() specify a -range of frequency bands to process; here we pass INT_MIN, -INT_MAX to process all of them. Similarly, the fourth and -fifth argument specify a time range to process, and we pass -INT64_MIN, INT64_MAX to process all the coefficients -in coefs regardless of time. -

-
-        process([&](int band, int64_t, std::complex<float> &coef) {
-                coef *= band_gains[band];
-            },
-            INT_MIN, INT_MAX,
-            INT64_MIN, INT64_MAX,
-            coefs);
-
- -

Resynthesis

-

We can now resynthesize audio from the filtered coefficients by -calling synthesize(). This is a mirror image of the call to -analyze(): now the coefficients are the input, and -the buffer of samples is the output. The channel -vector that originally contained the input samples for the channel -is now reused to hold the output samples.

-
-        analyzer.synthesize(coefs, 0, channel.size(), channel.data());
-
- -

Re-interleaving

-

The audio vector that contained the -original interleaved audio is reused for the interleaved -filtered audio. This concludes the loop over the channels. -

-
-        for (sf_count_t i = 0; i < n_frames; i++)
-            audio[i * sfinfo.channels + ch] = channel[i];
-    }
-
- -

Writing the Audio

-

The filtered audio is written using libsndfile, -using code that closely mirrors that for reading. -Note that we use SFC_SET_CLIPPING -to make sure that any samples too loud for the file format -will saturate; by default, libsndfile makes them -wrap around, which sounds really bad.

- -
-    SNDFILE *sf_out = sf_open(argv[2], SFM_WRITE, &sfinfo);
-    if (! sf_out) {
-        std::cerr << "could not open output audio file: "
-            << sf_strerror(sf_out) << "\n";
-        exit(1);
-    }
-    sf_command(sf_out, SFC_SET_CLIPPING, NULL, SF_TRUE);
-    sf_count_t n_written = sf_writef_float(sf_out, audio.data(), n_frames);
-    if (n_written != n_frames) {
-        std::cerr << "write error\n";
-        exit(1);
-    }
-    sf_close(sf_out);
-
-
- -

Postamble

-

-We need a couple more lines of boilerplate to make the example a -complete program: -

-
-    return 0;
-}
-
- -

Compiling

-

Like Example 1, this example -can be built using a one-line build command: -

-
-c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` filter.cc `pkg-config --libs sndfile` -o filter
-
-

Or using the vDSP FFT on macOS:

-
-c++ -std=c++11 -I.. -O3 -ffast-math -DGABORATOR_USE_VDSP `pkg-config --cflags sndfile` filter.cc `pkg-config --libs sndfile` -framework Accelerate -o filter
-
-

Or using PFFFT (see Example 1 for how to download and build PFFFT):

-
-c++ -std=c++11 -I.. -Ipffft -O3 -ffast-math -DGABORATOR_USE_PFFFT `pkg-config --cflags sndfile` filter.cc pffft/pffft.o pffft/fftpack.o `pkg-config --libs sndfile` -o filter
-
- -

Running

-

Running the following shell commands will download an example -audio file containing five seconds of white noise and filter it, -producing pink noise.

-
-wget http://download.gaborator.com/audio/white_noise.wav
-./filter white_noise.wav pink_noise.wav
-
- -

Frequency response

-

The following plot shows the actual measured frequency response of the -filter, with the expected 3 dB/octave slope above 20 Hz and minimal -ripple:

-Frequency response plot - - - - - diff --git a/lib/gaborator/doc/gen/allkernels_v1_bpo12_ffmin0.03_ffref0.5_anl_wob.png b/lib/gaborator/doc/gen/allkernels_v1_bpo12_ffmin0.03_ffref0.5_anl_wob.png deleted file mode 100644 index b75fa6c..0000000 Binary files a/lib/gaborator/doc/gen/allkernels_v1_bpo12_ffmin0.03_ffref0.5_anl_wob.png and /dev/null differ diff --git a/lib/gaborator/doc/gen/allkernels_v1_bpo12_ffmin0.03_ffref0.5_syn_wob.png b/lib/gaborator/doc/gen/allkernels_v1_bpo12_ffmin0.03_ffref0.5_syn_wob.png deleted file mode 100644 index 03ff1a1..0000000 Binary files a/lib/gaborator/doc/gen/allkernels_v1_bpo12_ffmin0.03_ffref0.5_syn_wob.png and /dev/null differ diff --git a/lib/gaborator/doc/gen/grid_v1_bpo12_ffmin0.03_ffref0.5_wob.png b/lib/gaborator/doc/gen/grid_v1_bpo12_ffmin0.03_ffref0.5_wob.png deleted file mode 100644 index 91fdf7b..0000000 Binary files a/lib/gaborator/doc/gen/grid_v1_bpo12_ffmin0.03_ffref0.5_wob.png and /dev/null differ diff --git a/lib/gaborator/doc/index.html b/lib/gaborator/doc/index.html deleted file mode 100644 index d7d3231..0000000 --- a/lib/gaborator/doc/index.html +++ /dev/null @@ -1,78 +0,0 @@ - - - - - - -The Gaborator - - -

The Gaborator

-

The Gaborator is a library that generates constant-Q spectrograms -for visualization and analysis of audio signals. It also supports a -fast and accurate inverse transformation of the spectrogram coefficients -back into audio for spectral effects and editing.

- -

The Gaborator implements the invertible constant-Q transform of -Velasco, Holighaus, Dörfler, and Grill, described in the papers - -Constructing an invertible constant-Q transform with nonstationary Gabor frames, 2011 -and -A Framework for invertible, real-time constant-Q transforms, 2012, -using Gaussian bandpass filters and an efficient multi-rate architecture. -

- -

The Gaborator is written in C++11 and compatible with C++14 and C++17. -It has been tested on macOS, Linux, NetBSD, FreeBSD, and iOS, on Intel -x86_64 and ARM processors.

- -

The Gaborator is open source under the GNU Affero General Public -License, version 3, and is also available for commercial licensing. -See the file LICENSE for details.

- -

Example Code

- -

The following examples demonstrate the use of the library in -various scenarios. They are presented in a "literate -programming" style, with the code embedded in the commentary -rather than the other way around. -Concatenating the code fragments in each example yields a complete C++ -program, which can also be found as a .cc file in -the examples/ directory.

- - -

API Reference

-

The following documents define the library API. -

- - -

How it Works

-

The following document outlines the operation of the library.

- - -

FAQ

- - -

Contact

-

Email questions and bug reports to the author at info@gaborator.com.

- - - diff --git a/lib/gaborator/doc/overview.html b/lib/gaborator/doc/overview.html deleted file mode 100644 index cfb3857..0000000 --- a/lib/gaborator/doc/overview.html +++ /dev/null @@ -1,135 +0,0 @@ - - - - - -Overview of Operation - - -

Overview of Operation

- -

The Gaborator performs three main functions:

-
    -
  • spectrum analysis, which turns a signal into a set -of spectrogram coefficients -
  • resynthesis (aka reconstruction), which turns a -set of coefficients back into a signal, and -
  • rendering, which -turns a set of coefficients into a rectangular array of -amplitude values that can be turned into pixels to display -a spectrogram. -
- -

The following sections give a high-level overview of each -of these functions.

- -

Analysis

- -

The first step of the analysis is to run the signal through -an analysis filter bank, to split it into a number of -overlapping frequency bands.

- -

The filter bank consists of a number of logarithmically spaced -Gaussian bandpass filters and a single lowpass filter. Each bandpass -filter has a bandwidth proportional to its center frequency, which -means they all have the same quality factor Q and form -a constant-Q filter bank. The highest-frequency bandpass -filter will have a center frequency close to half the sample rate; in -the graphs below, this is simple labeled 0.5 because all frequencies -in the Gaborator are in units of the sample rate. The -lowest-frequency bandpass filter should be centered at, or slightly -below, the lowest frequency of interest to the application at hand. -For example, when analyzing audio, this is often the lower limit of -human hearing; at a sample rate of 44100 Hz, this means 20 Hz / 44100 -Hz ≈ 0.00045. This lower frequency limit is referred to as -the minimum frequency or fmin. -

- -

Although frequencies below fmin are assumed to not be of -interest, they nonetheless need to be preserved to achieve perfect -reconstruction, and that is what the lowpass filter is for. Together, -the lowpass filter and the bandpass filters overlap to cover the full -frequency range from 0 to 0.5.

- -

The spacing of the bandpass filters is specified by the user as an -integer number of filters (or, equivalently, bands) per octave. For -example, when analyzing music, this is often 12 bands per octave (one -band per semitone in the equal-tempered scale), or if a finer -frequency resolution is needed, some multiple of 12.

- -

The following plot shows the frequency responses of the analysis -filters at 12 bands per octave and fmin = 0.03. A more -typical fmin for audio work would be 0.00045, but -that would make the plot hard to read because both the lowpass filter -and the lowest-frequency bandpass filters would be extremely narrow.

- -Analysis filters - -

The output of each bandpass filter is shifted down in frequency to -a complex quadrature baseband. The baseband signal is then resampled -at a reduced sample rate, lower than that of the orignal signal but -high enough that there is negligible aliasing given the bandwidth of -the filter in case. The Gaborator uses sample rates related to the -original signal sample rate by powers of two. This means some of -frequency bands are sampled a bit more often than strictly -necessary, but has the advantage that the sampling can be synchronized -to make the samples of many frequency bands coincide in time, which -can be convenient in later analysis or spectrogram rendering. The -complex samples resulting from this process are the spectrogram -coefficients.

- -

The center frequencies of the analysis filters and the points in -time at which they are sampled form a two-dimensional, -multi-resolution time-frequency grid, where high frequencies -are sampled sparsely in frequency but densely in time, and low -frequencies are sampled densely in frequency but sparsely in time.

- -

The following plot illustrates the time-frequency sampling grid -corresponding to the parameters used in the previous plot. Note that -frequency was the X axis in the previous plot, but is the Y axis -here. The plot covers a time range of 128 signal samples, but -conceptually, the grid extends arbitrarily far in time, in both the -positive and the negative direction.

- -Sampling grid - -

Resynthesis

- -

Resynthesizing a signal from the coefficients is more or less the -reverse of the analysis process. The coefficients are frequency -shifted from the complex baseband back to their original center -frequencies and run through a reconstruction filter bank -that is a dual of the analysis filter bank. The following -plot shows the frequency responses of the reconstruction filters -corresponding to the analysis filters shown earlier.

- -Reconstruction filters - -

Although the bandpass filters may look similar to the Gaussian -filters of the analysis filter bank, their shapes are actually subtly -different.

- -

Spectrogram Rendering

- -

Rendering a spectrogram image from the coefficients involves -taking the magnitude of each complex coefficient, and then -resampling the resulting multi-resolution grid of magnitudes -into an evenly spaced pixel grid.

- -

Because the coefficient sample rate varies by frequency band, the -resampling required in the horizontal (time) direction also varies. -Typically, the high-frequency bands of an audio spectrogram have more -than one coefficient per pixel and require downsampling (decimation), -some bands in the mid-range frequencies have a one-to-one relationship -between coefficients and pixels, and the low-frequency bands -have more than one pixel per coefficient and require upsampling -(interpolation).

- - - - - diff --git a/lib/gaborator/doc/realtime.html b/lib/gaborator/doc/realtime.html deleted file mode 100644 index 5ea5bec..0000000 --- a/lib/gaborator/doc/realtime.html +++ /dev/null @@ -1,126 +0,0 @@ - - - - - -Is it real-time? - - -

Is it real-time?

- -

Several people have asked whether the Gaborator is suitable for -real-time applications. There is no simple yes or no answer to -this question, because there are many different definitions of -"real-time", and the answer will depend the definition. -Below are some answers to the question "is it real-time?" -rephrased in terms of different definitions.

- -

Can it processes a recording in less time than its duration?

- -

Yes. For example, at 48 frequency bands per -octave, a single core of a 2.5 GHz Intel Core i5 CPU can analyze some -10 million samples per second, which is more than 200 times faster -than real time for a single channel of 44.1 kHz audio.

- -

Does it have bounded latency? -Can it start producing output before consuming the entire input? -Will it stream?

-

Yes. See the streaming example. - -

Does it have low latency?

- -

Probably not low enough for applications such as live musical -effects. The exact latency depends on factors such as the frequency -range and number of bands per octave, but tends to range between -"high" and "very high". For example, with the parameters used in the -online demo, 48 frequency bands per octave down to 20 Hz, the latency -of the analysis side alone is some 3.5 seconds, and if you do -analysis followed by resynthesis, the total latency will -be almost 13 seconds.

- -

This can be reduced by choosing the analysis parameters for low latency. -For example, if you decrease the number of frequency bands per octave to 12, -and increase the minimum frequency to 200 Hz, the latency -will be about 85 milliseconds for analysis only, and about -300 milliseconds for analysis + resynthesis, but this is -still too much for a live effect.

- -

Any constant-Q spectrum analysis involving low frequencies will -inherently have rather high latency (at least for musically useful -values of Q), because the lowest-frequency analysis filters will have -narrow bandwidths, which lead to long impulse responses. Furthermore, -the Gaborator uses symmetric Gaussian analysis filters that were -chosen for properties such as linear phase and accurate -reconstruction, not for low latency, so the latency will be higher -than what might be achievable with a constant-Q filter bank -specifically designed for low latency.

- -

The latency only affects causal applications, and -arises from the need to wait for the arrival of future input samples -needed to calculate the present output, and not from the time it takes -to perform the calculations. In a non-causal application, -such as applying an effect to a recording, the latency does not apply, -and performance is limited only by the speed of the calculations. -This can lead to the somewhat paradoxical situation that applying an -effect to a live stream causes a latency of several seconds, but -applying the same effect to an entire one-minute recording runs in a -fraction of a second.

- -

In analysis and visualization applications that don't need to -perform resynthesis, it is possible to partly hide the latency by -taking advantage of the fact that the coefficients for the higher -frequencies exhibit lower latency than those for low frequencies. -For example, a live spectrogram display could update the -high-frequency parts of the display before the corresponding -low-frequency parts. Alternatively, low-frequency parts of the -spectrogram may be drawn multiple times, effectively animating -the display of the low-frequency coefficients as they converge to -their final values. This approach can be seen in action in -the Spectrolite -iOS app.

- -

Does it support small blocks sizes?

- -

Yes, but there is a significant performance penalty. -The Gaborator works most efficiently when the signal is processed -in large blocks, preferably 217 samples or more, -corresponding to several seconds of signal at typical audio sample -rates.

- -

A real-time application aiming for low latency will want to -use smaller blocks, for examples 25 to 210 -samples, and processing these will be significantly slower. -For example, as of version 1.4, analyzing a signal in blocks of -210 samples takes roughly five times as much CPU as -analyzing it in blocks of 220 samples.

- -

For sufficiently small blocks, the processing time will exceed the -duration of the signal, at which point the system can no longer be -considered real-time. For example, analyzing a 48 kHz audio -stream on a 2.5 GHz Intel Core i5 CPU, this happens at block sizes -below about 24 = 16 samples.

- -

The resynthesis code is currently less optimized for small block -sizes than the analysis code, so the performance penalty for -resynthesizing small blocks is even greater than for analyzing small -blocks.

- -

Can it process a signal stream of any length?

- -

Not in practice — the length is limited by floating point -precision. At typical audio sample rates, roundoff errors start to -become significant after some hours.

- -

Does it avoid dynamic memory allocation in the audio processing path?

- -

Currently, no — it dynamically allocates both the coefficient data -structures and various temporary buffers.

- - - - - diff --git a/lib/gaborator/doc/ref/gaborator_h.html b/lib/gaborator/doc/ref/gaborator_h.html deleted file mode 100644 index 206901e..0000000 --- a/lib/gaborator/doc/ref/gaborator_h.html +++ /dev/null @@ -1,462 +0,0 @@ - - - - - -Gaborator reference: gaborator.h - - -

Gaborator reference: gaborator.h

- -

Spectrum Analysis Parameters

- -

A parameters object holds a set of parameters that -determine the frequency range and resolution of the spectrum -analysis.

- -
-class parameters {
-
- -
-

Constructor

-
-parameters(unsigned int bands_per_octave,
-           double ff_min,
-           double ff_ref = 1.0);
-
-
-
bands_per_octave
-
The number of frequency bands per octave. - Values from 4 to 384 (inclusive) are supported. -
-
ff_min
-
The lower limit of the analysis frequency range, in units of the - sample rate. The analysis filter bank will extend low enough in - frequency that ff_min falls between the two lowest - frequency bandpass filters. - Values from 0.001 to 0.13 are supported.
-
ff_ref
-
The reference frequency, in units of the sample rate. - This allows fine-tuning of the analysis and synthesis filter - banks such that the center frequency of one of the filters - is aligned with ff_ref. If ff_ref - falls outside the frequency range of the bandpass filter bank, this - works as if the range were extended to include - ff_ref. Must be positive. A typical value - when analyzing music is 440.0 / fs, where - fs is the sample rate in Hz. -
-
-

Comparison

-

-Comparison operators are provided for compatibility with -standard container classes. The ordering is arbitrary but consistent. -

-
-bool operator<(const parameters &rhs) const;
-bool operator==(const parameters &rhs) const;
-
- -
-
-};
-
- -

Spectrogram Coefficients

- -
-template<class T> class analyzer;
-
- -

-A coefs object stores a set of spectrogram coefficients. -It is a dynamic data structure and will be automatically grown to -accommodate new time ranges, for example as newly recorded audio is analyzed. -The template argument T -must match that of the analyzer (usually float). -The template argument C is the data type used to store each -coefficient value; there is usually no need to specify it explicitly as -it will default to std::complex<T>. -

- -
-template<class T, class C = std::complex<T>>
-class coefs {
-
-
-

Constructor

-
-coefs(analyzer<T> &a);
-
-

-Construct an empty set of coefficients for use with the spectrum -analyzer a. This represents a signal that is zero -at all points in time. -

- -
-
-};
-
- -

Spectrum Analyzer

- -

-The analyzer object performs spectrum analysis and/or resynthesis -according to the given parameters. The template argument T is -the floating-point type to use for the calculations. This is typically float; -alternatively, double can be used for increased accuracy at the -expense of speed and memory consumption.

- -
template<class T>
-class analyzer {
-
- -

Constructor

- -
-analyzer(const parameters &params);
-
-
-
params
-
The spectrum analysis parameters. -
- -

Analysis and synthesis

- -
-void
-analyze(const T *signal,
-        int64_t t0,
-        int64_t t1,
-        coefs<T> &coefs) const;
-
-

Spectrum analyze the samples at *signal and add the -resulting coefficients to coefs. -

-
signal
-
The signal samples to analyze, beginning with the sample from time t0 - and ending with the last sample before time t1, for a total of - t1 - t0 samples. -
t0
-
The point in time when the sample at signal[0] was taken, - in samples. For example, when analyzing an audio recording, this is typically - 0 for the first sample in the recording, but this reference point is arbitrary, - and negative times are valid. Accuracy begins to successively decrease - outside the range of about ±108 samples, so using - large time values should be avoided when they are not necessary because - of the length of the track. -
-
t1
-
The point in time of the sample one past the - end of the array of samples at signal, - in samples. -
-
coefs
The coefficient object that the results of the - spectrum analysis are added to. -
-

If the coefs object already contains some -coefficients, the new coefficients are summed to those already -present. Because the analysis is a linear operation, this allows a -signal to be analyzed in blocks, by making multiple calls -to analyze() with non-overlapping ranges that together -cover the entire signal. For efficiency, the blocks should -be large, as in -analyze(first_131072_samples, 0, 131072, coefs), -analyze(next_131072_samples, 131072, 262144, coefs), -etc. -

- -
-void
-synthesize(const coefs<T> &coefs,
-           uint64_t t0,
-           uint64_t t1,
-           T *signal) const;
-
-

Synthesize signal samples from the coefficients coef and store -them at *signal. -

-
-
coefs
The coefficients to synthesize the signal from.
-
t0
-
The point in time of the first sample to synthesize, - in samples, using the same time scale as in analyze().
-
t1
-
The point in time of the sample one past the last one to synthesize.
-
signal
-
The synthesized signal samples will be written here, - beginning with the sample from time t0 and - and ending with the last sample before time t1, - for a total of t1 - t0 samples.
-
-

The time range t0...t1 may extend outside -the range analyzed using analyze(), in which case the -signal is assumed to be zero in the un-analyzed range.

- -

A signal may be synthesized in blocks by making multiple calls to -analyze() with different sample ranges. For efficiency, -the blocks should be large, and each t0 should -be multiple of a large power of two.

- -

Frequency Band Numbering

- -

The frequency bands of the analysis filter bank are numbered by -nonnegative integers that increase towards lower (sic) frequencies. -There is a number of bandpass bands corresponding to the -logarithmically spaced bandpass analysis filters, from near 0.5 -(half the sample rate) to -near fmin, and a single lowpass band containing the -residual signal from frequencies below fmin. -The numbering can be examined using the following methods: -

- -
-int bandpass_bands_begin() const;
-
-

-Return the smallest valid bandpass band number, corresponding to the -highest-frequency bandpass filter.

-
-int bandpass_bands_end() const;
-
-

-Return the bandpass band number one past the highest valid bandpass -band number, corresponding to one past the lowest-frequency bandpass -filter. -

-
-int band_lowpass() const;
-
-

-Return the band number of the lowpass band. -

-
-int band_ref() const;
-
-

-Return the band number corresponding to the reference frequency -ff_ref. If ff_ref falls within -the frequency range of the bandpass filter bank, this will -be a valid bandpass band number, otherwise it will not. -

-
-double band_ff(int band) const;
-
-

-Return the center frequency of band number band, in units of the -sampling frequency. -

- -

Support

-
-double analysis_support() const;
-
-

Returns the one-sided worst-case time domain support of any of the -analysis filters. When calling analyze() with a sample at time t, -only spectrogram coefficients within the time range t ± support -will be significantly changed. Coefficients outside the range may change, -but the changes will sufficiently small that they may be ignored without -significantly reducing accuracy.

- -
-double synthesis_support() const;
-
-

Returns the one-sided worst-case time domain support of any of the -reconstruction filters. When calling synthesize() to -synthesize a sample at time t, the sample will only be -significantly affected by spectrogram coefficients in the time -range t ± support. Coefficients outside the range may -be used in the synthesis, but substituting zeroes for the actual -coefficient values will not significantly reduce accuracy.

- -
-
-};
-
- -

Functions

- -

Iterating Over Existing Coefficients

- -
-template <class T, class F, class C0, class... CI>
-void process(F f,
-             int b0,
-             int b1,
-             int64_t t0,
-             int64_t t1,
-             coefs<T, C0> &coefs0,
-             coefs<T, CI>&... coefsi);
-
- -

-Process one or more coefficient sets coefs0... by applying -the function f to each coefficient present in coefs0, -in an indeterminate order.

-

-

This can be optionally limited to coefficients whose -band number b and sample time t satisfy -b0b < b1 and -t0t < t1. -To process every coefficient present -in coefs0, pass INT_MIN, INT_MAX, INT64_MIN, INT64_MAX -for the arguments b0, b1, t0, -and t1, respectively. -

-

The function f should have the call signature

-
-
-template<class T>
-void f(int b, int64_t t, std::complex<T> &c0, std::complex<T> &ci...);
-
-

where

-
-
b
-
The band number of the frequency band the coefficients - c0 and ci... pertain to. - This may be either a bandpass band or the lowpass band.
-
t
-
The point in time the coefficients c0 and - ci... pertain to, in samples
-
c0
-
A reference to a complex coefficient from coefs0
-
ci...
-
Optional references to complex coefficients from the additional - coefficient sets coefsi....
-
-
- - -

The function f may read and/or modify each of the -coefficients passed through c0 and each -ci....

- -

The first coefficient set c0 is a special case when -it comes to the treatment of missing values. Coefficients missing -from c0 will not be iterated over at all, but when a -coefficient is iterated over and is missing from one of the additional -coefficient sets ci..., it will be automatically created -and initialized to zero in that additional coefficient set.

- -

Note: The template parameters C0 -and CI... exist to support the processing of coefficient -sets containing data of types other -than std::complex<T>, which is not currently part of the -documented API. In typical use, there is no need to specify them when -calling apply() because the template parameter list -can be deduced, but if they are expicitly specified, they should all -be std::complex<T>. -

- -

Creating New Coefficients

- -
-template <class T, class F, class C0, class... CI>
-void fill(F f,
-          int b0,
-          int b1,
-          int64_t t0,
-          int64_t t1,
-          coefs<T, C0> &coefs0,
-          coefs<T, CI>&... coefsi);
-
-

-Fill a region of the time-frequency plane with coefficients -and apply the function f to each. -

-

This works like process() except that it is not limited -to processing coefficients that already exist in coefs0; -instead, any missing coefficients in coefs0 as well as -any of the coefsi... are created and initialized to zero -before f is called.

- -

The t0 and t1 arguments must specify an -explicit, bounded time range — they must not be given as -INT64_MIN and/or INT64_MAX as that would mean creating coefficients -for an an astronomically large time range, requiring a correspondingly -astronomical amount of memory.

- -

Forgetting Coefficients

-
-template <class T>
-void forget_before(const analyzer<T> &a,
-                   coefs<T> &c,
-                   int64_t limit);
-
-

Allow the coefficients for points in time before limit -(a time in units of samples) to be forgotten. -Streaming applications can use this to free memory used by coefficients -that are no longer needed. Coefficients that have been forgotten will -read as zero. This does not guarantee that all coefficients before -limit are forgotten, only that ones for -limit or later are not, and that the amount of memory -consumed by any remaining coefficients before limit is -bounded.

- -

Legacy API For Iterating Over Existing Coefficients

- -

Prior to version 1.5, the only way to iterate over -coefficients was the apply() function. -It is similar to process(), except that it -

-
    -
  • requires an additional analyzer argument, -
  • takes arguments in a different order, -
  • applies a function f taking arguments in a different order, -
  • does not support restricting the processing to a range of band numbers, -
  • only supports iterating over a single coefficient set, and -
  • provides default values for t0 and t1. -
-

In new code, process() is preferred.

- -
-template <class T, class F>
-void apply(const analyzer<T> &a,
-           coefs<T> &c,
-           F f,
-           int64_t t0 = INT64_MIN,
-           int64_t t1 = INT64_MAX);
-
-

-Apply the function f to each coefficient in the coefficient -set c for points in time t that satisfy -t0t < t1. -If the t0 and t1 arguments are omitted, f -is applied to every coefficient. -

-
-
a
-
The spectrum analyzer that produced the coefficients c
-
c
-
A set of spectrogram coefficients
-
f
-
A function to apply to each coefficient in c, - with the call signature -
-template<class T>
-void f(std::complex<T> &coef, int band, int64_t t);
-
-
-
coef
-
A reference to a single complex coefficient. This may be read and/or modified.
-
band
-
The band number of the frequency band the coefficient coef0 pertains to. - This may be either a bandpass band or the lowpass band.
-
t
-
The point in time the coefficient c0 pertains to, in samples
-
t0
When not INT64_MIN, only apply f to the coefficients for time ≥ t0
-
t1
When not INT64_MAX, only apply f to the coefficients for time < t1
-
-
-
- - - - - diff --git a/lib/gaborator/doc/ref/intro.html b/lib/gaborator/doc/ref/intro.html deleted file mode 100644 index d45d775..0000000 --- a/lib/gaborator/doc/ref/intro.html +++ /dev/null @@ -1,41 +0,0 @@ - - - - - -Gaborator reference: API Introdution - - -

Gaborator reference: API Introduction

- -

The public API of the Gaborator library is defined in the HTML -documentation in the form of annotated C++ declarations. These are -similar to the actual declarations in the respective header files, but -simplified for clarity and omitting implementation details.

- -

The actual implementation in the header file may be different in a -number of ways but nonetheless compatible with the documented API. -For example, classes may be declared using the keyword struct -rather than class, function parameter names may be -different, types may be declared using different but equivalent -typedefs, and functions or templates in the header file may have -additional arguments with default values. Any classes, functions, and -other definitions not mentioned in the documentation should be -considered private and are subject to change or deletion without -notice. -

- -

All definitions are in the namespace gaborator. -Applications need to either prefix class names -with gaborator::, or use using namespace -gaborator;. -

- - - - - diff --git a/lib/gaborator/doc/ref/render_h.html b/lib/gaborator/doc/ref/render_h.html deleted file mode 100644 index 2ea3351..0000000 --- a/lib/gaborator/doc/ref/render_h.html +++ /dev/null @@ -1,74 +0,0 @@ - - - - - -Gaborator reference: render.h - - -

Gaborator reference: render.h

- -

Spectrogram Rendering with Power-of-Two Scaling

-
-template <class OI, class T>
-void render_p2scale(const analyzer<T> &a,
-                    const coefs<T> &c,
-                    int64_t xorigin, int64_t yorigin,
-                    int64_t xi0, int64_t xi1, int xe,
-                    int64_t yi0, int64_t yi1, int ye,
-                    OI output);
-
-

Render a rectangular array of pixel values representing signal -amplitudes in time-frequency space, optionally scaling up or -down by powers of two. -

-
-
a
-
The spectrum analyzer that produced the coefficients c
-
c
-
A set of spectrogram coefficients to render
-
xorigin
-
The point in time corresponding to pixel X coordinate 0, in samples
-
yorigin
-
The band number of the frequency band corresponding to pixel Y coordinate 0
-
xi0
-
The X coordinate of the first pixel to render
-
xi1
-
The X coordinate one past the last pixel to render
-
xe
-
The horizontal scaling exponent. One horizontal pixel corresponds to 2xe signal samples.
-
yi0
-
The Y coordinate of the first pixel to render
-
yi1
-
The Y coordinate one past the last pixel to render
-
ye
-
The vertical scaling exponent. One vertical pixel corresponds to 2ye frequency bands.
-
output
-
A random access iterator through which the output - pixel amplitude values will be written. This is - typically a float *. A total of - (xi1 - xi0) * (yi1 - yi0)) values will be written. -
-
- -

Utility Functions

-
-template <class T>
-unsigned int float2pixel_8bit(T amp);
-
-

Convert a normalized amplitude value to a 8-bit greyscale pixel value.

-
-
amp
-
A floating point value representing a signal amplitude, nominally ranging from 0 to 1
-
-

Returns an pixel value ranging from 0 to 255 (inclusive), using an -approximation of the sRGB gamma.

- - - - - diff --git a/lib/gaborator/doc/render.html b/lib/gaborator/doc/render.html deleted file mode 100644 index 981f954..0000000 --- a/lib/gaborator/doc/render.html +++ /dev/null @@ -1,395 +0,0 @@ - - - - - -Gaborator Example 1: Rendering a Spectrogram Image - - -

Example 1: Rendering a Spectrogram Image

- -

Introduction

- -

This example shows how to generate a greyscale constant-Q -spectrogram image from an audio file using the Gaborator library. -

- -

Preamble

- -

We start off with some boilerplate #includes.

- -
-#include <memory.h>
-#include <iostream>
-#include <fstream>
-#include <sndfile.h>
-
- -

The Gaborator is a header-only library — there are no C++ files -to compile, only header files to include. -The core spectrum analysis and resynthesis code is in -gaborator/gaborator.h, and the code for rendering -images from the spectrogram coefficients is in -gaborator/render.h.

- -
-#include <gaborator/gaborator.h>
-#include <gaborator/render.h>
-
- -

The program takes the names of the input audio file and output spectrogram -image file as command line arguments, so we check that they are present:

- -
-int main(int argc, char **argv) {
-    if (argc < 3) {
-        std::cerr << "usage: render input.wav output.pgm\n";
-        exit(1);
-    }
-
- -

Reading the Audio

- -

The audio file is read using the libsndfile library -and stored in a std::vector<float>. -Note that although libsndfile is used in this example, -the Gaborator library itself does not depend on or -use libsndfile.

-
-    SF_INFO sfinfo;
-    memset(&sfinfo, 0, sizeof(sfinfo));
-    SNDFILE *sf_in = sf_open(argv[1], SFM_READ, &sfinfo);
-    if (! sf_in) {
-        std::cerr << "could not open input audio file: "
-            << sf_strerror(sf_in) << "\n";
-        exit(1);
-    }
-    double fs = sfinfo.samplerate;
-    sf_count_t n_frames = sfinfo.frames;
-    sf_count_t n_samples = sfinfo.frames * sfinfo.channels;
-    std::vector<float> audio(n_samples);
-    sf_count_t n_read = sf_readf_float(sf_in, audio.data(), n_frames);
-    if (n_read != n_frames) {
-        std::cerr << "read error\n";
-        exit(1);
-    }
-    sf_close(sf_in);
-
-

In case the audio file is a stereo or multi-channel one, -mix down the channels to mono, into a new std::vector<float>: -

-    std::vector<float> mono(n_frames);
-    for (size_t i = 0; i < (size_t)n_frames; i++) {
-        float v = 0;
-        for (size_t c = 0; c < (size_t)sfinfo.channels; c++)
-            v += audio[i * sfinfo.channels + c];
-        mono[i] = v;
-    }
-
- -

The Spectrum Analysis Parameters

- -

Next, we need to choose some parameters for the spectrum analysis: -the frequency resolution, the frequency range, and optionally a -reference frequency.

- -

The frequency resolution is specified as a number of frequency -bands per octave. A typical number for analyzing music signals is 48 -bands per octave, or in other words, four bands per semitone -in the 12-note equal tempered scale.

- -

The frequency range is specified by giving a minimum frequency; -this is the lowest frequency that will be included in the spectrogram -display. -For audio signals, a typical minimum frequency is 20 Hz, -the lower limit of human hearing. In the Gaborator library, -all frequencies are given in units of the sample rate rather -than in Hz, so we need to divide the 20 Hz by the sample -rate of the input audio file: 20.0 / fs.

- -

Unlike the minimum frequency, the maximum frequency is not given -explicitly — instead, the analysis always produces coefficients -for frequencies all the way up to half the sample rate -(a.k.a. the Nyquist frequency). If you don't need the coefficients -for the highest frequencies, you can simply ignore them.

- -

If desired, one of the frequency bands can be exactly aligned with -a reference frequency. When analyzing music signals, this is -typically 440 Hz, the standard tuning of the note A4. -Like the minimum frequency, it is given in -units of the sample rate, so we pass 440.0 / fs.

- -

The parameters are held in an object of type -gaborator::parameters: -

-    gaborator::parameters params(48, 20.0 / fs, 440.0 / fs);
-
- -

The Spectrum Analyzer

- -

Next, we create an object of type gaborator::analyzer; -this is the workhorse that performs the actual spectrum analysis -(and/or resynthesis, but that's for a later example). -It is a template class, parametrized by the floating point type to -use for the calculations; this is typically float. -Constructing the gaborator::analyzer involves allocating and -precalculating all the filter coefficients and other auxiliary data needed -for the analysis and resynthesis, and this takes considerable time and memory, -so when analyzing multiple pieces of audio with the same -parameters, creating a single gaborator::analyzer -and reusing it is preferable to destroying and recreating it.

-
-    gaborator::analyzer<float> analyzer(params);
-
- -

The Spectrogram Coefficients

- -

The result of the spectrum analysis will be a set of spectrogram -coefficients. To store them, we will use a gaborator::coefs -object. Like the analyzer, this is a template class parametrized -by the data type. Because the layout of the coefficients is determined by -the spectrum analyzer, it must be passed as an argument to the constructor:

-
-    gaborator::coefs<float> coefs(analyzer);
-
- -

Running the Analysis

- -

Now we are ready to do the actual spectrum analysis, -by calling the analyze method of the spectrum -analyzer object. -The first argument to analyze is a float pointer -pointing to the first element in the array of samples to analyze. -The second and third arguments are of type -int64_t and indicate the time range covered by the -array, in units of samples. Since we are passing the whole file at -once, the beginning of the range is sample number zero, and the end is -sample number mono.size(). The fourth argument is a -reference to the set of coefficients that the results of the spectrum -analysis will be stored in. -

-
-    analyzer.analyze(mono.data(), 0, mono.size(), coefs);
-
- -

Rendering an Image

- -

Now there is a set of spectrogram coefficients in coefs. -To render them as an image, we will use the function -gaborator::render_p2scale(). -

- -

Rendering involves two different coordinate -spaces: the time-frequency coordinates of the spectrogram -coefficients, and the x-y coordinates of the image. -The two spaces are related by an origin and a scale factor, -each with an x and y component.

- -

The origin specifies the point in time-frequency space that -corresponds to the pixel coordinates (0, 0). Here, we will -use an origin where the x (time) component -is zero (the beginning of the signal), and the y (frequency) component -is the band number of the first (highest frequency) band:

- -
-    int64_t x_origin = 0;
-    int64_t y_origin = analyzer.bandpass_bands_begin();
-
- -

render_p2scale() supports scaling the spectrogram in -both the time (horizontal) and frequency (vertical) dimension, but only -by power-of-two scale factors. These scale factors are specified -relative to a reference scale of one vertical pixel per frequency band -and one horizontal pixel per signal sample. - -

Although a horizontal scale of one pixel per signal sample is a -mathematically pleasing reference point, this reference scale is not -used in practice because it would result in a spectrogram that is much -too stretched out horizontally. A more typical scale factor might be -210 = 1024, yielding one pixel for every 1024 signal -samples, which is about one pixel per 23 milliseconds of signal at a -sample rate of 44.1 kHz.

-
-    int x_scale_exp = 10;
-
- -

To ensure that the spectrogram will fit on the screen even in the -case of a long audio file, let's auto-scale it down further until -it is no more than 1000 pixels wide:

-
-    while ((n_frames >> x_scale_exp) > 1000)
-        x_scale_exp++;
-
- -

In the vertical, the reference scale factor of one pixel per -frequency band is reasonable, so we will use it as-is. In other words, -the vertical scale factor will be 20.

-
-    int y_scale_exp = 0;
-
- -

Next, we need to define the rectangular region of the image -coordinate space to render. Since we are rendering the entire -spectrogram rather than a tile, the top left corner of the -rectangle will have an origin of (0, 0). -

- -
-    int64_t x0 = 0;
-    int64_t y0 = 0;
-
- -

The coordinates of the bottom right corner are determined by the -length of the signal and the number of bands, respectively, taking the -scale factors into account. -The length of the signal in samples is n_frames, -and we get the number of bands as the difference of the end points of -the range of band numbers: -analyzer.bandpass_bands_end() - analyzer.bandpass_bands_begin(). -The scale factor is taken into account by right shifting by the -scale exponent. -

- -
-    int64_t x1 = n_frames >> x_scale_exp;
-    int64_t y1 = (analyzer.bandpass_bands_end() - analyzer.bandpass_bands_begin()) >> y_scale_exp;
-
- -

The right shift by y_scale_exp above doesn't actually -do anything because y_scale_exp is zero, but it would be -needed if, for example, you were to change y_scale_exp to -1 to get a spectrogram scaled to half the height. You could also make a -double-height spectrogram by setting y_scale_exp to -1, -but then you also need to change the ->> y_scale_exp to -<< -y_scale_exp since you can't shift by -a negative number. -

- -

We are now ready to render the spectrogram, producing -a vector of floating-point amplitude values, one per pixel. -Although this is stored as a 1-dimensional vector of floats, its -contents should be interpreted as a 2-dimensional rectangular array of -(y1 - y0) rows of (x1 - x0) columns -each, with the row indices increasing towards lower -frequencies and column indices increasing towards later -sampling times. -

-
-    std::vector<float> amplitudes((x1 - x0) * (y1 - y0));
-    gaborator::render_p2scale(
-        analyzer,
-        coefs,
-        x_origin, y_origin,
-        x0, x1, x_scale_exp,
-        y0, y1, y_scale_exp,
-        amplitudes.data());
-
- -

Writing the Image File

- -

To keep the code simple and to avoid additional library -dependencies, the image is stored in -pgm (Portable GreyMap) format, which is simple -enough to be generated with just a few lines of inline code. -Each amplitude value in amplitudes is converted into an 8-bit -gamma corrected pixel value by calling gaborator::float2pixel_8bit(). -To control the brightness of the resulting image, each -amplitude value is multiplied by a gain; this may have to be adjusted -depending on the type of signal and the amount of headroom in the -recording, but a gain of about 15 often works well for typical music -signals.

-
-    float gain = 15;
-    std::ofstream f;
-    f.open(argv[2], std::ios::out | std::ios::binary);
-    f << "P5\n" << (x1 - x0) << ' ' << (y1 - y0) << "\n255\n";
-    for (size_t i = 0; i < amplitudes.size(); i++)
-        f.put(gaborator::float2pixel_8bit(amplitudes[i] * gain));
-    f.close();
-
- -

Postamble

-

-To make the example code a complete program, -we just need to finish main(): -

-
-    return 0;
-}
-
- -

Compiling

-

-If you are using macOS, Linux, NetBSD, or a similar system, you can build -the example by running the following command in the examples -subdirectory. -You need to have libsndfile is installed and supported by -pkg-config. -

-
-c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` render.cc `pkg-config --libs sndfile` -o render
-
- -

Compiling for Speed

-

The above build command uses the Gaborator's built-in FFT implementation, -which is simple and portable but rather slow. Performance can be -significantly improved by using a faster FFT library. On macOS, you -can use the FFT from Apple's vDSP library by defining -GABORATOR_USE_VDSP and linking with the Accelerate -framework: -

-
-c++ -std=c++11 -I.. -O3 -ffast-math -DGABORATOR_USE_VDSP `pkg-config --cflags sndfile` render.cc `pkg-config --libs sndfile` -framework Accelerate -o render
-
- -

On Linux and NetBSD, you can use the PFFFT (Pretty Fast FFT) library. -You can get the latest version from -https://bitbucket.org/jpommier/pffft, -or the exact version that was used for testing from gaborator.com: -

- -
-wget http://download.gaborator.com/mirror/pffft/29e4f76ac53b.zip
-unzip 29e4f76ac53b.zip
-mv jpommier-pffft-29e4f76ac53b pffft
-
-

Then, compile it:

-
-cc -c -O3 -ffast-math pffft/pffft.c -o pffft/pffft.o
-
-

(If you are building for ARM, you will need to add -mfpu=neon to -both the above compilation command and the ones below.)

-

PFFFT is single precision only, but it comes with a copy of FFTPACK which can -be used for double-precision FFTs. Let's compile that, too:

-
-cc -c -O3 -ffast-math -DFFTPACK_DOUBLE_PRECISION pffft/fftpack.c -o pffft/fftpack.o
-
-

Then build the example and link it with both PFFFT and FFTPACK:

-
-c++ -std=c++11 -I.. -Ipffft -O3 -ffast-math -DGABORATOR_USE_PFFFT `pkg-config --cflags sndfile` render.cc pffft/pffft.o pffft/fftpack.o `pkg-config --libs sndfile` -o render
-
- -

Running

-

Running the following shell commands will download a short example -audio file (of picking each string on an acoustic guitar), generate -a spectrogram from it as a .pgm image, and then convert -the .pgm image into a JPEG image: -

-wget http://download.gaborator.com/audio/guitar.wav
-./render guitar.wav guitar.pgm
-cjpeg <guitar.pgm >guitar.jpg
-
- -

Example Output

-

The JPEG file produced by the above will look like this:

-Spectrogram - - - - - diff --git a/lib/gaborator/doc/snr.html b/lib/gaborator/doc/snr.html deleted file mode 100644 index 632c320..0000000 --- a/lib/gaborator/doc/snr.html +++ /dev/null @@ -1,121 +0,0 @@ - - - - - -Gaborator Example 4: Measuring the Signal-to-Noise Ratio - - -

Example 4: Measuring the Signal-to-Noise Ratio

- -

Introduction

- -

This example measures the signal-to-noise ratio (SNR) of the -resynthesis by analyzing and resynthesizing a test signal -and comparing the resynthesis result to the original. -

- -

Since it does not involve any audio file I/O, this example -does not require the sndfile library, making it the shortest -and simplest one by far.

- -

Preamble

- -
-#include <iostream>
-#include <iomanip>
-#include <random>
-#include <gaborator/gaborator.h>
-
- -

Amplitude Measurement

-

To calculate the signal-to-noise ratio, we need to measure the -amplitude of the orignal signal and the error residue. We will use -the root-mean-square amplitude, which is calculcated by the -function rms(). -

-
-double rms(const std::vector<float> &v) {
-    double sqsum = 0;
-    for (size_t i = 0; i < v.size(); i++) {
-        sqsum += v[i] * v[i];
-    }
-    return sqrt(sqsum);
-}
-
- -

Main Program

- -

For the test signal, we use a million samples of white noise with a -uniform amplitude distribution between -1 and +1.

-
-int main(int argc, char **argv) {
-    size_t len = 1000000;
-    std::vector<float> signal_in(len);
-    std::minstd_rand rand;
-    std::uniform_real_distribution<> uniform(-1.0, 1.0);
-    for (size_t i = 0; i < len; i++)
-        signal_in[i] = uniform(rand);
-
-

Then we create a spectrum analyzer with 48 bands per octave -and a frequency range of 3 decades (0.0005 to 0.5 times the sample rate):

-
-    gaborator::parameters params(48, 5e-4);
-    gaborator::analyzer<float> analyzer(params);
-
-

...and run the spectrum analyzis:

-
-    gaborator::coefs<float> coefs(analyzer);
-    analyzer.analyze(signal_in.data(), 0, len, coefs);
-
-

...resynthesize the signal into signal_out: -

-    std::vector<float> signal_out(len);
-    analyzer.synthesize(coefs, 0, len, signal_out.data());
-
-

...measure the resynthesis error:

-
-    std::vector<float> error(len);
-    for (size_t i = 0; i < len; i++)
-         error[i] = signal_out[i] - signal_in[i];
-
-

...calculate the signal-to-noise ratio:

-
-    double snr = rms(signal_in) / rms(error);
-
-

...and print it in decibels:

-
-    std::cout << std::fixed << std::setprecision(1) << 20 * log10(snr) << " dB\n";
-}
-
-

Compiling

-

Like Example 1, this example -can be built using a one-line build command: -

-
-c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` snr.cc `pkg-config --libs sndfile` -o snr
-
-

Or using the vDSP FFT on macOS:

-
-c++ -std=c++11 -I.. -O3 -ffast-math -DGABORATOR_USE_VDSP `pkg-config --cflags sndfile` snr.cc `pkg-config --libs sndfile` -framework Accelerate -o snr
-
-

Or using PFFFT (see Example 1 for how to download and build PFFFT):

-
-c++ -std=c++11 -I.. -Ipffft -O3 -ffast-math -DGABORATOR_USE_PFFFT `pkg-config --cflags sndfile` snr.cc pffft/pffft.o pffft/fftpack.o `pkg-config --libs sndfile` -o snr
-
- -

Running

-

The program is run with no arguments:

-
-./snr
-
-

This will print the SNR which should be more than 100 dB if the library is working correctly.

- - - - - diff --git a/lib/gaborator/doc/spectrogram.jpg b/lib/gaborator/doc/spectrogram.jpg deleted file mode 100644 index f854519..0000000 Binary files a/lib/gaborator/doc/spectrogram.jpg and /dev/null differ diff --git a/lib/gaborator/doc/stream.html b/lib/gaborator/doc/stream.html deleted file mode 100644 index f270ecc..0000000 --- a/lib/gaborator/doc/stream.html +++ /dev/null @@ -1,279 +0,0 @@ - - - - - -Gaborator Example 3: Streaming - - -

Example 3: Streaming

- -

Introduction

- -

This example shows how to process streaming audio a block at a time, -rather than operating on a complete recording at once as in the previous -examples.

- -

This program doesn't do anything particulary useful — it just -inverts the phase of the signal, but not using the obvious method of -changing the sign of each sample but by changing the sign of each -spectrogram coefficient. Consider the expression coef = --coef a placeholder for your own streaming filter or effect -code.

- -

Preamble

- -
-#include <memory.h>
-#include <iostream>
-#include <sndfile.h>
-#include <gaborator/gaborator.h>
-
-int main(int argc, char **argv) {
-    if (argc < 3) {
-        std::cerr << "usage: stream input.wav output.wav\n";
-        exit(1);
-    }
-
- -

Opening the Streams

- -

We again use libsndfile to read the input and write the output. -To keep it simple, this example only handles mono files.

-
-    SF_INFO sfinfo;
-    memset(&sfinfo, 0, sizeof(sfinfo));
-    SNDFILE *sf_in = sf_open(argv[1], SFM_READ, &sfinfo);
-    if (! sf_in) {
-        std::cerr << "could not open input audio file: "
-            << sf_strerror(sf_in) << "\n";
-        exit(1);
-    }
-    if (sfinfo.channels != 1) {
-        std::cerr << "only mono files are supported\n";
-        exit(1);
-    }
-    double fs = sfinfo.samplerate;
-
-    SNDFILE *sf_out = sf_open(argv[2], SFM_WRITE, &sfinfo);
-    if (! sf_out) {
-        std::cerr << "could not open output audio file: "
-            << sf_strerror(sf_out) << "\n";
-        exit(1);
-    }
-
-

The next couple of lines work around a design flaw in -libsndfile. By default, when reading a 16-bit -audio file as floating point data and then writing them -as another 16-bit audio file, libsndfile will use slightly -different scale factors on input and output, and the output will -not be bit-identical to the input. To make it easier to verify -that this example actually yields correct results to within -the full 16-bit precision, we select a non-normalized floating -point representation, which does not suffer from this flaw.

-
-    sf_command(sf_in, SFC_SET_NORM_FLOAT, NULL, SF_FALSE);
-    sf_command(sf_out, SFC_SET_NORM_FLOAT, NULL, SF_FALSE);
-
- -

Spectrum Analysis Parameters

- -

As in Example 1, the parameters are chosen for analyzing music, but -to reduce the latency, the number of frequency bands per octave is reduced -from 48 to 12 (one per semitone), and the lower frequency limit of -the bandpass filter bank is raised from 20 Hz to 200 Hz.

-
-    gaborator::parameters params(12, 200.0 / fs, 440.0 / fs);
-    gaborator::analyzer<float> analyzer(params);
-
- -

Calculating Latency

-

The spectrogram coefficients are calculated by applying symmetric -FIR filters to the audio signal. This means a spectrogram coefficient -for any given point in time t is a weighted average of samples -from both before and after t, representing both past and future -signal. The width of the filter impulse response depends on the -bandwidth, which in turn depends on the center frequency of its band. -The lowest-frequency filters have the narrowest bandwidths, and -therefore the widest impulses response, and need the greatest amount -of past and future signal. The width of the filter impulse response -is called its support, and the worst-case (widest) support of -any analysis filter can be found by calling the function -gaborator::analyzer::analysis_support(). This returns -the one-sided support, the width of the impulse -response to each side of its center, as a floating point number. -To be on the safe side, let's round this up to the next integer:

-
-    size_t analysis_support = ceil(analyzer.analysis_support());
-
-

Similarly, when resynthesizing audio from coefficients, calculating -a sample at time t involves applying symmetric FIR -reconstruction filters, calculating a weighted average of both past and -future spectrogram coefficients. The support of the widest reconstruction -filter can be calculated by calling -gaborator::analyzer::synthesis_support(): -

-
-    size_t synthesis_support = ceil(analyzer.synthesis_support());
-
- -

In a real-time application, the need to access future signal -samples and/or coefficients causes latency. A real-time audio -analysis application that needs to examine the coefficients for -time t can only do so when it has received the input samples up -to time t + analysis_support, and therefore has a minimum latency of -analysis_support. A real-time filtering or effect -application, such as the present example, -incurs latency from both analysis and reconstruction -filters, and can only produce the output sample for time t once -it has received the input samples up to -t + analysis_support + synthesis_support, -for a minimum latency of analysis_support + synthesis_support. -Let's print this total latency to standard output: -

-
-    std::cerr << "latency: " << analysis_support + synthesis_support << " samples\n";
-
- -

In a practical real-time system, there will be additional latency -caused by processing the signal in blocks of samples rather than a -sample at a time. Since the block size is a property of the overall -system, and causes latency even if the Gaborator is not involved, that -latency is considered outside the scope of this discussion. -

- -

Streaming

-

To mimic a typical real-time system, the audio is processed -in fixed-size blocks (here, 1024 samples). If the size -of the input file is not divisible by the block size, the last block -is padded with zeroes. -The variable t_in keeps track of time, indicating -the sampling time of the first sample of the current input block, -in units of samples. -

-
-    gaborator::coefs<float> coefs(analyzer);
-    const size_t blocksize = 1024;
-    std::vector<float> buf(blocksize);
-    int64_t t_in = 0;
-    for (;;) {
-        sf_count_t n_read = sf_readf_float(sf_in, buf.data(), blocksize);
-        if (n_read == 0)
-            break;
-        if (n_read < blocksize)
-            std::fill(buf.data() + n_read, buf.data() + blocksize, 0);
-
-

Now we can spectrum analyze the current block of samples. Note how -the time range, -t_in...t_in + blocksize, -is explicitly passed to analyze(). -

-
-        analyzer.analyze(buf.data(), t_in, t_in + blocksize, coefs);
-
-

The call to analyze() updates the coefficients -for the time range from t_in - analysis_support to -t_in + blocksize + analysis_support. The oldest -blocksize samples of this time range, -that is, from t_in - analysis_support to -t_in - analysis_support + blocksize, were now updated for -the last time and will not be affected by future input blocks. -Therefore, it is now safe to examine and/or modify these -coefficients as required by your application. Here, by way -of example, we simply change their signs to invert the phase of the signal. -Note that unlike the earlier filter example where prorcess() -applied a function to all the coefficients, here it is applied only to -the coefficients within a limited time range. -

-
-        process(
-            [&](int, int64_t, std::complex<float> &coef) {
-                 coef = -coef;
-            },
-            INT_MIN, INT_MAX,
-            t_in - (int)analysis_support,
-            t_in - (int)analysis_support + (int)blocksize,
-            coefs);
-
- -

Next, we will generate a block of output samples. To get correct results, -we can only generate output when the coefficients that the output samples -depend on will no longer change. Specifically, a resynthesized audio -sample for time t will depend on the coefficients of the -time range t - synthesis_support...t + -synthesis_support. To ensure that the resynthesis uses only -coefficients that have already been processed by -the process() call above, the most recent block of samples -that can safely be resynthesized ranges from t_out = t_in - -analysis_support - synthesis_support to t_out + -blocksize.

-
-        int64_t t_out = t_in - analysis_support - synthesis_support;
-        analyzer.synthesize(coefs, t_out, t_out + blocksize, buf.data());
-
-

The synthesized audio can now be written to the output file:

-
-        sf_count_t n_written = sf_writef_float(sf_out, buf.data(), blocksize);
-        if (n_written != blocksize) {
-            std::cerr << "write error\n";
-            exit(1);
-        }
-
-

Coefficients older than t_out + blocksize - synthesis_support -will no longer be needed to synthesize the next block of output signal, so -it's now OK to forget them and free the memory they used: -

-
-        forget_before(analyzer, coefs, t_out + blocksize - synthesis_support);
-
-

This concludes the block-by-block processing loop.

-
-        t_in += blocksize;
-    }
-
- -

Postamble

-
-    sf_close(sf_in);
-    sf_close(sf_out);
-    return 0;
-}
-
- -

Compiling

-

Like the previous ones, this example can also be built using a one-line build command: -

-
-c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` stream.cc `pkg-config --libs sndfile` -o stream
-
-

Or using the vDSP FFT on macOS:

-
-c++ -std=c++11 -I.. -O3 -ffast-math -DGABORATOR_USE_VDSP `pkg-config --cflags sndfile` stream.cc `pkg-config --libs sndfile` -framework Accelerate -o stream
-
-

Or using PFFFT (see Example 1 for how to download and build PFFFT):

-
-c++ -std=c++11 -I.. -Ipffft -O3 -ffast-math -DGABORATOR_USE_PFFFT `pkg-config --cflags sndfile` stream.cc pffft/pffft.o pffft/fftpack.o `pkg-config --libs sndfile` -o stream
-
- -

Running

-

Running the following shell commands will download an example -audio file containing an impulse (a single sample of maximum amplitude) -padded with silence to a total of 65536 samples, and process it.

-
-wget http://download.gaborator.com/audio/impulse.wav
-./stream impulse.wav impulse_streamed.wav
-
- -

The file impulse_streamed.wav will be identical to -impulse.wav except that the impulse will be of -opposite polarity, and delayed by the latency of -analysis_support + synthesis_support samples.

- - - - - diff --git a/lib/gaborator/doc/synth.html b/lib/gaborator/doc/synth.html deleted file mode 100644 index edfa6f2..0000000 --- a/lib/gaborator/doc/synth.html +++ /dev/null @@ -1,213 +0,0 @@ - - - - - -Gaborator Example 5: Synthesis from Scratch - - -

Example 5: Synthesis from Scratch

- -

Introduction

- -

This example demonstrates how to synthesize a signal by creating -spectrogram coefficients from scratch rather than by analyzing an -existing signal. It creates a random pentatonic melody of decaying -sine waves as spectrogram coefficients and then synthesizes audio -from them. -

- -

Preamble

- -

This example program takes a single command line argument, the name -of the output file.

-
-#include <memory.h>
-#include <iostream>
-#include <sndfile.h>
-#include <gaborator/gaborator.h>
-
-int main(int argc, char **argv) {
-    if (argc < 2) {
-        std::cerr << "usage: synth output.wav\n";
-        exit(1);
-    }
-
- -

Synthesis Parameters

- -

Although this example does not perform any analysis, we nonetheless -need to create an analyzer object, as it is used for both -analysis and synthesis purposes. To generate the frequencies of the -12-note equal-tempered scale, we need 12 bands per octave; a multiple -of 12 would also work, but here we don't need the added frequency -resolution that would bring, and the time resolution would be -worse.

- -

To simplify converting MIDI note numbers to band numbers, we choose -the frequency of MIDI note 0 as the reference frequency; this is -8.18 Hz, which happens to be outside the frequency range of the -bandpass filter bank, but that doesn't matter.

- -
-    double fs = 44100;
-    gaborator::parameters params(12, 20.0 / fs, 8.18 / fs);
-    gaborator::analyzer<float> analyzer(params);
-
- -

Melody Parameters

- -

-We will use the A minor pentatonic scale, which contains the -following notes (using the MIDI note numbering):

-
-    static int pentatonic[] = { 57, 60, 62, 64, 67 };
-
- -

-The melody will consist of 64 notes, at a tempo of 120 beats per -minute: -

-
-    int n_notes = 64;
-    double tempo = 120.0;
-    double beat_duration = 60.0 / tempo;
-
- -

-The variable volume determines the amplitude of -each note, and has been chosen such that there will be no clipping -of the final output. -

-
-    float volume = 0.2;
-
- -

Composition

- -

We start with an empty coefficient set:

-
-    gaborator::coefs<float> coefs(analyzer);
-
- -

Each note is chosen randomly from the pentatonic scale and added -to the coefficient set by calling the function fill(). -The fill() function is similar to the process() -function used in previous examples, except that it can be used to -create new coefficients rather than just modifying existing ones.

- -

Each note is created by calling fill() on a region of -the time-frequency plane that covers a single band in the frequency -dimension and the duration of the note in the time dimension. Each -coefficient within this region is set to a complex number whose -magnitude decays exponentially over time, like the amplitude of a -plucked string. The phase is arbitrarily set to zero by using an -imaginary part of zero. Since notes can overlap, the new coefficients -are added to any existing ones using the += operator -rather than overwriting them.

- -

Note that band numbers increase towards lower frequencies but MIDI -note numbers increase towards higher frequencies, hence the minus sign -in front of midi_note. -

- -
-    for (int i = 0; i < n_notes; i++) {
-        int midi_note = pentatonic[rand() % 5];
-        double note_start_time = beat_duration * i;
-        double note_end_time = note_start_time + 3.0;
-        int band = analyzer.band_ref() - midi_note;
-        fill([&](int, int64_t t, std::complex<float> &coef) {
-                float amplitude =
-                    volume * expf(-2.0f * (float)(t / fs - note_start_time));
-                coef += std::complex<float>(amplitude, 0.0f);
-            },
-            band, band + 1,
-            note_start_time * fs, note_end_time * fs,
-            coefs);
-    }
-
- -

Synthesis

- -

We can now synthesize audio from the coefficients by -calling synthesize(). Audio will be generated -starting half a second before the first note to allow for the pre-ringing -of the synthesis filter, and ending a few seconds after the -last note to allow for its decay. -

-
-    double audio_start_time = -0.5;
-    double audio_end_time = beat_duration * n_notes + 5.0;
-    int64_t start_frame = audio_start_time * fs;
-    int64_t end_frame = audio_end_time * fs;
-    size_t n_frames = end_frame - start_frame;
-    std::vector<float> audio(n_frames);
-    analyzer.synthesize(coefs, start_frame, end_frame, audio.data());
-
- -

Writing the Audio

- -

Since there is no input audio file to inherit a file format from, -we need to choose a file format for the output file by filling in the -sfinfo structure:

-
-    SF_INFO sfinfo;
-    memset(&sfinfo, 0, sizeof(sfinfo));
-    sfinfo.samplerate = fs;
-    sfinfo.channels = 1;
-    sfinfo.format = SF_FORMAT_WAV | SF_FORMAT_PCM_16;
-
- -

The rest is identical to -Example 2: -

-
-    SNDFILE *sf_out = sf_open(argv[1], SFM_WRITE, &sfinfo);
-    if (! sf_out) {
-        std::cerr << "could not open output audio file: "
-            << sf_strerror(sf_out) << "\n";
-        exit(1);
-    }
-    sf_command(sf_out, SFC_SET_CLIPPING, NULL, SF_TRUE);
-    sf_count_t n_written = sf_writef_float(sf_out, audio.data(), n_frames);
-    if (n_written != n_frames) {
-        std::cerr << "write error\n";
-        exit(1);
-    }
-    sf_close(sf_out);
-    return 0;
-}
-
- -

Compiling

-

Like Example 1, this example -can be built using a one-line build command: -

-
-c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` synth.cc `pkg-config --libs sndfile` -o synth
-
-

Or using the vDSP FFT on macOS:

-
-c++ -std=c++11 -I.. -O3 -ffast-math -DGABORATOR_USE_VDSP `pkg-config --cflags sndfile` synth.cc `pkg-config --libs sndfile` -framework Accelerate -o synth
-
-

Or using PFFFT (see Example 1 for how to download and build PFFFT):

-
-c++ -std=c++11 -I.. -Ipffft -O3 -ffast-math -DGABORATOR_USE_PFFFT `pkg-config --cflags sndfile` synth.cc pffft/pffft.o pffft/fftpack.o `pkg-config --libs sndfile` -o synth
-
- -

Running

-

The example program can be run using the command

-
-./synth melody.wav
-
-

The resulting audio will be in melody.wav.

- - - - - diff --git a/lib/gaborator/examples/filter.cc b/lib/gaborator/examples/filter.cc deleted file mode 100644 index ad3891d..0000000 --- a/lib/gaborator/examples/filter.cc +++ /dev/null @@ -1,69 +0,0 @@ -// See ../doc/filter.html for commentary - -#include -#include -#include -#include - -int main(int argc, char **argv) { - if (argc < 3) { - std::cerr << "usage: filter input.wav output.wav\n"; - exit(1); - } - SF_INFO sfinfo; - memset(&sfinfo, 0, sizeof(sfinfo)); - SNDFILE *sf_in = sf_open(argv[1], SFM_READ, &sfinfo); - if (! sf_in) { - std::cerr << "could not open input audio file: " - << sf_strerror(sf_in) << "\n"; - exit(1); - } - double fs = sfinfo.samplerate; - sf_count_t n_frames = sfinfo.frames; - sf_count_t n_samples = sfinfo.frames * sfinfo.channels; - std::vector audio(n_samples); - sf_count_t n_read = sf_readf_float(sf_in, audio.data(), n_frames); - if (n_read != n_frames) { - std::cerr << "read error\n"; - exit(1); - } - sf_close(sf_in); - gaborator::parameters params(100, 20.0 / fs); - gaborator::analyzer analyzer(params); - std::vector band_gains(analyzer.bands_end()); - for (int band = analyzer.bandpass_bands_begin(); band < analyzer.bandpass_bands_end(); band++) { - double f_hz = analyzer.band_ff(band) * fs; - band_gains[band] = 1.0 / sqrt(f_hz / 20.0); - } - band_gains[analyzer.band_lowpass()] = band_gains[analyzer.bandpass_bands_end() - 1]; - for (sf_count_t ch = 0; ch < sfinfo.channels; ch++) { - std::vector channel(n_frames); - for (sf_count_t i = 0; i < n_frames; i++) - channel[i] = audio[i * sfinfo.channels + ch]; - gaborator::coefs coefs(analyzer); - analyzer.analyze(channel.data(), 0, channel.size(), coefs); - process([&](int band, int64_t, std::complex &coef) { - coef *= band_gains[band]; - }, - INT_MIN, INT_MAX, - INT64_MIN, INT64_MAX, - coefs); - analyzer.synthesize(coefs, 0, channel.size(), channel.data()); - for (sf_count_t i = 0; i < n_frames; i++) - audio[i * sfinfo.channels + ch] = channel[i]; - } - SNDFILE *sf_out = sf_open(argv[2], SFM_WRITE, &sfinfo); - if (! sf_out) { - std::cerr << "could not open output audio file: " - << sf_strerror(sf_out) << "\n"; - exit(1); - } - sf_command(sf_out, SFC_SET_CLIPPING, NULL, SF_TRUE); - sf_count_t n_written = sf_writef_float(sf_out, audio.data(), n_frames); - if (n_written != n_frames) { - std::cerr << "write error\n"; - exit(1); - } - sf_close(sf_out); - return 0; -} diff --git a/lib/gaborator/examples/render.cc b/lib/gaborator/examples/render.cc deleted file mode 100644 index 72cd7d3..0000000 --- a/lib/gaborator/examples/render.cc +++ /dev/null @@ -1,69 +0,0 @@ -// See ../doc/render.html for commentary - -#include -#include -#include -#include -#include -#include -int main(int argc, char **argv) { - if (argc < 3) { - std::cerr << "usage: render input.wav output.pgm\n"; - exit(1); - } - SF_INFO sfinfo; - memset(&sfinfo, 0, sizeof(sfinfo)); - SNDFILE *sf_in = sf_open(argv[1], SFM_READ, &sfinfo); - if (! sf_in) { - std::cerr << "could not open input audio file: " - << sf_strerror(sf_in) << "\n"; - exit(1); - } - double fs = sfinfo.samplerate; - sf_count_t n_frames = sfinfo.frames; - sf_count_t n_samples = sfinfo.frames * sfinfo.channels; - std::vector audio(n_samples); - sf_count_t n_read = sf_readf_float(sf_in, audio.data(), n_frames); - if (n_read != n_frames) { - std::cerr << "read error\n"; - exit(1); - } - sf_close(sf_in); - std::vector mono(n_frames); - for (size_t i = 0; i < (size_t)n_frames; i++) { - float v = 0; - for (size_t c = 0; c < (size_t)sfinfo.channels; c++) - v += audio[i * sfinfo.channels + c]; - mono[i] = v; - } - gaborator::parameters params(48, 20.0 / fs, 440.0 / fs); - gaborator::analyzer analyzer(params); - gaborator::coefs coefs(analyzer); - analyzer.analyze(mono.data(), 0, mono.size(), coefs); - int64_t x_origin = 0; - int64_t y_origin = analyzer.bandpass_bands_begin(); - int x_scale_exp = 10; - while ((n_frames >> x_scale_exp) > 1000) - x_scale_exp++; - int y_scale_exp = 0; - int64_t x0 = 0; - int64_t y0 = 0; - int64_t x1 = n_frames >> x_scale_exp; - int64_t y1 = (analyzer.bandpass_bands_end() - analyzer.bandpass_bands_begin()) >> y_scale_exp; - std::vector amplitudes((x1 - x0) * (y1 - y0)); - gaborator::render_p2scale( - analyzer, - coefs, - x_origin, y_origin, - x0, x1, x_scale_exp, - y0, y1, y_scale_exp, - amplitudes.data()); - float gain = 15; - std::ofstream f; - f.open(argv[2], std::ios::out | std::ios::binary); - f << "P5\n" << (x1 - x0) << ' ' << (y1 - y0) << "\n255\n"; - for (size_t i = 0; i < amplitudes.size(); i++) - f.put(gaborator::float2pixel_8bit(amplitudes[i] * gain)); - f.close(); - return 0; -} diff --git a/lib/gaborator/examples/snr.cc b/lib/gaborator/examples/snr.cc deleted file mode 100644 index b4d65d8..0000000 --- a/lib/gaborator/examples/snr.cc +++ /dev/null @@ -1,32 +0,0 @@ -// See ../doc/snr.html for commentary - -#include -#include -#include -#include -double rms(const std::vector &v) { - double sqsum = 0; - for (size_t i = 0; i < v.size(); i++) { - sqsum += v[i] * v[i]; - } - return sqrt(sqsum); -} -int main(int argc, char **argv) { - size_t len = 1000000; - std::vector signal_in(len); - std::minstd_rand rand; - std::uniform_real_distribution<> uniform(-1.0, 1.0); - for (size_t i = 0; i < len; i++) - signal_in[i] = uniform(rand); - gaborator::parameters params(48, 5e-4); - gaborator::analyzer analyzer(params); - gaborator::coefs coefs(analyzer); - analyzer.analyze(signal_in.data(), 0, len, coefs); - std::vector signal_out(len); - analyzer.synthesize(coefs, 0, len, signal_out.data()); - std::vector error(len); - for (size_t i = 0; i < len; i++) - error[i] = signal_out[i] - signal_in[i]; - double snr = rms(signal_in) / rms(error); - std::cout << std::fixed << std::setprecision(1) << 20 * log10(snr) << " dB\n"; -} diff --git a/lib/gaborator/examples/stream.cc b/lib/gaborator/examples/stream.cc deleted file mode 100644 index fb2da1b..0000000 --- a/lib/gaborator/examples/stream.cc +++ /dev/null @@ -1,72 +0,0 @@ -// See ../doc/stream.html for commentary - -#include -#include -#include -#include - -int main(int argc, char **argv) { - if (argc < 3) { - std::cerr << "usage: stream input.wav output.wav\n"; - exit(1); - } - SF_INFO sfinfo; - memset(&sfinfo, 0, sizeof(sfinfo)); - SNDFILE *sf_in = sf_open(argv[1], SFM_READ, &sfinfo); - if (! sf_in) { - std::cerr << "could not open input audio file: " - << sf_strerror(sf_in) << "\n"; - exit(1); - } - if (sfinfo.channels != 1) { - std::cerr << "only mono files are supported\n"; - exit(1); - } - double fs = sfinfo.samplerate; - - SNDFILE *sf_out = sf_open(argv[2], SFM_WRITE, &sfinfo); - if (! sf_out) { - std::cerr << "could not open output audio file: " - << sf_strerror(sf_out) << "\n"; - exit(1); - } - sf_command(sf_in, SFC_SET_NORM_FLOAT, NULL, SF_FALSE); - sf_command(sf_out, SFC_SET_NORM_FLOAT, NULL, SF_FALSE); - gaborator::parameters params(12, 200.0 / fs, 440.0 / fs); - gaborator::analyzer analyzer(params); - size_t analysis_support = ceil(analyzer.analysis_support()); - size_t synthesis_support = ceil(analyzer.synthesis_support()); - std::cerr << "latency: " << analysis_support + synthesis_support << " samples\n"; - gaborator::coefs coefs(analyzer); - const size_t blocksize = 1024; - std::vector buf(blocksize); - int64_t t_in = 0; - for (;;) { - sf_count_t n_read = sf_readf_float(sf_in, buf.data(), blocksize); - if (n_read == 0) - break; - if (n_read < blocksize) - std::fill(buf.data() + n_read, buf.data() + blocksize, 0); - analyzer.analyze(buf.data(), t_in, t_in + blocksize, coefs); - process( - [&](int, int64_t, std::complex &coef) { - coef = -coef; - }, - INT_MIN, INT_MAX, - t_in - (int)analysis_support, - t_in - (int)analysis_support + (int)blocksize, - coefs); - int64_t t_out = t_in - analysis_support - synthesis_support; - analyzer.synthesize(coefs, t_out, t_out + blocksize, buf.data()); - sf_count_t n_written = sf_writef_float(sf_out, buf.data(), blocksize); - if (n_written != blocksize) { - std::cerr << "write error\n"; - exit(1); - } - forget_before(analyzer, coefs, t_out + blocksize - synthesis_support); - t_in += blocksize; - } - sf_close(sf_in); - sf_close(sf_out); - return 0; -} diff --git a/lib/gaborator/examples/synth.cc b/lib/gaborator/examples/synth.cc deleted file mode 100644 index b5aff70..0000000 --- a/lib/gaborator/examples/synth.cc +++ /dev/null @@ -1,62 +0,0 @@ -// See ../doc/synth.html for commentary - -#include -#include -#include -#include - -int main(int argc, char **argv) { - if (argc < 2) { - std::cerr << "usage: synth output.wav\n"; - exit(1); - } - double fs = 44100; - gaborator::parameters params(12, 20.0 / fs, 8.18 / fs); - gaborator::analyzer analyzer(params); - static int pentatonic[] = { 57, 60, 62, 64, 67 }; - int n_notes = 64; - double tempo = 120.0; - double beat_duration = 60.0 / tempo; - float volume = 0.2; - gaborator::coefs coefs(analyzer); - for (int i = 0; i < n_notes; i++) { - int midi_note = pentatonic[rand() % 5]; - double note_start_time = beat_duration * i; - double note_end_time = note_start_time + 3.0; - int band = analyzer.band_ref() - midi_note; - fill([&](int, int64_t t, std::complex &coef) { - float amplitude = - volume * expf(-2.0f * (float)(t / fs - note_start_time)); - coef += std::complex(amplitude, 0.0f); - }, - band, band + 1, - note_start_time * fs, note_end_time * fs, - coefs); - } - double audio_start_time = -0.5; - double audio_end_time = beat_duration * n_notes + 5.0; - int64_t start_frame = audio_start_time * fs; - int64_t end_frame = audio_end_time * fs; - size_t n_frames = end_frame - start_frame; - std::vector audio(n_frames); - analyzer.synthesize(coefs, start_frame, end_frame, audio.data()); - SF_INFO sfinfo; - memset(&sfinfo, 0, sizeof(sfinfo)); - sfinfo.samplerate = fs; - sfinfo.channels = 1; - sfinfo.format = SF_FORMAT_WAV | SF_FORMAT_PCM_16; - SNDFILE *sf_out = sf_open(argv[1], SFM_WRITE, &sfinfo); - if (! sf_out) { - std::cerr << "could not open output audio file: " - << sf_strerror(sf_out) << "\n"; - exit(1); - } - sf_command(sf_out, SFC_SET_CLIPPING, NULL, SF_TRUE); - sf_count_t n_written = sf_writef_float(sf_out, audio.data(), n_frames); - if (n_written != n_frames) { - std::cerr << "write error\n"; - exit(1); - } - sf_close(sf_out); - return 0; -} diff --git a/lib/gaborator/gaborator/affine_transform.h b/lib/gaborator/gaborator/affine_transform.h deleted file mode 100644 index c2ea997..0000000 --- a/lib/gaborator/gaborator/affine_transform.h +++ /dev/null @@ -1,42 +0,0 @@ -// -// A class for affine transforms (ax + b) of scalar values -// -// Copyright (C) 2020-2021 Andreas Gustafsson. This file is part of -// the Gaborator library source distribution. See the file LICENSE at -// the top level of the distribution for license information. -// - -#ifndef _GABORATOR_AFFINE_TRANSFORM_H -#define _GABORATOR_AFFINE_TRANSFORM_H - -namespace gaborator { - -struct affine_transform { - affine_transform(): a(0), b(0) { } - affine_transform(double a_, double b_): a(a_), b(b_) { } - affine_transform(const affine_transform &rhs): a(rhs.a), b(rhs.b) { } - double operator()(double x) const { return a * x + b; } - affine_transform inverse() const { - return affine_transform(1.0 / a, -b / a); - } - static affine_transform identity() { return affine_transform(1, 0); } - double a, b; -}; - -// Composition - -static inline affine_transform -operator *(const affine_transform &a, const affine_transform &b) { - return affine_transform(a.a * b.a, a.a * b.b + a.b); -} - -// Equality - -static inline bool -operator ==(const affine_transform &a, const affine_transform &b) { - return a.a == b.a && a.b == b.b; -} - -} // namespace - -#endif diff --git a/lib/gaborator/gaborator/fft.h b/lib/gaborator/gaborator/fft.h deleted file mode 100644 index 66921e9..0000000 --- a/lib/gaborator/gaborator/fft.h +++ /dev/null @@ -1,29 +0,0 @@ -// -// Fast Fourier transform -// -// Copyright (C) 2016-2021 Andreas Gustafsson. This file is part of -// the Gaborator library source distribution. See the file LICENSE at -// the top level of the distribution for license information. -// - -#ifndef _GABORATOR_FFT_H -#define _GABORATOR_FFT_H - -#include "gaborator/fft_naive.h" - -#if GABORATOR_USE_VDSP -#include "gaborator/fft_vdsp.h" -#define GABORATOR_USE_REAL_FFT 1 -#define GABORATOR_MIN_FFT_SIZE 1 -#elif GABORATOR_USE_PFFFT -#include "gaborator/fft_pffft.h" -#define GABORATOR_USE_REAL_FFT 1 -#define GABORATOR_MIN_FFT_SIZE 32 -#else -// Use the naive FFT -// Do not define GABORATOR_USE_REAL_FFT as it is slower than -// using the complex code. -#define GABORATOR_MIN_FFT_SIZE 1 -#endif - -#endif diff --git a/lib/gaborator/gaborator/fft_naive.h b/lib/gaborator/gaborator/fft_naive.h deleted file mode 100644 index 875d99e..0000000 --- a/lib/gaborator/gaborator/fft_naive.h +++ /dev/null @@ -1,191 +0,0 @@ -// -// Fast Fourier transform, naive reference implementations -// -// Copyright (C) 1992-2021 Andreas Gustafsson. This file is part of -// the Gaborator library source distribution. See the file LICENSE at -// the top level of the distribution for license information. -// - -// Based on the module "fft" used in audsl/test, audsl/mls, -// scope/core, whitesig - -#ifndef _GABORATOR_FFT_NAIVE_H -#define _GABORATOR_FFT_NAIVE_H - -#include -#include -#include -#include - -#include - -namespace gaborator { - -template -struct fft { - typedef typename std::iterator_traits::value_type C; // complex - typedef typename C::value_type T; // float/double - typedef typename std::vector twiddle_vector; - - fft(unsigned int n): n_(n), wtab(n / 2) { init_wtab(); } - ~fft() { } - - unsigned int size() { return n_; } - - // Transform the contents of the array "a", leaving results in - // bit-reversed order. - void - br_transform(I a) { - unsigned int i, j, m, n; - typename twiddle_vector::iterator wp; // twiddle factor pointer - I p, q; - - // n is the number of points in each subtransform (butterfly group) - // m is the number of subtransforms (butterfly groups), = n_ / n - // i is the index of the first point in the current butterfly group - // j is the number of the butterfly within the group - - for (n = 2, m = n_ / 2; n <= n_; n *= 2 , m /= 2) // each stage - for (i = 0; i < n_; i += n) // each butterfly group - for (j = 0, wp = wtab.begin(), p = a + i, q = a + i + n / 2; - j < n / 2; - j++, wp += m, p++, q++) // each butterfly - { - C temp((*q) * (*wp)); - *q = *p - temp; - *p += temp; - } - } - - void - bit_reverse(I a) { - unsigned int i, j; - for (i = 0, j = 0; i < n_; i++, j = bitrev_inc(j)) { - if (i < j) - std::swap(*(a + i), *(a + j)); - } - } - - void - reverse(I a) { - for (unsigned int i = 1; i < n_ / 2; i++) - std::swap(*(a + i), *(a + n_ - i)); - } - - // in-place - void - transform(I a) { - bit_reverse(a); - br_transform(a); - } - - void - itransform(I a) { - reverse(a); - transform(a); - } - - // out-of-place - // XXX const - void - transform(I in, I out) { - std::copy(in, in + n_, out); - transform(out); - } - - void - itransform(I in, I out) { - std::copy(in, in + n_, out); - itransform(out); - } - -private: - // Initialize twiddle factor array - void init_wtab() { - size_t wt_size = wtab.size(); - for (size_t i = 0; i < wt_size; ++i) { - double arg = (-2.0 * M_PI / n_) * i; - wtab[i] = C(cos(arg), sin(arg)); - } - } - - unsigned int - bitrev_inc(unsigned int i) { - unsigned int carry = n_; - do { - carry >>= 1; - unsigned int new_i = i ^ carry; - carry &= i; - i = new_i; - } while(carry); - return i; - } - - // Size of the transform - unsigned int n_; - - // Twiddle factor array (size n / 2) - twiddle_vector wtab; -}; - -// Real FFT -// -// This is a trivial implementation offering no performance advantage -// over a complex FFT. It is intended as a placeholder to be -// overridden with a specialization, and as a reference implementation -// to compare the results of specializations against. -// - -template -struct rfft { - typedef typename std::iterator_traits::value_type C; // complex - typedef typename C::value_type T; // float/double - typedef T *RI; // Real iterator - typedef const T *CONST_RI; - - rfft(unsigned int n): cf(n) { } - ~rfft() { } - - void - transform(CONST_RI in, CI out) { - size_t n = cf.size(); - C *tmp = new C[n]; - C *out_tmp = new C[n]; - std::copy(in, in + cf.size(), tmp); // Real to complex - cf.transform(tmp, out_tmp); - delete [] tmp; -#if GABORATOR_REAL_FFT_NEGATIVE_FQS - std::copy(out_tmp, out_tmp + n, out); -#else - std::copy(out_tmp, out_tmp + n / 2 + 1, out); -#endif - delete [] out_tmp; - } - - void - itransform(CI in, RI out) { - size_t n = cf.size(); - // Make sure not to use the negative frequency part of "in", - // because it may not be valid. - C *in_tmp = new C[n]; - for (size_t i = 0; i < n / 2 + 1; i++) { - in_tmp[i] = in[i]; - } - for (size_t i = 1; i < n / 2; i++) { - in_tmp[n - i] = conj(in[i]); - } - C *tmp = new C[n]; - cf.itransform(in_tmp, tmp); - for (size_t i = 0; i < n; i++) { - *out++ = tmp[i].real(); - } - delete [] tmp; - delete [] in_tmp; - } - - fft cf; -}; - -} // Namespace - -#endif diff --git a/lib/gaborator/gaborator/fft_pffft.h b/lib/gaborator/gaborator/fft_pffft.h deleted file mode 100644 index b3e64b9..0000000 --- a/lib/gaborator/gaborator/fft_pffft.h +++ /dev/null @@ -1,214 +0,0 @@ -// -// Fast Fourier transform using PFFFT -// -// Copyright (C) 2017-2021 Andreas Gustafsson. This file is part of -// the Gaborator library source distribution. See the file LICENSE at -// the top level of the distribution for license information. -// - -#ifndef _GABORATOR_FFT_PFFFT_H -#define _GABORATOR_FFT_PFFFT_H - -#include -#include - -#include -#include - - -#include "pffft.h" - -// XXX disable in production -#ifdef __x86_64__ -#define GABORATOR_PFFFT_CHECK_ALIGN(p) assert((((uint64_t)(p)) & 0xF) == 0) -#else -#define GABORATOR_PFFFT_CHECK_ALIGN(p) do {} while (0) -#endif - -namespace gaborator { - -template <> -struct fft *> { - typedef std::complex *I; - typedef const std::complex *CONST_I; - typedef std::iterator_traits::value_type C; // complex - typedef C::value_type T; // float/double - - fft(unsigned int n_): n(n_) { - setup = pffft_new_setup(n, PFFFT_COMPLEX); - assert(setup); - } - ~fft() { - pffft_destroy_setup(setup); - } - - unsigned int size() { return n; } - - // in-place - void - transform(I a) { - pffft_transform_ordered(setup, (float *)a, (float *)a, NULL, PFFFT_FORWARD); - } - - void - itransform(I a) { - pffft_transform_ordered(setup, (float *)a, (float *)a, NULL, PFFFT_BACKWARD); - } - - // out-of-place - void - transform(CONST_I in, I out) { - GABORATOR_PFFFT_CHECK_ALIGN(in); - GABORATOR_PFFFT_CHECK_ALIGN(out); - pffft_transform_ordered(setup, (const float *)in, (float *)out, NULL, PFFFT_FORWARD); - } - - void - itransform(CONST_I in, I out) { - GABORATOR_PFFFT_CHECK_ALIGN(in); - GABORATOR_PFFFT_CHECK_ALIGN(out); - pffft_transform_ordered(setup, (const float *)in, (float *)out, NULL, PFFFT_BACKWARD); - } - -private: - // Size of the transform - unsigned int n; - PFFFT_Setup *setup; -}; - -// Support transforming std::vector>::iterator - -template <> -struct fft>::iterator>: - public fft *> -{ - typedef fft *> base; - typedef std::vector>::iterator I; - fft(unsigned int n_): fft *>(n_) { } - void - transform(I a) { - base::transform(&(*a)); - } - void - itransform(I a) { - base::itransform(&(*a)); - } - void - transform(I in, I out) { - base::transform(&(*in), &(*out)); - } - void - itransform(I in, I out) { - base::itransform(&(*in), &(*out)); - } -}; - -// Use fftpack for double precision - -#define FFTPACK_DOUBLE_PRECISION 1 -#include "fftpack.h" -#undef FFTPACK_DOUBLE_PRECISION - -template <> -struct fft *> { - typedef std::complex *I; - typedef const std::complex *CONST_I; - typedef std::iterator_traits::value_type C; // complex - typedef C::value_type T; // float/double - - fft(unsigned int n_): n(n_), wsave(4 * n_ + 15) { - cffti(n, wsave.data()); - } - ~fft() { - } - - unsigned int size() { return n; } - - // in-place - void - transform(I a) { - cfftf(n, (double *)a, wsave.data()); - } - - void - itransform(I a) { - cfftb(n, (double *)a, wsave.data()); - } - - // out-of-place - void - transform(CONST_I in, I out) { - std::copy(in, in + n, out); - transform(out); - } - - void - itransform(CONST_I in, I out) { - std::copy(in, in + n, out); - itransform(out); - } - -private: - // Size of the transform - unsigned int n; - std::vector wsave; -}; - -// Real FFT - -template <> -struct rfft *> { - typedef std::complex *CI; // Complex iterator - typedef const std::complex *CONST_CI; - typedef typename std::iterator_traits::value_type C; // complex - typedef typename C::value_type T; // float/double - typedef T *RI; // Real iterator - typedef const T *CONST_RI; - - rfft(unsigned int n_): n(n_) { - setup = pffft_new_setup(n, PFFFT_REAL); - assert(setup); - } - ~rfft() { - pffft_destroy_setup(setup); - } - - unsigned int size() { return n; } - - // out-of-place only - void - transform(CONST_RI in, CI out) { - GABORATOR_PFFFT_CHECK_ALIGN(in); - GABORATOR_PFFFT_CHECK_ALIGN(out); - pffft_transform_ordered(setup, in, (float *) out, NULL, PFFFT_FORWARD); - C tmp = out[0]; -#if GABORATOR_REAL_FFT_NEGATIVE_FQS - for (unsigned int i = 1; i < (n >> 1); i++) - out[n - i] = conj(out[i]); -#endif - out[0] = C(tmp.real(), 0); - out[n >> 1] = C(tmp.imag(), 0); - } - - // Note: this temporarily modifies in[0], in spite of the const - void - itransform(CONST_CI in, RI out) { - GABORATOR_PFFFT_CHECK_ALIGN(in); - GABORATOR_PFFFT_CHECK_ALIGN(out); - C tmp = in[0]; - const_cast(in)[0] = C(tmp.real(), in[n >> 1].real()); - pffft_transform_ordered(setup, (const float *) in, out, NULL, PFFFT_BACKWARD); - const_cast(in)[0] = tmp; - } - -private: - // Size of the transform - unsigned int n; - PFFFT_Setup *setup; -}; - -#undef GABORATOR_PFFFT_CHECK_ALIGN - -} // namespace - -#endif diff --git a/lib/gaborator/gaborator/fft_vdsp.h b/lib/gaborator/gaborator/fft_vdsp.h deleted file mode 100644 index f0fe232..0000000 --- a/lib/gaborator/gaborator/fft_vdsp.h +++ /dev/null @@ -1,209 +0,0 @@ -// -// Fast Fourier transform using the Apple vDSP framework -// -// Copyright (C) 2013-2021 Andreas Gustafsson. This file is part of -// the Gaborator library source distribution. See the file LICENSE at -// the top level of the distribution for license information. -// - -#ifndef _GABORATOR_FFT_VDSP_H -#define _GABORATOR_FFT_VDSP_H - -#include - -#include -#include - -#include - -#include -#include -#include -#include - -#include - -namespace gaborator { - -static inline int log2_int_exact(int n) { - // n must be a power of two - assert(n != 0 && ((n & (n >> 1)) == 0)); - int r = 0; - for (;;) { - n >>= 1; - if (n == 0) - break; - r++; - } - return r; -} - -template <> -struct fft *> { - typedef std::complex *I; - typedef typename std::iterator_traits::value_type C; // complex - typedef typename C::value_type T; // float/double - - fft(unsigned int n_): n(n_), log2n(log2_int_exact(n)) { - setup = vDSP_create_fftsetup(log2n, kFFTRadix2); - } - ~fft() { - vDSP_destroy_fftsetup(setup); - } - - unsigned int size() { return n; } - - // in-place - void - transform(I a) { - DSPSplitComplex s; - // XXX this result in disoptimal alignment - s.realp = (float *) a; - s.imagp = (float *) a + 1; - vDSP_fft_zip(setup, &s, 2, log2n, kFFTDirection_Forward); - } - - void - itransform(I a) { - DSPSplitComplex s; - s.realp = (float *) a; - s.imagp = (float *) a + 1; - vDSP_fft_zip(setup, &s, 2, log2n, kFFTDirection_Inverse); - } - - // out-of-place - // XXX const - void - transform(I in, I out) { - DSPSplitComplex si; - si.realp = (float *) in; - si.imagp = (float *) in + 1; - DSPSplitComplex so; - so.realp = (float *) out; - so.imagp = (float *) out + 1; - vDSP_fft_zop(setup, - &si, 2, - &so, 2, - log2n, kFFTDirection_Forward); - } - - void - itransform(I in, I out) { - DSPSplitComplex si; - si.realp = (float *) in; - si.imagp = (float *) in + 1; - DSPSplitComplex so; - so.realp = (float *) out; - so.imagp = (float *) out + 1; - vDSP_fft_zop(setup, - &si, 2, - &so, 2, - log2n, kFFTDirection_Inverse); - } - -private: - // Size of the transform - unsigned int n; - unsigned int log2n; - FFTSetup setup; -}; - -// Real FFT - -template <> -struct rfft *> { - typedef std::complex *CI; // Complex iterator - typedef const std::complex *CONST_CI; - typedef typename std::iterator_traits::value_type C; // complex - typedef typename C::value_type T; // float/double - typedef T *RI; // Real iterator - typedef const T *CONST_RI; - - rfft(unsigned int n_): n(n_), log2n(log2_int_exact(n)) { - setup = vDSP_create_fftsetup(log2n, kFFTRadix2); - } - ~rfft() { - vDSP_destroy_fftsetup(setup); - } - - unsigned int size() { return n; } - - // out-of-place only - void - transform(CONST_RI in, CI out) { - DSPSplitComplex si; - si.realp = (float *) in; - si.imagp = (float *) in + 1; - DSPSplitComplex so; - so.realp = (float *) out; - so.imagp = (float *) out + 1; - vDSP_fft_zrop(setup, - &si, 2, - &so, 2, - log2n, kFFTDirection_Forward); - // Undo vDSP scaling - for (unsigned int i = 0; i < (n >> 1); i++) - out[i] *= (T)0.5; - C tmp = out[0]; -#if GABORATOR_REAL_FFT_NEGATIVE_FQS - for (unsigned int i = 1; i < (n >> 1); i++) - out[n - i] = conj(out[i]); -#endif - out[0] = C(tmp.real(), 0); - out[n >> 1] = C(tmp.imag(), 0); - } - - void - itransform(CONST_CI in, RI out) { - C tmp = in[0]; - const_cast(in)[0] = C(tmp.real(), in[n >> 1].real()); - DSPSplitComplex si; - si.realp = (float *) in; - si.imagp = (float *) in + 1; - DSPSplitComplex so; - so.realp = (float *) out; - so.imagp = (float *) out + 1; - vDSP_fft_zrop(setup, - &si, 2, - &so, 2, - log2n, kFFTDirection_Inverse); - const_cast(in)[0] = tmp; - } - -private: - // Size of the transform - unsigned int n; - unsigned int log2n; - FFTSetup setup; -}; - -// Support transforming std::vector>::iterator - -template <> -struct fft>::iterator>: - public fft *> -{ - typedef fft *> base; - typedef typename std::vector>::iterator I; - fft(unsigned int n_): fft *>(n_) { } - void - transform(I a) { - base::transform(&(*a)); - } - void - itransform(I a) { - base::itransform(&(*a)); - } - void - transform(I in, I out) { - base::transform(&(*in), &(*out)); - } - void - itransform(I in, I out) { - base::itransform(&(*in), &(*out)); - } -}; - -} // Namespace - -#endif diff --git a/lib/gaborator/gaborator/gaborator.h b/lib/gaborator/gaborator/gaborator.h deleted file mode 100644 index 358f2e2..0000000 --- a/lib/gaborator/gaborator/gaborator.h +++ /dev/null @@ -1,2966 +0,0 @@ -// -// Constant Q spectrum analysis and resynthesis -// -// Copyright (C) 2015-2021 Andreas Gustafsson. This file is part of -// the Gaborator library source distribution. See the file LICENSE at -// the top level of the distribution for license information. -// - -#ifndef _GABORATOR_GABORATOR_H -#define _GABORATOR_GABORATOR_H - -#define __STDC_LIMIT_MACROS - -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include - -#include "gaborator/fft.h" -#include "gaborator/gaussian.h" -#include "gaborator/affine_transform.h" -#include "gaborator/pod_vector.h" -#include "gaborator/pool.h" -#include "gaborator/ref.h" -#include "gaborator/vector_math.h" - - -namespace gaborator { - -using std::complex; - -// An integer identifying an audio sample -typedef int64_t sample_index_t; - -// An integer identifying a coefficient -typedef int64_t coef_index_t; - -// An integer identifying a slice -typedef int64_t slice_index_t; - -// See https://tauday.com/tau-manifesto -static const double tau = 2.0 * M_PI; - -// Round up to next higher or equal power of 2 - -inline int -next_power_of_two(int x) { - --x; - x |= x >> 1; - x |= x >> 2; - x |= x >> 4; - x |= x >> 8; - x |= x >> 16; - return x + 1; -} - -// Determine if x is a power of two. -// Note that this considers 0 to be a power of two. - -static inline bool -is_power_of_two(unsigned int x) { - return (x & (x - 1)) == 0; -} - -// Given a power of two v, determine log2(v) -// https://graphics.stanford.edu/~seander/bithacks.html#DetermineIfPowerOf2 - -static inline unsigned int whichp2(unsigned int v) { - assert(is_power_of_two(v)); - unsigned int r = (v & 0xAAAAAAAA) != 0; - r |= ((v & 0xCCCCCCCC) != 0) << 1; - r |= ((v & 0xF0F0F0F0) != 0) << 2; - r |= ((v & 0xFF00FF00) != 0) << 3; - r |= ((v & 0xFFFF0000) != 0) << 4; - return r; -} - -// Floor division: return the integer part of a / b -// rounded down (not towards zero). For positive b only. - -inline int64_t floor_div(int64_t a, int64_t b) { - assert(b > 0); - if (a >= 0) - return a / b; - else - return (a - b + 1) / b; -} - -// Floating point modulus, the remainder r of a / b -// satisfying 0 <= r < b even for negative a. -// For positive b only. - -static inline double -sane_fmod(double a, double b) { - assert(b > 0); - double m = fmod(a, b); - if (m < 0) - m += b; - return m; -} - -// Do an arithmetic left shift of a 64-bit signed integer. This is -// what a << b ought to do, but according to the C++11 draft (n3337), -// section 5.8, that invokes undefined behavior when a is negative. -// GCC is actually smart enough to optimize this into a single shlq -// instruction. -// -// No corresponding kludge is needed for right shifts, because a right -// shift of a negative signed integer is implementation-defined, not -// undefined, and we trust implementations to define it sanely. - -static inline int64_t -shift_left(int64_t a, unsigned int b) { - if (a < 0) - return -(((uint64_t) -a) << b); - else - return (((uint64_t) a) << b); -} - -// Convert between complex types - -template -complex c2c(complex c) { return complex(c.real(), c.imag()); } - -// Convert a sequence of complex values to real - -template -O complex2real(I b, I e, O o) { - while (b != e) { - *o++ = (*b++).real(); - } - return o; -} - -// A vector-like object that allows arbitrary integer indices -// (positive or negative, but excluding the largest possible integer) -// and automatically resizes the storage. Uses storage proportional -// to the difference between the smallest and largest index value (for -// example, if indices range from -102 to -100 (inclusive), memory use -// is on the order of 3 elements). -// -// T is the element type -// I is the integer index type - -template -struct range_vector { - range_vector() { - init_bounds(); - } - range_vector(const range_vector &) = default; - range_vector &operator=(const range_vector &rhs) = default; - range_vector(range_vector &&rhs): - v(std::move(rhs.v)), - lower(rhs.lower), - upper(rhs.upper) - { - rhs.init_bounds(); - } - range_vector &operator=(range_vector &&rhs) { - if (this == &rhs) - return *this; - v = std::move(rhs.v); - lower = rhs.lower; - upper = rhs.upper; - rhs.init_bounds(); - return *this; - } -private: - void init_bounds() { - lower = std::numeric_limits::max(); - upper = std::numeric_limits::min(); - } - T *unchecked_get(I i) { - return &v[(size_t)(i & ((I)v.size() - 1))]; - } - const T *unchecked_get(I i) const { - return &v[i & ((I)v.size() - 1)]; - } - -public: - // Get a pointer to an existing element, or null if out of range - const T * - get(I i) const { - if (! has_index(i)) - return 0; - return unchecked_get(i); - } - - // Note: Reference returned becomes invalid when range_vector - // is changed - T & - get_or_create(I i) { - if (! has_index(i)) - extend(i); - return *unchecked_get(i); - - } - - // Get a reference to the element at index i, which must be valid - T & - get_existing(I i) { - assert(has_index(i)); - return *unchecked_get(i); - } - - // Const version of the above - const T & - get_existing(I i) const { - assert(has_index(i)); - return *unchecked_get(i); - } - -private: - void extend(I i) { - I new_lower = lower; - I new_upper = upper; - if (i < lower) - new_lower = i; - if (i + 1 > upper) - new_upper = i + 1; - I old_size = v.size(); - I new_need = new_upper - new_lower; - if (new_need > old_size) { - if (old_size == 0) { - v.resize(1); - } else { - I new_size = old_size; - while (new_size < new_need) - new_size *= 2; - v.resize(new_size); - if (old_size) { - for (I j = lower; j < upper; j++) { - I jo = j & (old_size - 1); - I jn = j & (new_size - 1); - if (jo != jn) - std::swap(v[jo], v[jn]); - } - } - } - } - lower = new_lower; - upper = new_upper; - } - -public: - // Erase the elements whose index is less than "limit" - void erase_before(I limit) { - I i = lower; - for (; i < upper && i < limit; i++) - *unchecked_get(i) = T(); - lower = i; - } - - void clear() { - v.clear(); - init_bounds(); - } - - I begin_index() const { return lower; } - I end_index() const { return upper; } - bool empty() const { return lower >= upper; } - bool has_index(I i) const { return i >= lower && i < upper; } - -private: - std::vector v; - I lower, upper; -}; - -// Calculate the size of the alias-free part (the "filet") -// of a signal slice of size "fftsize" - -static inline unsigned int filet_part(unsigned int fftsize) { - return fftsize >> 1; -} - -// Calculate the size of the padding (the "fat") at each -// end of a signal slice of size "fftsize" - -static inline unsigned fat_part(unsigned int fftsize) { - return fftsize >> 2; -} - -// Per-band, per-plan data - -template -struct band_plan { - typedef complex C; - unsigned int sftsize; // Size of "short FFT" spanning the band - unsigned int sftsize_log2; // log2(sftsize) - fft *sft; // Fourier transform for windows, of size sftsize - std::vector kernel; // Frequency-domain filter kernel - std::vector dual_kernel; // Dual of the above - pod_vector shift_kernel; // Complex exponential for fractional frequency compensation - pod_vector shift_kernel_conj; // Conjugate of the above - int fq_offset_int; // Frequency offset in bins (big-FFT bin of left window edge) - double center; // Center frequency in units of FFT bins - int icenter; // Center frequency rounded to nearest integer FFT bin -}; - -// Frequency band parameters shared between octaves - -template -struct band_params: public refcounted { - typedef complex C; - bool dc; // True iff this is the lowpass (DC) band - double ff; // Center (bp) or corner (lp) frequency in units of the sampling frequency - double ffsd; // Standard deviation of the bandpass Gaussian, as fractional frequency - unsigned int step; // Signal samples per coefficient sample - unsigned int step_log2; // log2(step) - double ff_support; // Filter support in frequency domain - double time_support; // Filter support in time domain, in octave subsamples - std::vector> anl_plans; - std::vector> syn_plans; -}; - -// Downsampling parameters. These have some similarity to band -// parameters, but only some. For example, these may use a real -// rather than complex FFT for the "short FFT". - -template -struct downsampling_params { - typedef complex C; - unsigned int sftsize; - std::vector kernel; // Frequency-domain filter kernel - std::vector dual_kernel; -#if GABORATOR_USE_REAL_FFT - rfft *rsft; -#else - fft *sft; -#endif -}; - -// Forward declarations -template struct analyzer; -template struct zone; -template> struct coefs; -template struct sliced_coefs; -template *, class C = complex> - struct row_source; -template *, class C = complex> - struct row_dest; -template *, class C = complex> - struct row_add_dest; - -// Abstract class for tracking changes to a coefficient set. -// This may be used for updating a resolution pyramid of -// magnitude data. - -template -struct shadow { - virtual ~shadow() { } - virtual void update(const coefs> &msc, - sample_index_t i0, sample_index_t i1) = 0; - // Deprecated - virtual void update(int oct, int ze, const sliced_coefs> &sc, - slice_index_t sli0, slice_index_t sli1) = 0; -}; - -// Per-band coefficient metadata, shared between octaves. - -struct band_coefs_meta { - unsigned int slice_len; // Number of coefficients per slice - unsigned short slice_len_log2; // Log2 of the above - // Log2 of the downsampling factor of the coefficients in - // this band relative to the signal samples. The value - // applies as such in the top octave; in other octaves, - // the octave number needs to be added. - unsigned short step_log2; - // Offset of the beginning of this band in the data array - unsigned int band_offset; -}; - -// Per-zone coefficient metadata. This contains information -// describing an "oct_coefs" structure that is common to many -// instances and should not be duplicated in each one. - -struct zone_coefs_meta { - typedef std::vector band_vector; - void init(const band_vector &bands_) { - bands = bands_; - unsigned int offset = 0; - for (band_coefs_meta &b: bands) { - b.band_offset = offset; - offset += b.slice_len; - } - total_size = offset; - } - band_vector bands; - // Total size of data array (in elements, not bytes) - unsigned int total_size; -}; - -// Per-octave coefficient metadata. -// Cf. struct octave - -struct oct_coefs_meta { - zone_coefs_meta *z; - unsigned int n_bands_above; // Total number of bands in higher octaves -}; - -// Coefficient metadata for multirate coefficients. - -struct coefs_meta: public refcounted { - coefs_meta() = default; - coefs_meta(const coefs_meta &) = delete; - unsigned int n_octaves; - unsigned int n_bands_total; - unsigned int bands_per_octave; - unsigned int slice_len; // octave subsamples per slice - std::vector zones; - std::vector octaves; -}; - -// Split a "global band number" gbno into an octave and band -// number within octave ("obno"). -// -// Global band numbers start at 0 for the band at or close to -// fs/2, and increase towards lower frequencies. -// -// Include the DC band if "dc" is true. -// Returns true iff gbno is valid. - -static inline bool -bno_split(const coefs_meta &meta, int gbno, int &oct, unsigned int &obno, bool dc) { - if (gbno < 0) { - // Above top octave - return false; - } else if (gbno >= (int)meta.n_bands_total - 1) { - // At or below DC - if (gbno == (int)meta.n_bands_total - 1) { - // At DC - if (dc) { - oct = meta.n_octaves - 1; - obno = 0; - return true; - } else { - return false; - } - } else { - // Below DC - return false; - } - } else { - // Within bandpass region - // Start by determining the octave - int n_bands_top_octave = (int)meta.octaves[0].z->bands.size(); - if (gbno < n_bands_top_octave) { - // Top octave - oct = 0; - obno = n_bands_top_octave - 1 - gbno; - return true; - } - gbno -= n_bands_top_octave; - int oct_tmp = 1 + gbno / meta.bands_per_octave; - int obno_tmp = gbno % meta.bands_per_octave; - oct = oct_tmp; - // Now determine the band within the octave. - // obno_tmp counts down, but obno counts up. - obno = (unsigned int)meta.octaves[oct_tmp].z->bands.size() - 1 - obno_tmp; - return true; - } -} - -// The inverse of bno_split(). Returns a gbno. The arguments must -// be valid. - -static inline -int bno_merge(const coefs_meta &meta, int oct, unsigned int obno) { - unsigned int n_bands = (unsigned int)meta.octaves[oct].z->bands.size(); - assert(obno < n_bands); - int bno_from_end = n_bands - 1 - obno; - return bno_from_end + meta.octaves[oct].n_bands_above; -} - -// Coefficients of a single octave for a single input signal slice. -// C is the coefficient type, typically complex but can also -// be e.g. unsigned int to store cluster numbers, or float to store -// magnitudes. - -template -struct oct_coefs: public refcounted { - oct_coefs(const zone_coefs_meta &zmeta_, bool clear_ = true): - zmeta(zmeta_), - data(zmeta.total_size), - bands(*this) - { - if (clear_) - clear(); - } - oct_coefs(const oct_coefs &) = delete; - uint64_t estimate_memory_usage() const { - return zmeta.total_size * sizeof(C) + sizeof(*this); - } - void clear() { - std::fill(data.begin(), data.end(), C()); - } - // Deep copy - oct_coefs &operator=(const oct_coefs &rhs) { - assert(data.size() == rhs.data.size()); - memcpy(data.data(), rhs.data.data(), data.size() * sizeof(C)); - return *this; - } - - const zone_coefs_meta &zmeta; - - // The data for all the bands are allocated together - // as a single vector to reduce the number of allocations - pod_vector data; - // Vector-like collection of pointers into "data", one for each band - struct band_array { - band_array(oct_coefs &outer_): outer(outer_) { } - C *operator[](size_t i) const { - return outer.data.data() + outer.zmeta.bands[i].band_offset; - } - size_t size() const { return outer.zmeta.bands.size(); } - oct_coefs &outer; - } bands; -}; - -// Add the oct_coefs "b" to the oct_coefs "a" - -template -void add(oct_coefs &a, const oct_coefs &b) { - size_t n_bands = a.bands.size(); - assert(n_bands == b.bands.size()); - for (size_t obno = 0; obno < n_bands; obno++) { - unsigned int len = a.zmeta.bands[obno].slice_len; - C *band_a = a.bands[obno]; - C *band_b = b.bands[obno]; - for (unsigned int j = 0; j < len; j++) { - band_a[j] += band_b[j]; - } - } -} - -// Sliced coefficients. These cover an arbitrary time range, but only -// a single octave. Template argument is as for struct oct_coefs. -// This is default constructible so that we can create an array of -// them, but not usable until "meta" has been set up. - -template -struct sliced_coefs { - typedef range_vector>, slice_index_t> slices_t; - uint64_t estimate_memory_usage() const { - unsigned int n = 0; - size_t size_each = 0; - for (slice_index_t sl = slices.begin_index(); sl < slices.end_index(); sl++) { - const ref> &t = slices.get_existing(sl); - if (t) { - if (! size_each) - size_each = (size_t)t->estimate_memory_usage(); - n++; - } - } - return n * size_each; - } - void clear() { - slices.clear(); - } - zone_coefs_meta *meta; - slices_t slices; -}; - -// Get a pointer to an existing existing coefficient slice, -// or null if one does not exist. Like get_or_create_coefs(), -// this hides the distinction between the two types of nonexistence. - -template -oct_coefs *get_existing_coefs(const sliced_coefs &sc, - slice_index_t i) -{ - const ref> *p = sc.slices.get(i); - if (! p) - return 0; - return p->get(); -} - -// Get an existing coefficient slice, or create a new one. Note that -// this hides the distinction between two types of nonexistence: that -// of slices outside the range of the range_vector, and that of -// missing slices within the range (having a null ref). CT is the -// coefficient type, which is typically C aka complex, but can be -// different, for example float to represent magnitudes. - -template -oct_coefs &get_or_create_coefs(sliced_coefs &sc, slice_index_t i) { - ref> &p(sc.slices.get_or_create(i)); - if (! p) - p.reset(new oct_coefs(*sc.meta)); - return *p; -} - -// Return the signal sample time corresponding to coefficient sample 0 -// of coefficient slice 0, for slices of length len. It would be nice -// if this were zero, but for historical reasons, it's offset by half -// a slice (corresponding to the analysis fat). - -static inline int coef_offset(int len) { - return len >> 1; -} - -// Get the base 2 logarithm of the downsampling factor of -// band "obno" in octave "oct" - -static inline int -band_scale_exp(const zone_coefs_meta &meta, int oct, unsigned int obno) { - return meta.bands[obno].step_log2 + oct; -} - -// Return the coefficient index (the time in terms of coefficient -// subsamples) of the first cofficient of slice "sli" of band -// "obno" in octave "oct" - -static inline coef_index_t -coef_time(const zone_coefs_meta &meta, slice_index_t sli, int oct, int obno) { - int len = meta.bands[obno].slice_len; - return coef_offset(len) + sli * len; -} - -// Return the sample index (the time in terms of samples) time of -// coefficient "i" in slice "sli" of band "obno" in octave "oct" - -static inline sample_index_t -sample_time(const zone_coefs_meta &meta, slice_index_t sli, int i, int oct, int obno) { - coef_index_t sst = coef_time(meta, sli, oct, obno) + i; - return shift_left(sst, band_scale_exp(meta, oct, obno)); -} - -// Multirate sliced coefficients. These cover an arbitrary time -// range and the full frequency range (all octaves). -// Template arguments: -// T analyzer sample data type -// C coefficient data type -// Note default for template argument C defined in forward declaration. - -template -struct coefs { - typedef C value_type; - coefs(const analyzer &anl_, shadow *shadow_ = 0): - octaves(anl_.n_octaves), shadow0(shadow_) - { - meta = anl_.cmeta_any.get(); - // Set up shortcut pointer to zone metadata in each octave - for (unsigned int oct = 0; oct < octaves.size(); oct++) - octaves[oct].meta = meta->octaves[oct].z; - } - uint64_t estimate_memory_usage() const { - uint64_t s = 0; - for (unsigned int oct = 0; oct < octaves.size(); oct++) - s += octaves[oct].estimate_memory_usage(); - return s; - } - void clear() { - for (unsigned int oct = 0; oct < octaves.size(); oct++) - octaves[oct].clear(); - } - coefs_meta *meta; - std::vector> octaves; - shadow *shadow0; -}; - -// Read coefficients i0..i1 of band gbno in msc into buf. - -template -void read(const coefs &msc, int gbno, - coef_index_t i0, coef_index_t i1, C *buf) -{ - int oct; - unsigned int obno; // Band number within octave - bool valid = gaborator::bno_split(*msc.meta, gbno, oct, obno, true); - assert(valid); - row_source(msc, oct, obno)(i0, i1, buf); -} - -template -void write(coefs &msc, int gbno, - coef_index_t i0, coef_index_t i1, C *buf) -{ - int oct; - unsigned int obno; // Band number within octave - bool valid = gaborator::bno_split(*msc.meta, gbno, oct, obno, true); - assert(valid); - row_dest(msc, oct, obno)(i0, i1, buf); -} - -template -void add(coefs &msc, int gbno, - coef_index_t i0, coef_index_t i1, C *buf) -{ - int oct; - unsigned int obno; // Band number within octave - bool valid = gaborator::bno_split(*msc.meta, gbno, oct, obno, true); - assert(valid); - row_add_dest(msc, oct, obno)(i0, i1, buf); -} - -// Return the base 2 logarithm of the time step (aka downsampling -// factor) of band "gbno". - -static inline -unsigned int band_step_log2(const coefs_meta &meta, int gbno) { - int oct; - unsigned int obno; - bool valid = bno_split(meta, gbno, oct, obno, true); - assert(valid); - return band_scale_exp(*meta.octaves[oct].z, oct, obno); -} - - -// Convert a signal time t into a coefficient sample -// index. t must coincide with a coefficient sample time. - -static inline -coef_index_t t2i_exact(const coefs_meta &meta, int gbno, sample_index_t t) { - int shift = band_step_log2(meta, gbno); - int64_t mask = ((sample_index_t)1 << shift) - 1; - assert((t & mask) == 0); - return t >> shift; -} - -// Read a single coefficient sample at signal time t, -// which must coincide with a coefficient sample time - -template -C read1t(const coefs &msc, int gbno, sample_index_t t) { - coef_index_t i = t2i_exact(*msc.meta, gbno, t); - C c; - read(msc, gbno, i, i + 1, &c); - return c; -} - -// Read a single coefficient sample at signal time t, -// which must coincide with a coefficient sample time - -template -void write1t(coefs &msc, int gbno, sample_index_t t, C c) { - coef_index_t i = t2i_exact(*msc.meta, gbno, t); - write(msc, gbno, i, i + 1, &c); -} - -// Perform an fftshift of the range between iterators a and b. -// Not optimized - not for use in inner loops. - -template -void fftshift(I b, I e) { - size_t len = e - b; - assert(len % 2 == 0); - for (size_t i = 0; i < len / 2; i++) - std::swap(*(b + i), *(b + len / 2 + i)); -} - -// Given an unsigned index i into an FFT of a power-of-two size -// "size", return the corresponding signed index, ranging from -size/2 -// to size/2-1. This is equivalent to sign extension of an integer of -// log2(size) bits; see Hacker's Delight, page 18. This can be used -// to convert an FFT index into a frequency (positive or negative; -// note that Nyquist is considered negatitive = -fs/2). - -static inline int signed_index(unsigned int i, unsigned int size) { - unsigned int t = size >> 1; - return (i ^ t) - t; -} - -// Evaluate a Gaussian windowed lowpass filter frequency response. -// This is the convolution of a rectangle centered at f=0 and a Gaussian, -// and corresponds to a Gaussian windowed sinc in the time domain. -// The -6 dB cutoff freqency is ff_cutoff (a fractional frequency), -// the standard deviation of the Gaussian is ff_sd, and the frequency -// response is evaluated at ff. The frequency response is smooth at -// f=0 even if the transition bands overlap. - -inline double -gaussian_windowed_lowpass_1(double ff_cutoff, double ff_sd, double ff) { - return - // A rectangle is the sum of a rising step and a later falling - // step, or the difference between a rising step and a later - // rising step. By linearity, a Gaussian filtered rectangle - // is the difference between two Gaussian filtered rising - // steps. - gaussian_edge(ff_sd, -ff + ff_cutoff) - - gaussian_edge(ff_sd, -ff - ff_cutoff); -} - -// Fill a sequence with a frequency-ddomain lowpass filter as above. -// The returned filter covers the full frequency range from 0 to fs -// (with negative frequencies at the end, the usual convention for FFT -// spectra). -// -// When center=true, construct a time-domain window instead, -// passing the center of the time-domain signal. -// -// The result is stored between iterators b and e, which must have a -// real value_type. - -template -inline void gaussian_windowed_lowpass(double ff_cutoff, double ff_sd, - I b, I e, bool center = false) -{ - size_t len = e - b; - double inv_len = 1.0 / len; - for (I it = b; it != e; ++it) { - size_t i = it - b; - double thisff; - if (center) - // Symmetric around center - thisff = std::abs(i - (len * 0.5)) * inv_len; - else - // Symmetric around zero - thisff = (i > len / 2 ? len - i : i) * inv_len; - *it = gaussian_windowed_lowpass_1(ff_cutoff, ff_sd, thisff); - } -} - -// A set of octaves having identical parameters form a "zone", -// and their shared parameters are stored in a "struct zone". - -template -struct zone: public refcounted { - zone(): n_bands(0) { } - ~zone() { } - // Zone number, 0..3 - unsigned int zno; - // Total number of bands, including DC band if lowest octave - unsigned int n_bands; - unsigned int max_step_log2; - // Band parameters by increasing frequency; DC band is index 0 if - // present - std::vector>> bandparams; - std::vector>> mock_bandparams; -}; - -template -struct octave { - zone *z; -}; - -// Helper function for pushing parameters onto the vectors in struct zone - -template -void push(std::vector>> &v, band_params *p) { - v.push_back(ref>(p)); -} - -// Phase conventions: coef_phase::absolute means the phase of a -// coefficient at time tc is relative to e^(i tau f t), and -// coef_phase::local means it is relative to -// e^(i tau f (t - tc)) - -enum class coef_phase { global, local }; - -// A set of spectrum analysis parameters - -struct parameters { - parameters(unsigned int bands_per_octave_, - double ff_min_, - double ff_ref_ = 1.0, - double overlap_ = 0.7, - double max_error_ = 1e-5): - bands_per_octave(bands_per_octave_), - ff_min(ff_min_), - ff_ref(ff_ref_), - overlap(overlap_), - max_error(max_error_), - coef_scale(1.0), - synthesis(true), - multirate(true) - { - init_v1(); - } - // Pseudo-constructor with version 1 defaults - static parameters v1(unsigned int bands_per_octave_, - double ff_min_, - double ff_ref_ = 1.0, - double overlap_ = 0.7, - double max_error_ = 1e-5) - { - parameters p(bands_per_octave_, ff_min_, ff_ref_, overlap_, max_error_); - p.init_v1(); - return p; - } - // Pseudo-constructor with version 2 defaults - static parameters v2(unsigned int bands_per_octave_, - double ff_min_, - double ff_ref_ = 1.0, - double overlap_ = 0.7, - double max_error_ = 1e-5) - { - parameters p(bands_per_octave_, ff_min_, ff_ref_, overlap_, max_error_); - p.init_v2(); - return p; - } - void init_v1() { - phase = coef_phase::global; - bandwidth_version = 1; - lowpass_version = 1; - } - void init_v2() { - phase = coef_phase::local; - bandwidth_version = 2; - lowpass_version = 2; - } - // Provide an operator< so that we can create a set or map of parameters - bool operator<(const parameters &b) const { -#define GABORATOR_COMPARE_LESS(member) do { \ - if (member < b.member) \ - return true; \ - if (member > b.member) \ - return false; \ - } while(0) - GABORATOR_COMPARE_LESS(bands_per_octave); - GABORATOR_COMPARE_LESS(ff_min); - GABORATOR_COMPARE_LESS(ff_ref); - GABORATOR_COMPARE_LESS(overlap); - GABORATOR_COMPARE_LESS(max_error); - GABORATOR_COMPARE_LESS(phase); - GABORATOR_COMPARE_LESS(bandwidth_version); - GABORATOR_COMPARE_LESS(lowpass_version); - GABORATOR_COMPARE_LESS(coef_scale); - GABORATOR_COMPARE_LESS(synthesis); - GABORATOR_COMPARE_LESS(multirate); -#undef GABORATOR_COMPARE_LESS - // Equal - return false; - } - bool operator==(const parameters &b) const { - return !((*this < b) || (b < *this)); - } - // The frequency increases by a factor of band_spacing from - // one bandpass band to the next. - double band_spacing_log2() const { - return 1.0 / bands_per_octave; - } - double band_spacing() const { - return exp2(band_spacing_log2()); - } - // The standard deviation of the Gaussian in units of the mean - double sd() const { - return overlap * - (bandwidth_version == 1 ? - band_spacing() - 1 : - log(2) / bands_per_octave); - } - - // Defining Q as the frequency divided by the half-power bandwidth, - // we get - // - // norm_gaussian(sd, hbw) = sqrt(2) - // - // (%i1) e1: exp(-(hbw * hbw) / (2 * sd * sd)) = 1 / sqrt(2); - // (%i2) solve(e1, hbw); - // (%o2) [hbw = - sqrt(log(2)) sd, hbw = sqrt(log(2)) sd] - double q() const { - return 1.0 / (2 * sqrt(log(2)) * sd()); - } - - template friend class analyzer; - unsigned int bands_per_octave; - double ff_min; - double ff_ref; - double overlap; - double max_error; - coef_phase phase; - int bandwidth_version; - int lowpass_version; - double coef_scale; - bool synthesis; // Synthesis is supported - bool multirate; -}; - -// Like std::fill, but returns the end iterator - -template -I fill(I b, I e, T v) { - std::fill(b, e, v); - return e; -} - -// Multiply a vector by a scalar, in-place. -// Used only at the setup stage, so performance is not critical. - -template -void scale_vector(V &v, S s) { - for (auto &e: v) - e *= s; -} - -// Zero-padding source wrapper. This returns data from the underlying -// source within the interval src_i0 to src_i1, and zero elsewhere. - -template -struct zeropad_source { - typedef typename std::iterator_traits::value_type T; - zeropad_source(const S &source_, int64_t src_i0_, int64_t src_i1_): - source(source_), src_i0(src_i0_), src_i1(src_i1_) - { } - OI operator()(int64_t i0, int64_t i1, OI output) const { - int64_t overlap_begin = std::max(i0, src_i0); - int64_t overlap_end = std::min(i1, src_i1); - if (overlap_end <= overlap_begin) { - // No overlap - output = gaborator::fill(output, output + (i1 - i0), (T)0); - } else { - // Some overlap - if (overlap_begin != i0) { - output = gaborator::fill(output, output + (overlap_begin - i0), (T)0); - } - output = source(overlap_begin, overlap_end, output); - if (overlap_end != i1) { - output = gaborator::fill(output, output + (i1 - overlap_end), (T)0); - } - } - return output; - } - const S &source; - int64_t src_i0, src_i1; -}; - -template -struct pointer_source { - pointer_source(const T *p_, int64_t buf_i0_, int64_t buf_i1_): - p(p_), buf_i0(buf_i0_), buf_i1(buf_i1_) { } - T *operator()(int64_t i0, int64_t i1, T *output) const { - assert(i1 >= i0); - assert(i0 >= buf_i0); - assert(i1 <= buf_i1); - return std::copy(p + (i0 - buf_i0), p + (i1 - buf_i0), output); - } - const T *p; - int64_t buf_i0, buf_i1; -}; - -// Fill the buffer at dst, of length dstlen, with data from src where -// available, otherwise with zeroes. The data in src covers dst indices -// from i0 (inclusive) to i1 (exclusive). - -template -void copy_overlapping_zerofill(T *dst, size_t dstlen, const T *src, - int64_t i0, int64_t i1) -{ - pointer_source ps(src, i0, i1); - zeropad_source, T *> zs(ps, i0, i1); - zs(0, dstlen, dst); -} - -// Given a set of FFT coefficients "coefs" of a real sequence, where -// only positive-frequency coefficients (including DC and Nyquist) are -// valid, return the coefficient for an arbitrary frequency index "i" -// which may correspond to a negative frequency, or even an alias -// outside the range (0..fftsize-1). - -template -complex get_real_spectrum_coef(complex *coefs, int i, unsigned int fftsize) { - i &= fftsize - 1; - // Note that this is >, not >=, becase fs/2 is considered nonnegative - bool neg_fq = (i > (int)(fftsize >> 1)); - if (neg_fq) { - i = fftsize - i; - } - complex c = coefs[i]; - if (neg_fq) { - c = conj(c); - } - return c; -} - -// A set of buffers of various sizes used for temporary vectors during -// analysis. These are allocated as a single block to reduce the -// number of dynamic memory allocations. - -template -struct buffers { - static const size_t maxbufs = 10; - typedef complex C; - buffers(unsigned int fftsize_max, - unsigned int sftsize_max): - n(0) - { - offset[0] = 0; - // Define the size of each buffer - def(fftsize_max * sizeof(C)); // 0 - def(fftsize_max * sizeof(C)); // 1 - def(sftsize_max * sizeof(C)); // 2 - def(sftsize_max * sizeof(C)); // 3 - def(sftsize_max * sizeof(C)); // 4 - def(fftsize_max * sizeof(T)); // 5 - assert(n <= maxbufs); - data = ::operator new(offset[n]); - } - ~buffers() { - ::operator delete(data); - } - void def(size_t size) { - size_t o = offset[n++]; - offset[n] = o + size; - } - // A single buffer of element type E - template - struct buffer { - typedef E *iterator; - buffer(void *b_, void *e_): - b((E *)b_), e((E *)e_) - { } - iterator begin() const { return b; } - iterator end() const { return e; } - E *data() { return b; } - const E *data() const { return b; } - E &operator[](size_t i) { return b[i]; } - const E &operator[](size_t i) const { return b[i]; } - size_t size() const { return e - b; } - private: - E *b; - E *e; - }; - // Get buffer number "i" as a vector-like object with element type "E" - // and a length of "len" elements. - template - buffer get(size_t i, size_t len) { - len *= sizeof(E); - size_t o = offset[i]; - assert(len <= offset[i + 1] - o); - return buffer((char *)data + o, (char *)data + o + len); - } -private: - void *data; - size_t n; - size_t offset[maxbufs + 1]; -}; - -// Get the bounds of the range of existing coefficients for a -// given band, in units of coefficient samples. - -template -void get_band_coef_bounds(const coefs &msc, int oct, unsigned int obno, - coef_index_t &ci0_ret, coef_index_t &ci1_ret) -{ - const sliced_coefs &sc = msc.octaves[oct]; - const typename sliced_coefs::slices_t &slices = sc.slices; - if (slices.empty()) { - // Don't try to convert int64t_min/max slices to coef time - ci0_ret = 0; - ci1_ret = 0; - return; - } - // Convert from slices to coefficient samples - ci0_ret = coef_time(*sc.meta, slices.begin_index(), oct, obno); - ci1_ret = coef_time(*sc.meta, slices.end_index(), oct, obno); -} - -template -void get_band_coef_bounds(const coefs &msc, int gbno, - coef_index_t &ci0_ret, coef_index_t &ci1_ret) -{ - int oct; - unsigned int obno; // Band number within octave - bool r = gaborator::bno_split(*msc.meta, gbno, oct, obno, true); - assert(r); - get_band_coef_bounds(msc, oct, obno, ci0_ret, ci1_ret); -} - -// Evaluate the frequency-domain analysis filter kernel of band "bp" -// at frequency "ff" - -template -double eval_kernel(parameters *, band_params *bp, double ff) { - if (bp->dc) { - return gaussian_windowed_lowpass_1(bp->ff, bp->ffsd, ff); - } else { - return norm_gaussian(bp->ffsd, ff - bp->ff); - } -} - -// Evaluate the frequency-domain synthesis filter kernel of band "bp" -// at frequency "ff" - -template -double eval_dual_kernel(parameters *params, band_params *bp, double ff) { - double gain = 1.0; - if (params->lowpass_version == 2) { - if (bp->dc) { - // Adjust the gain of the reconstruction lowpass - // filter to make the overall gain similar to to - // the bandpass region. - double adjusted_overlap = params->sd() / - (log(2) / params->bands_per_octave); - double avg_bandpass_gain = adjusted_overlap * sqrt(M_PI); - gain = avg_bandpass_gain * 0.5; - } - } - return eval_kernel(params, bp, ff) * gain; -} - -template -struct analyzer: public refcounted { - typedef complex C; - struct plan; - - analyzer(const parameters ¶ms_): - params(params_), - max_step_log2(0), - fftsize_max(0), - sftsize_max(0) - { - // Sanity check - assert(params.ff_min < 0.5); - - band_spacing_log2 = params.band_spacing_log2(); - band_spacing = params.band_spacing(); - - // The tuning adjustment, as a log2ff. This is a number between - // 0 and band_spacing_log2, corresponding to a frequency at - // or slightly above the sampling frequency where a band - // center would fall if they actually went that high. - // Tuning is done by increasing the center frequencies of - // all bands by this amount relative to the untuned case - // where one band would fall on fs exactly. - tuning_log2ff = sane_fmod(log2(params.ff_ref), band_spacing_log2); - - // Calculate the total number of bands needed so that - // the lowest band has a frequency <= params.ff_min. - // end_log2ff = the log2ff of the band after the last (just past fs/2); - // the -1 below is the log of the /2 above - double end_log2ff = tuning_log2ff - 1; - n_bandpass_bands_total = - (unsigned int)ceil((end_log2ff - log2(params.ff_min)) / - band_spacing_log2); - n_bands_total = n_bandpass_bands_total + 1; - - top_band_log2ff = end_log2ff - band_spacing_log2; - - ffref_gbno = (int)rint((top_band_log2ff - log2(params.ff_ref)) / - band_spacing_log2); - - // Establish affine transforms for converting between - // log-frequencies (log2(ff)) and bandpass band numbers. - // Derivation: - //ff = exp2(tuning_log2ff - 1 - (gbno + 1) * band_spacing_log2) - //log2(ff) = tuning_log2ff - 1 - (gbno + 1) * band_spacing_log2 - //tuning_log2ff - 1 - (gbno + 1) * band_spacing_log2 = log2(ff) - //-(gbno + 1) * band_spacing_log2 = log2(ff) - tuning_log2ff + 1 - //-(gbno + 1) * band_spacing_log2 = log2(ff) - tuning_log2ff + 1 - //-(gbno + 1) = (log2(ff) - tuning_log2ff + 1) / band_spacing_log2 - //-gbno - 1 = (log2(ff) - tuning_log2ff + 1) / band_spacing_log2 - //-gbno = ((log2(ff) - tuning_log2ff + 1) / band_spacing_log2) + 1 - //gbno = -(((log2(ff) - tuning_log2ff + 1) / band_spacing_log2) + 1) - //gbno = a log2(ff) + b, - // where a = -1 / band_spacing_log2 = -params.bands_per_octave - // and b = -a * tuning_log2ff + a - 1 - // The cast to double is necessary because we can't take the - // negative of an unsigned int. - double a = -(double)params.bands_per_octave; - double b = -a * tuning_log2ff + a - 1; - log2ff_bandpass_band = affine_transform(a, b); - bandpass_band_log2ff = log2ff_bandpass_band.inverse(); - - { - // Precalculate the parameters of the downsampling filter. - // These are the same for all plans, and need to be - // calculated before creating the plans; in particular, we - // need to know the support before we can create the - // plans, because in low-bpo cases, it can determine the - // minimum amount of fat needed. The filter kernel is - // specific to the plan as it depends on the FFT size, - // and will be calculated later. - - // When operating at a high Q, the we will need to use - // large FFTs in any case, and it makes sense to use a - // narrow transition band because we can get that - // essentially for free, and the passband will be - // correspondingly wider, which will allow processing more - // bands at the lower sample rate. Conversely, at low - // Q, we should use a wide transition band so that the - // FFTs can be kept short. - - // The filter is defined in terms of the lower - // (downsampled) sample rate. - - // Make the transition band the same width as the width - // (two-sided support) of a band at ff=0.25, but don't let - // the low edge go below 0.25 to make sure we have a - // reasonable amount of passband left. - double f1 = 0.5; - double f0 = - std::max(f1 - 2 * gaussian_support(ff_sd(0.25), params.max_error), - 0.25); - assert(f0 < f1); - - // The cutoff frequency is in the center of the transition band - double ff = (f0 + f1) * 0.5; - double support = (f1 - f0) * 0.5; - double ff_sd = gaussian_support_inv(support, params.max_error); - - // Calculate and save the time-domain support of the - // downsampling lowpass filter for use in analyze_sliced(). - double time_sd = sd_f2t(ff_sd); - - // Set members - ds_passband = f0; - ds_ff = ff; - ds_ff_sd = ff_sd; - // Since the filter is designed at the lower sample rate, - // ds_time_support is in the unit of lower octave samples - ds_time_support = gaussian_support(time_sd, params.max_error * 0.5); - } - - // Determine the octave structure, packing each band into the - // lowest octave possible. For now, while bpo is restricted - // to integer values, this just means determining how many - // bands need to go in the top octave, and the remaining ones - // will be divided into groups of bpo bands (except possibly - // the last). - int gbno; - for (gbno = bandpass_bands_begin(); gbno < bandpass_bands_end(); gbno++) { - double ff = bandpass_band_ff(gbno); - double ffsd = ff_sd(ff); - double ff_support = gaussian_support(ffsd, params.max_error * 0.5); - // If the bandpass support falls within the downsampling filter - // passband of the next octave, we can switch octaves. - if (params.multirate && ff + ff_support <= ds_passband / 2) - break; - } - n_bands_top_octave = gbno; - - // Figure out the number of octaves, keeping in mind that - // the top octave is of variable size, and DC band is added - // to the bottom octave even if that makes it larger than - // the others. - n_octaves = 1 // The top octave - + (n_bandpass_bands_total - n_bands_top_octave + - (params.bands_per_octave - 1)) / params.bands_per_octave; - - // Calculate the kernel support needed for the lowest-frequency - // bandpass band to use as a basis for an initial estimate of - // the FFT size needed. This duplicates some code at the - // beginning of make_band(). - assert(n_bands_top_octave >= 1); - int low_bp_band = n_bands_top_octave - 1; - double low_bp_band_time_sd = time_sd(band_ff(low_bp_band)); - - double low_bp_band_time_analysis_support = - gaussian_support(low_bp_band_time_sd, params.max_error); - double low_bp_band_time_synthesis_support = - low_bp_band_time_analysis_support * synthesis_support_multiplier(); - - make_zones(); - - // Make analysis plans - // Since ds_time_support is in the unit of lower octave samples, - // we need to multiply it by two to get upper octave samples. - unsigned int max_support = - std::max(ceil(low_bp_band_time_analysis_support), - ds_time_support * 2); - unsigned int size = next_power_of_two(max_support * 2); - ref p; - for (;;) { - p = new plan(this, false, size, max_support); - if (p->ok) - break; - size *= 2; - } - anl_plans.push_back(p); // Smallest possible plan - p = new plan(this, false, size * 2, max_support); - assert(p->ok); - anl_plans.push_back(p); // Next larger plan - - if (params.synthesis) { - // Make synthesis plan (only one for now) - max_support = std::max(ceil(low_bp_band_time_synthesis_support), - ds_time_support * 2); - // Room for at at least the two fats + as much filet - size = next_power_of_two(max_support * 2) * 2; - p = new plan(this, true, size, max_support); - assert(p->ok); - syn_plans.push_back(p); - } - - for (int i = 0; i < (int)anl_plans.size(); i++) - make_band_plans(i, false); - for (int i = 0; i < (int)syn_plans.size(); i++) - make_band_plans(i, true); - - // Find the largest fftsize and sftsize of any plan - for (size_t i = 0; i < anl_plans.size(); i++) { - fftsize_max = std::max(fftsize_max, anl_plans[i]->fftsize); - sftsize_max = std::max(sftsize_max, anl_plans[i]->sftsize_max); - } - for (size_t i = 0; i < syn_plans.size(); i++) { - fftsize_max = std::max(fftsize_max, syn_plans[i]->fftsize); - sftsize_max = std::max(sftsize_max, syn_plans[i]->sftsize_max); - } - - // Lay out the coefficient structures according to the - // synthesis plan if we have one, or the largest analysis - // plan if not. - std::vector> *cmeta_source = - params.synthesis ? &syn_plans : &anl_plans; - ref &largest_plan(((*cmeta_source)[cmeta_source->size() - 1])); - cmeta_any = make_meta(filet_part(largest_plan->fftsize)); - } - - void make_zones() { - // Band number starting at 0 close to fs/2 and increasing - // with decreasing frequency - int tbno = 0; - int oct = 0; - int zno = 0; - // Loop over the octaves, from high to low frequencies, - // creating new zones where needed - for (;;) { - int max_bands_this_octave = (zno == 0) ? - n_bands_top_octave : params.bands_per_octave; - int bands_remaining = n_bandpass_bands_total - tbno; - int bands_this_octave = std::min(max_bands_this_octave, bands_remaining); - int bands_below = bands_remaining - bands_this_octave; - bool dc_zone = (bands_below == 0); - bool dc_adjacent_zone = (bands_below < (int)params.bands_per_octave); - if (zno < 2 || dc_zone || dc_adjacent_zone || - params.bands_per_octave < 6) - { - make_zone(oct, zno, tbno, tbno + bands_this_octave, - dc_zone, bands_below); - zno++; - } - octaves.push_back(octave()); - octaves.back().z = zones[zno - 1].get(); - oct++; - tbno += bands_this_octave; - if (dc_zone) - break; - } - assert(octaves.size() == n_octaves); - } - - // Create a zone consisting of the bandpass bands band0 - // (inclusive) to band1 (exclusive), using the usual gbno - // numbering going from high to low frequencies, and - // possible a lowpass band band1. - - void make_zone(int oct, unsigned int zno, - int band0, int band1, - bool dc_zone, int bands_below) - { - assert(zones.size() == zno); - zone *z = new zone(); - z->zno = zno; - zones.push_back(ref>(z)); - - pod_vector power; - - // Create the real (non-mock) bands, from low to high - // frequency. - if (dc_zone) - // This zone has a lowpass band - push(z->bandparams, make_band(oct, band1, true, false)); - // The actual (non-mock) bandpass bands of this zone - for (int i = band1 - 1; i >= band0; i--) - push(z->bandparams, make_band(oct, i, false, false)); - - if (! dc_zone) { - // There are other zones below this; add mock bands - // to simulate them for purposes of calculating the dual. - - // Identify the lowest frequency of interest in the zone - assert(z->bandparams.size() >= 1); - band_params *low_band = z->bandparams[0].get(); - double zone_bottom_ff = low_band->ff - low_band->ff_support; - - int i = band1; - for (; i < band1 + bands_below; i++) { - band_params *mock_band = make_band(oct, i, false, true); - push(z->mock_bandparams, mock_band); - // There's no point in creating further mock bands - // once they no longer overlap with the current zone. - // The condition used here may cause the creation of - // one more mock band than is actually needed, as it - // is easier to create the band first and check for - // overlap later than the other way round. - if (mock_band->ff + mock_band->ff_support < zone_bottom_ff) { - i++; - break; - } - } - // Create a mock lowpass band. This may correspond to the - // actual lowpass band, or if the loop above exited early, - // it may merely be be a placeholder that serves no real - // purpose other than making the power vector look better. - push(z->mock_bandparams, make_band(oct, i, true, true)); - } - - // If there are other zones above this, add mock bands - // to simulate them for purposes of calculating the dual, - // but only up to the Nyquist frequency of the current - // octave. - int nyquist_band = oct * (int)params.bands_per_octave; - for (int i = band0 - 1; i >= nyquist_band; i--) - push(z->mock_bandparams, make_band(oct, i, false, true)); - - z->n_bands = (unsigned int)z->bandparams.size(); - - // Find the largest coefficient step in the zone, as this will - // determine the necessary alignment of signal slices in time, - // but make it at least two (corresponding max_step_log2 = 1) - // because the downsampling code requires alignement to even - // indices. - unsigned int m = 1; - for (unsigned int obno = 0; obno < z->bandparams.size(); obno++) { - m = std::max(m, z->bandparams[obno]->step_log2); - } - z->max_step_log2 = m; - - max_step_log2 = std::max(max_step_log2, m); - } - - // Calculate band parameters for a single band. - // - // If dc is true, this is the DC band, and gbno indicates - // the cutoff frequency; it is one more than the gbno of - // the lowest-frequency bandpass band. - - band_params * - make_band(int oct, double gbno, bool dc, bool mock) { - band_params *bp = new band_params; - if (dc) - // Make the actual DC band cutoff frequency a bit higher, - // by an empirically chosen fraction of a band, to reduce - // power fluctuations. - gbno -= 0.8750526596806952; - - // For bandpass bands, the center frequency, or for the - // lowpass band, the lowpass cutoff frequency, as a - // fractional frequency, in terms of the octave's sample - // rate. - double ff = ldexp(bandpass_band_ff(gbno), oct); - - // Standard deviation of the bandpass Gaussian, as a - // fractional frequency - double ffsd = ff_sd(ff); - double time_sd = sd_f2t(ffsd); - double time_support = - gaussian_support(time_sd, params.max_error); - - // The support of the Gaussian, i.e., the smallest standard - // deviation at which it can be truncated on each side - // without the error exceeding our part of the error budget, - // which is some fraction of params.max_error. Note - // that this is one-sided; the full width of the support - // is 2 * ff_support. - double bp_ff_support = gaussian_support(ffsd, params.max_error * 0.5); - // Additional support for the flat portion of the DC band lowpass - double dc_support = dc ? ff : 0; - // Total frequency-domain support for this band, one-sided - double band_support = bp_ff_support + dc_support; - // Total support needed for this band, two-sided - double band_2support = band_support * 2; - - // Determine the downsampling factor for this band. - int exp = 0; - while (band_2support <= 0.5) { - band_2support *= 2; - exp++; - } - bp->step_log2 = exp; - bp->step = 1U << bp->step_log2; - - bp->dc = dc; - bp->ff = ff; - bp->ff_support = band_support; - bp->time_support = time_support; - bp->ffsd = ffsd; - - return bp; - } - - // Given a fractional frequency, return the standard deviation - // of the frequency-domain window as a fractional frequency - double ff_sd(double ff) const { return params.sd() * ff; } - - // Given a fractional frequency, return the standard deviation - // of the time-domain window in samples. - // - // ff_sd = 1.0 / (tau * t_sd) - // per http://users.ece.gatech.edu/mrichard/ - // Gaussian%20FT%20and%20random%20process.pdf - // and python test program gaussian-overlap.py - // => (tau * t_sd) * ff_sd = 1.0 - // => t_sd = 1.0 / (tau * f_sd) - double time_sd(double ff) const { return 1.0 / (tau * ff_sd(ff)); } - - double q() const { return params.q(); } - - // Find the worst-case time support of the analysis filters, i.e., - // the largest distance in time between a signal sample and a - // coefficient affected by that sample. - double analysis_support() const { - // The lowpass filter is the steepest one - return analysis_support(band_lowpass()); - } - - // Find the time support of the analysis filter for bandpass band - // number gbno - double analysis_support(double gbno) const { - return gaussian_support(time_sd(bandpass_band_ff(gbno)), - params.max_error); - } - - // Ditto for the resynthesis filters. - double synthesis_support() const { - return analysis_support() * synthesis_support_multiplier(); - } - - double synthesis_support(double gbno) const { - return analysis_support(gbno) * synthesis_support_multiplier(); - } - - // The synthesis support multiplier, the factor by which the - // synthesis filters are wider than the analysis filters in - // the time domain. - double synthesis_support_multiplier() const { - if (! params.synthesis) - return 1.0; - // At high ff_min, the top band may be the one with the widest support. - if (params.ff_min > 0.25) - return 4; - if (params.bands_per_octave <= 4) - return 2.5; - return 2.3; - } - - // Get the center frequency of bandpass band "gbno", which - // need not be a valid bandpass band number; out-of-range - // arguments will return extrapolated frequencies based on - // the logarithmic band spacing. - double bandpass_band_ff(double gbno) const { - return exp2(bandpass_band_log2ff(gbno)); - } - - // Get the band number of the bandpass band corresponding - // to the fractional frequency "ff", as a floating point - // number. This is the inverse of bandpass_band_ff(). - double ff_bandpass_band(double ff) const { - return log2ff_bandpass_band(log2(ff)); - } - - int choose_plan(const std::vector> &plans, int64_t size) const - { - unsigned int i = 0; - while (i < plans.size() - 1 && plans[i]->filet_size < size) - i++; - return i; - } - - void - synthesize_one_slice(int oct, int pno, const coefs &msc, - const pod_vector &downsampled, - sample_index_t t0, - T *signal_out, - pod_vector &buf0, // fftsize - pod_vector &buf2, // largest sftsize - pod_vector &buf3 // largest sftsize - ) const - { - const plan &plan(*syn_plans[pno]); - zone &z = *octaves[oct].z; - pod_vector &signal(buf0); - std::fill(signal.begin(), signal.end(), (T)0); - - pod_vector &coefbuf(buf3); - - for (unsigned int obno = 0; obno < z.bandparams.size(); obno++) { - band_params *bp = z.bandparams[obno].get(); - band_plan *bpl = &bp->syn_plans[pno]; - - // log2 of the coefficient downsampling factor - int coef_shift = bp->step_log2; - coef_index_t ii = t0 >> coef_shift; - - read(msc, bno_merge(oct, obno), ii, ii + bpl->sftsize, coefbuf.data()); - C *indata = coefbuf.data(); - - pod_vector &sdata(buf2); - - T scale_factor = (T)1 / ((T)params.coef_scale * bpl->sftsize); - - // Apply phase correction, adjust for non-integer center - // frequency, and apply scale factor. Note that phase - // must be calculated in double precision. - - // We can't use bp->ff here because in the case of the - // lowpass band, it's the cutoff rather than the center. - double ff = bpl->center * plan.inv_fftsize_double; - double arg = (params.phase == coef_phase::global) ? - tau * t0 * ff : 0; - C phase_times_scale = C(cos(arg), sin(arg)) * scale_factor; - elementwise_product_times_scalar(sdata.data(), indata, - bpl->shift_kernel_conj.data(), - phase_times_scale, bpl->sftsize); - - // Switch to frequency domain - bpl->sft->transform(sdata.data()); - - // Multiply signal spectrum by frequency-domain dual window, - // accumulating result in signal. - - for (unsigned int i = 0; i < bpl->sftsize; i++) { - int iii = (bpl->fq_offset_int + i) & (plan.fftsize - 1); - // Note the ifftshift of the input index, as f=0 - // appears in the middle of the window - C v = sdata[i ^ (bpl->sftsize >> 1)] * bpl->dual_kernel[i]; - // Frequency symmetry - signal[iii] += v; - if (params.lowpass_version == 2 || ! bp->dc) - signal[-iii & (plan.fftsize - 1)] += conj(v); - } - } - - if (oct + 1 < (int) n_octaves) { - // Upsample the downsampled data from the lower octaves - pod_vector &sdata(buf2); - assert(downsampled.size() == plan.dsparams.sftsize); - assert(sdata.size() >= plan.dsparams.sftsize); -#if GABORATOR_USE_REAL_FFT - plan.dsparams.rsft->transform(downsampled.data(), sdata.begin()); -#else - // Real to complex - std::copy(downsampled.begin(), downsampled.end(), sdata.begin()); - plan.dsparams.sft->transform(sdata.data()); -#endif - for (unsigned int i = 0; i < plan.dsparams.sftsize; i++) { - sdata[i] *= plan.dsparams.dual_kernel - [i ^ (plan.dsparams.sftsize >> 1)]; - } - - // This implicitly zero pads the spectrum, by not adding - // anything to the middle part. The splitting of the - // Nyquist band is per http://dsp.stackexchange.com/ - // questions/14919/upsample-data-using-ffts-how-is-this- - // exactly-done but should not really matter because there - // should be no energy there to speak of thanks to the - // windowing above. - - assert(plan.dsparams.sftsize == plan.fftsize / 2); - unsigned int i; - for (i = 0; i < plan.dsparams.sftsize / 2; i++) - signal[i] += sdata[i]; - //C nyquist = sdata[i] * (T)0.5; - C nyquist = sdata[i] * (T)0.5; - signal[i] += nyquist; - signal[i + plan.fftsize / 2] += nyquist; - i++; - for (;i < plan.dsparams.sftsize; i++) - signal[i + plan.fftsize / 2] += sdata[i]; - } - - // Switch to time domain -#if GABORATOR_USE_REAL_FFT - plan.rft->itransform(signal.data(), signal_out); -#else - plan.ft->itransform(signal.data()); - // Copy real part to output - complex2real(signal.begin(), signal.end(), signal_out); -#endif - } - -private: - - // Analyze a signal segment consisting of any number of samples. - // oct is the octave; this is 0 except in recursive calls - // real_signal points to the first sample - // t0 is the sample time of the first sample - // t1 is the sample time of the sample after the last sample - // coefficients are added to msc - - void - analyze_sliced(buffers &buf, int oct, const T *real_signal, - sample_index_t t0, sample_index_t t1, - double included_ds_support, - coefs &msc) const - { - assert(t1 >= t0); - int pno = choose_plan(anl_plans, t1 - t0); - const plan &plan(*anl_plans[pno]); - - // Even though we don't align the FFTs to full filet-size - // slices in this code path, we still need to align them to - // coefficient samples so that we don't have to do expensive - // sub-sample time shifts. Specifically, we need to align - // them to the largest coefficient time step of the octave. - // slice_t0 is the sample time of the first sample in the - // filet (not the FFT as a whole). - zone &z = *octaves[oct].z; - sample_index_t slice_t0 = t0 & ~((1 << z.max_step_log2) - 1); - //printf("slice t0 = %d\n", (int)slice_t0); - - // Find the number of slices we need to divide the signal into. - // We need it ahead of time so that we can size the "downsampled" - // array accordingly. - uint64_t n_slices = ((t1 - slice_t0) + (plan.filet_size - 1)) / plan.filet_size; - - // The "downsampled" array needs to fit one filet per slice + - // the fat on each side, all of half size thanks to downsampling. - // Length of each downsampled slice (including padding) - uint64_t dstotlen = ((n_slices * plan.filet_size) + (2 * plan.fat_size)) >> 1; - - // The range of sample times covered by the "downsampled" array - sample_index_t tmp = slice_t0 - (int)plan.fat_size; - assert((tmp & 1) == 0); - sample_index_t dst0 = tmp >> 1; - sample_index_t dst1 = dst0 + (int64_t)dstotlen; - - // Not all of the "downsampled" array actually contains - // nonzero data. Calculate adjusted bounds to use in the - // recursive analysis so that we don't needlessly analyze - // zeroes. - int ds_support = (int)ceil(ds_time_support - included_ds_support); - sample_index_t dst0a = std::max(dst0, (t0 >> 1) - ds_support); - sample_index_t dst1a = std::min(dst1, (t1 >> 1) + 1 + ds_support); - - // Buffer for the downsampled signal. Since the size depends - // on the total amount of signal analyzed in this call (being - // about half of it), it can't be preallocated, but has to be - // dynamically allocated in each call. - pod_vector downsampled(dstotlen); - // "downsampled" will be added to, not assigned to, so we need - // to set it to zero initially. - std::fill(downsampled.begin(), downsampled.end(), 0); - - auto slice(buf.template get(5, plan.fftsize)); - // Clear the fat on both ends (once) - std::fill(slice.data(), slice.data() + plan.fat_size, 0); - std::fill(slice.data() + slice.size() - plan.fat_size, - slice.data() + slice.size(), 0); - - // For each slice. Note that slice_i counts from 0, not from - // the slice index of the first slice. - for (uint64_t slice_i = 0; slice_i < n_slices; slice_i++) { - if (slice_t0 >= t1) - break; - sample_index_t slice_t1 = std::min(slice_t0 + plan.filet_size, t1); - // Copy into filet part of aligned buffer, possibly zero padding - // if the remaining signal is shorter than a full slice. - copy_overlapping_zerofill(slice.data() + (int)plan.fat_size, - plan.filet_size, - real_signal, - t0 - slice_t0, - t1 - slice_t0); - // Analyze the slice - auto spectrum(buf.template get(1, plan.fftsize)); -#if GABORATOR_USE_REAL_FFT - plan.rft->transform(slice.data(), spectrum.data()); -#else - // Real to complex - auto signal(buf.template get(0, plan.fftsize)); - std::copy(slice.data(), slice.data() + plan.fftsize, signal.begin()); - plan.ft->transform(signal.data(), spectrum.data()); -#endif - auto tmp(buf.template get(2, plan.sftsize_max)); - auto coefbuf(buf.template get(4, plan.sftsize_max)); - - T scale_factor = (T)params.coef_scale * plan.inv_fftsize_t; - - for (unsigned int obno = 0; obno < z.bandparams.size(); obno++) { - band_params *bp = z.bandparams[obno].get(); - band_plan *bpl = &bp->anl_plans[pno]; - C *sdata = tmp.data(); - - // Multiply a slice of the spectrum by the frequency- - // domain window and store in sdata. - // - // We need to take care not to overrun the beginning or - // end of the spectrum - for the dc band, we always - // need to wrap around to negative frequencies, and - // potentially it could happen with other bands, too, - // if they are really wide. To avoid the overhead of - // checking in the inner loop, use a separate slow path - // for the rare cases where wrapping happens. - - int start_index = bpl->fq_offset_int; - int end_index = bpl->fq_offset_int + bpl->sftsize; - if (start_index >= 0 && end_index < (int)((plan.fftsize >> 1) + 1)) { - // Fast path: the slice lies entirely within the - // positive-frequency half of the spectrum (including - // DC and Nyquist). - elementwise_product(sdata, - spectrum.data() + start_index, - bpl->kernel.data(), - bpl->sftsize); - } else { - // Slow path - for (size_t i = 0; i < bpl->sftsize; i++) - sdata[i] = get_real_spectrum_coef(spectrum.data(), - (int)(start_index + i), plan.fftsize) * bpl->kernel[i]; - } - // The band center frequency is at the center of the - // spectrum slice and at the center of the window, so - // it also ends up at the center of sdata. The center - // frequency of the band is considered f=0, so for the - // ifft, it should be at index 0, not the center. - // Therefore, in principle we should perform an - // ifftshift of sdata here before the ifft, but since - // the time-domain data are going to be multiplied by - // the shift kernel anyway, the ifftshift is baked - // into the shift kernel by flipping the sign of every - // other element so that it is effectively free. - - // Switch to time domain - auto band(buf.template get(3, plan.sftsize_max)); - bpl->sft->itransform(sdata, band.data()); - - // Apply ifftshift, adjust for non-integer center - // frequency, correct phase, scale amplitude, and add - // to the output coefficients. - double ff = bpl->center * plan.inv_fftsize_double; - double arg; - if (params.phase == coef_phase::global) - arg = -tau * (slice_t0 - plan.fat_size) * ff; - else - arg = 0; - C phase_times_scale = C(cos(arg), sin(arg)) * scale_factor; - elementwise_product_times_scalar(coefbuf.data(), band.data(), - bpl->shift_kernel.data(), - phase_times_scale, - bpl->sftsize); - // log2 of the coefficient downsampling factor - int coef_shift = bp->step_log2; - assert(((slice_t0 - (int)plan.fat_size) & ((1 << coef_shift) - 1)) == 0); - coef_index_t ii = (slice_t0 - (int)plan.fat_size) >> coef_shift; - - // Only part of coefbuf contains substantially nonzero - // data: that corresponding to the signal interval - // t0..t1 + the actual support of the filter for this band. - // There's no point adding the zeros to the coefficients, - // so trim. - int support = (int)ceil(bp->time_support); - coef_index_t ii0 = std::max(ii, (t0 - support) >> coef_shift); - coef_index_t ii1 = std::min(ii + bpl->sftsize, - ((t1 + support) >> coef_shift) + 1); - add(msc, bno_merge(oct, obno), ii0, ii1, coefbuf.data() + (ii0 - ii)); - } - - // Downsample - if (oct + 1 < (int) n_octaves) { - T *downsampled_dst = downsampled.data() + - slice_i * (plan.filet_size >> 1); - - auto sdata(buf.template get(2, plan.sftsize_max)); - // This is using a larger buffer than we actually need - auto ddata(buf.template get(0, plan.sftsize_max)); - assert(ddata.size() >= plan.dsparams.sftsize); - // Extract the low-frequency part of "spectrum" into "sdata" - // and multiply it by the lowpass filter frequency response. - // This means both positive and negative low frequencies. - size_t half_size = plan.dsparams.sftsize >> 1; - assert(plan.fftsize - half_size == 3 * half_size); -#if GABORATOR_USE_REAL_FFT - // Positive frequencies - elementwise_product(sdata.data(), spectrum.data(), - plan.dsparams.kernel.data() + half_size, - half_size); - // Nyquist - sdata[half_size] = 0; - // Use the same buffer as the complex FFT, but as floats - T *real_ddata = reinterpret_cast(ddata.data()); - plan.dsparams.rsft->itransform(sdata.data(), real_ddata); - // Only accumulate nonzero part - size_t n = ((slice_t1 - slice_t0) + 2 * (int)plan.fat_size + 1) >> 1; - for (size_t i = 0; i < n; i++) - downsampled_dst[i] += real_ddata[i]; -#else - // Positive frequencies - elementwise_product(sdata.data(), spectrum.data(), - plan.dsparams.kernel.data() + half_size, - half_size); - // Negative requencies - elementwise_product(sdata.data() + half_size, - spectrum.data() + plan.fftsize - half_size, - plan.dsparams.kernel.data(), half_size); - // Convert to time domain - plan.dsparams.sft->itransform(sdata.data(), ddata.data()); - for (unsigned int i = 0; i < plan.dsparams.sftsize; i++) - downsampled_dst[i] += ddata[i].real(); -#endif - } - - // Next slice - slice_t0 = slice_t1; - } - - // Recurse - if (oct + 1 < (int)n_octaves) - analyze_sliced(buf, oct + 1, downsampled.data() + (dst0a - dst0), - dst0a, dst1a, ds_time_support / 2, msc); - } - - // Resynthesize audio from the coefficients in "msc". The audio will - // cover samples from t0 (inclusive) to t1 (exclusive), and is stored - // starting at *real_signal, which must have room for (t1 - t0) - // samples. The octave "oct" is 0 except in recursive calls. - - void - synthesize_sliced(int oct, const coefs &msc, - sample_index_t t0, sample_index_t t1, - T *real_signal) const - { - int pno = choose_plan(syn_plans, t1 - t0); - const plan &plan(*syn_plans[pno]); - - // XXX clean up - no need to pass support arg - slice_index_t si0 = plan.affected_slice_b(t0, plan.oct_support); - slice_index_t si1 = plan.affected_slice_e(t1, plan.oct_support); - - // sub_signal holds the reconstructed subsampled signal from - // the lower octaves, for the entire time interval covered by - // the slices - int sub_signal_len = ((si1 - si0) * plan.filet_size + 2 * plan.fat_size) / 2; - pod_vector sub_signal(sub_signal_len); - std::fill(sub_signal.begin(), sub_signal.end(), 0); - if (oct + 1 < (int)n_octaves) { - int64_t sub_t0 = si0 * (plan.filet_size / 2); - int64_t sub_t1 = sub_t0 + sub_signal_len; - // Recurse - assert(sub_t1 - sub_t0 == (int64_t)sub_signal.size()); - synthesize_sliced(oct + 1, msc, sub_t0, sub_t1, sub_signal.data()); - } - - // Allocate buffers for synthesize_one_slice(), to be shared - // between successive calls to avoid repeated allocation - pod_vector buf0(plan.fftsize); - //pod_vector buf1(fftsize); - pod_vector buf2(plan.sftsize_max); - pod_vector buf3(plan.sftsize_max); - - pod_vector downsampled(plan.dsparams.sftsize); - pod_vector signal_slice(plan.fftsize); - - // For each slice - for (slice_index_t si = si0; si < si1; si++) { - sample_index_t slice_t0 = si * plan.filet_size; - - // Copy downsampled signal to "downsampled" for upsampling - if (oct + 1 < (int) n_octaves) { - int bi = (si - si0) * filet_part(plan.dsparams.sftsize); - int ei = bi + plan.dsparams.sftsize; - assert(bi >= 0); - assert(ei <= (int)sub_signal.size()); - std::copy(sub_signal.begin() + bi, - sub_signal.begin() + ei, - downsampled.begin()); - } - - synthesize_one_slice(oct, pno, msc, downsampled, slice_t0, - signal_slice.data(), buf0, buf2, buf3); - - // Copy overlapping part - sample_index_t b = std::max(slice_t0 + plan.fat_size, t0); - sample_index_t e = std::min(slice_t0 + plan.fftsize - plan.fat_size, t1); - for (sample_index_t i = b; i < e; i++) - real_signal[i - t0] = signal_slice[i - slice_t0]; - } - } - -public: - // The main analysis entry point. - // The resulting coefficients are added to any existing - // coefficients in "msc". - - void analyze(const T *real_signal, sample_index_t t0, sample_index_t t1, - coefs &msc, int n_threads = 1) const - { - analyze1(real_signal, t0, t1, msc, n_threads, 1); - } - - void analyze1(const T *real_signal, sample_index_t t0, sample_index_t t1, - coefs &msc, int n_threads, int level) const - { - assert(msc.octaves.size() == n_octaves); - (void)n_threads; - buffers buf(fftsize_max, sftsize_max); - analyze_sliced(buf, 0, real_signal, t0, t1, 0, msc); - } - - // The main synthesis entry point - - void - synthesize(const coefs &msc, sample_index_t t0, sample_index_t t1, - T *real_signal, int n_threads = 1) const - { - assert(params.synthesis); - (void)n_threads; - synthesize_sliced(0, msc, t0, t1, real_signal); - } - - bool bno_split(int gbno, int &oct, unsigned int &obno, bool dc) const { - return gaborator::bno_split(*cmeta_any, gbno, oct, obno, dc); - } - - int bno_merge(int oct, unsigned int obno) const { - return gaborator::bno_merge(*cmeta_any, oct, obno); - } - - // Get the bounds of the range of existing coefficients for all bands, - // in units of signal samples. - void get_coef_bounds(const coefs &msc, - sample_index_t &si0_ret, sample_index_t &si1_ret) - const - { - // The greatest coefficient range typically occurs in the - // lowest bandpass band, but this is not always the case, - // so to be certain, check them all. - sample_index_t min_si0 = INT64_MAX; - sample_index_t max_si1 = INT64_MIN; - for (int band = bands_begin(); band != bands_end(); band++) { - coef_index_t ci0, ci1; - get_band_coef_bounds(msc, band, ci0, ci1); - // Convert from coefficient samples to signal samples - int exp = band_scale_exp(band); - sample_index_t si0 = shift_left(ci0, exp); - sample_index_t si1 = shift_left(ci1 - 1, exp) + 1; - min_si0 = std::min(min_si0, si0); - max_si1 = std::max(max_si1, si1); - } - si0_ret = min_si0; - si1_ret = max_si1; - } - - unsigned int band_step_log2(int gbno) const { - return gaborator::band_step_log2(*cmeta_any, gbno); - } - - int bandpass_bands_begin() const { return 0; } - int bandpass_bands_end() const { return n_bands_total - 1; } - - int bands_begin() const { return 0; } - int bands_end() const { return n_bands_total; } - - int band_lowpass() const { return n_bands_total - 1; } - int band_ref() const { return ffref_gbno; } - - // Get the center frequency of band number gbno as a fractional - // frequency. gbno must be a valid band number. For the lowpass - // band, this returns zero. - double band_ff(int gbno) const { - if (gbno == band_lowpass()) - return 0; - return bandpass_band_ff(gbno); - } - - ~analyzer() { - } - - // Get the base 2 logarithm of the downsampling factor of - // band "obno" in octave "oct" - int band_scale_exp(int oct, unsigned int obno) const { - return gaborator::band_scale_exp(*cmeta_any->octaves[oct].z, oct, obno); - } - - // Get the base 2 logarithm of the downsampling factor of - // band "gbno" - int band_scale_exp(int gbno) const { - int oct; - unsigned int obno; // Band number within octave - bool r = bno_split(gbno, oct, obno, true); - assert(r); - return band_scale_exp(oct, obno); - } - - // Get the base 2 logarithm of the highest downsampling factor of - // any band - int band_scale_exp_max() const { - return band_scale_exp(bandpass_bands_end() - 1); - } - - - // Find the sample time of the band "gbno" coefficient closest to - // time "t". "gbno" must be a valid band number. - sample_index_t nearest_coef_sample(int gbno, double t) const { - int shift = band_step_log2(gbno); - return shift_left((sample_index_t) round(ldexp(t, -shift)), shift); - } - // Find the highest coefficient sample time less than or equal to - // "t" for band "gbno". "gbno" must be a valid band number. - sample_index_t floor_coef_sample(int gbno, double t) const { - int shift = band_step_log2(gbno); - return shift_left((sample_index_t) floor(ldexp(t, -shift)), shift); - } - // Find the lowestt coefficient sample time greater than or equal - // to "t" for band "gbno". "gbno" must be a valid band number. - sample_index_t ceil_coef_sample(int gbno, double t) const { - int shift = band_step_log2(gbno); - return shift_left((sample_index_t) ceil(ldexp(t, -shift)), shift); - } - - // Members initialized in the constructor, and listed in - // order of initialization - parameters params; - double band_spacing_log2; - double band_spacing; - double tuning_log2ff; - affine_transform log2ff_bandpass_band; - affine_transform bandpass_band_log2ff; - unsigned int n_bandpass_bands_total; - unsigned int n_bands_top_octave; - - struct plan: public refcounted { - plan(const plan &) = delete; - plan(analyzer *anl, bool synthesis_, unsigned int fftsize_, double support_): - ok(false), - synthesis(synthesis_), - fftsize(fftsize_), - oct_support(support_), - sftsize_max(0) - { - fftsize_log2 = whichp2(fftsize); - - inv_fftsize_double = 1.0 / fftsize; - inv_fftsize_t = (T) inv_fftsize_double; - -#if GABORATOR_USE_REAL_FFT - rft = pool, int>::shared.get(fftsize); -#else - ft = pool, int>::shared.get(fftsize); -#endif - // Set up the downsampling parameters in dsparams. - - // Downsampling is always by a factor of two. - // dsparams.sftsize is the size of the FFT used to go back to - // the time domain after discarding the top half of the - // spectrum. - dsparams.sftsize = fftsize >> 1; - dsparams.kernel.resize(dsparams.sftsize); - if (synthesis) - dsparams.dual_kernel.resize(dsparams.sftsize); - - // Use the convolution of a rectangle and a Gaussian. - // A piecewise function composed from two half-gaussians - // joined by a horizontal y=1 segment is not quite smooth - // enough. Put the passband in the middle. - for (int i = 0; i < (int)dsparams.sftsize; i++) - dsparams.kernel[i] = - gaussian_windowed_lowpass_1(anl->ds_ff, anl->ds_ff_sd, - ((double)i / dsparams.sftsize) - 0.5); - - if (synthesis) { - // The dual_kernel field of the downsampling pseudo-band holds - // the upsampling filter, identical to the downsampling filter - // except for amplitude scaling. - std::copy(dsparams.kernel.begin(), dsparams.kernel.end(), - dsparams.dual_kernel.begin()); - } - // Prescale the downsampling filter - scale_vector(dsparams.kernel, inv_fftsize_double); - if (synthesis) { - // Prescale the upsampling filter - scale_vector(dsparams.dual_kernel, 1.0 / dsparams.sftsize); - } -#if GABORATOR_USE_REAL_FFT - dsparams.rsft = pool, int>::shared.get(dsparams.sftsize); -#else - dsparams.sft = pool, int>::shared.get(dsparams.sftsize); -#endif - // It may be possible to reduce the size of the fat from 1/4 - // of the fftsize, but we need to keep things aligned with the - // coefficients, and it needs to be even for downsampling. - if (! synthesis) { - unsigned int align = 1 << std::max(anl->max_step_log2, 2U); - fat_size = (oct_support + (align - 1)) & ~(align - 1); - // There must be room for at least one signal sample in each - // half of the FFT; it can't be all fat - if (!(fat_size < (fftsize >> 1))) - return; // fail - } else { - fat_size = fftsize >> 2; - } - filet_size = fftsize - 2 * fat_size; - - // Constructor was successful - ok = true; - } - - // Index of first slice affected by sample at t0 - - // fft number i covers the sample range - // t = (i * filetsize .. i * filetsize + (fftsize - 1)) - // t >= i * filetsize and t < i * filetsize + fftsize - // A sample at t affects ffts i where - // i <= t / filetsize and - // i > (t - fftsize) / filetsize - // the filet of fft number i covers the sample range - // (fat + (i * filetsize) .. fat + (i * filetsize) + (filetsize - 1)) - // - // However, due to the FFT size being rounded up to a power of two, - // the outermost parts have near-zero weights and can be ignored; - // this is done by adjusting the time values by the width of that - // outermost part, which is (fatsize - support) - - slice_index_t affected_slice_b(sample_index_t t0, unsigned int support) const { - return floor_div(t0 - fftsize + (fat_part(fftsize) - support), filet_part(fftsize)) + 1; - } - - // Index of first slice not affected by sample at t1 - slice_index_t affected_slice_e(sample_index_t t1, unsigned int support) const { - return floor_div(t1 - 1 - (fat_part(fftsize) - support), filet_part(fftsize)) + 1; - } - - bool ok; - bool synthesis; - - unsigned int fftsize_log2; // log2(fftsize) - unsigned int fftsize; // The size of the main FFT, a power of two. - unsigned int fat_size; - unsigned int filet_size; - - // The width of the widest filter in the time domain, in - // octave subsamples - unsigned int oct_support; - - double inv_fftsize_double; // 1.0 / fftsize - T inv_fftsize_t; // 1.0f / fftsize (if using floats) - - unsigned int sftsize_max; // The size of the largest band FFT, a power of two - downsampling_params dsparams; - - // Fourier transform object for transforming a full slice -#if GABORATOR_USE_REAL_FFT - rfft *rft; -#else - fft *ft; -#endif - }; - - // Calculate per-plan, per-band coefficients for plan "pno", - // a synthesis plan if "syn" is true, otherwise an analysis plan. - - void make_band_plans(int pno, bool syn) { - std::vector::plan>> &plans - (syn ? syn_plans : anl_plans); - plan &plan(*plans[pno].get()); - - for (int zno = 0; zno < (int)zones.size(); zno++) { - zone *z = zones[zno].get(); - - make_band_plans_2(z->bandparams, pno, syn, false); - make_band_plans_2(z->mock_bandparams, pno, syn, true); - - if (plan.synthesis) { - // Accumulate window power for calculating dual - std::vector power(plan.fftsize); - // Real bands - for (unsigned int i = 0; i < z->bandparams.size(); i++) { - band_params *bp = z->bandparams[i].get(); - band_plan *bpl = syn ? &bp->syn_plans[pno] : &bp->anl_plans[pno]; - accumulate_power(plan, bp, bpl, power.data()); - } - // Mock bands - for (unsigned int i = 0; i < z->mock_bandparams.size(); i++) { - band_params *bp = z->mock_bandparams[i].get(); - band_plan *bpl = syn ? &bp->syn_plans[pno] : &bp->anl_plans[pno]; - accumulate_power(plan, bp, bpl, power.data()); - } - - // Calculate duals - for (unsigned int obno = 0; obno < z->bandparams.size(); obno++) { - band_params *bp = z->bandparams[obno].get(); - band_plan *bpl = syn ? &bp->syn_plans[pno] : &bp->anl_plans[pno]; - for (unsigned int i = 0; i < bpl->sftsize; i++) { - // ii = large-FFT bin number - int ii = i + bpl->fq_offset_int; - bpl->dual_kernel[i] /= power[ii & (plan.fftsize - 1)]; - } - // The analysis kernels are no longer needed - bpl->kernel = std::vector(); - bpl->shift_kernel = pod_vector(); - } - z->mock_bandparams.clear(); - } - } - } - - void make_band_plans_2(std::vector>> &bv, int pno, - bool syn, bool mock) - { - std::vector::plan>> &plans - (syn ? syn_plans : anl_plans); - plan &plan(*plans[pno].get()); - - for (unsigned int obno = 0; obno < bv.size(); obno++) { - band_params *bp = bv[obno].get(); - std::vector> *bplv = syn ? &bp->syn_plans : &bp->anl_plans; - // XXX redundant resizes - bplv->resize(plans.size()); - band_plan *bpl = &(*bplv)[pno]; - - // Note that bp->step_log2 cannot be negative, meaning - // that the bands can only be subsampled, not oversampled. - unsigned int sftsize = plan.fftsize >> bp->step_log2; - - // PFFFT has a minimum size - sftsize = std::max(sftsize, (unsigned int)GABORATOR_MIN_FFT_SIZE); - - bpl->sftsize = sftsize; - bpl->sftsize_log2 = whichp2(bpl->sftsize); - - if (! mock) { - plan.sftsize_max = std::max(plan.sftsize_max, bpl->sftsize); - bpl->sft = pool, int>::shared.get(bpl->sftsize); - } - - bpl->kernel.resize(bpl->sftsize); - bpl->shift_kernel.resize(bpl->sftsize); - if (plan.synthesis) { - bpl->dual_kernel.resize(bpl->sftsize); - bpl->shift_kernel_conj.resize(bpl->sftsize); - } - - if (bp->dc) - bpl->center = 0; - else - bpl->center = bp->ff * plan.fftsize; - bpl->icenter = (int)rint(bpl->center); - bpl->fq_offset_int = bpl->icenter - (bpl->sftsize >> 1); - - // Calculate frequency-domain window kernel, possibly with - // wrap-around - for (unsigned int i = 0; i < bpl->sftsize; i++) - bpl->kernel[i] = 0; - // i loops over the kernel, with i=0 at the center. - // The range is twice the support on each side so that - // any excess space in the kernel due to rounding up - // the size to a power of two is filled in with actual - // Gaussian values rather than zeros. - int fq_support = (int)ceil(bp->ff_support * plan.fftsize); - for (int i = - 2 * fq_support; i < 2 * fq_support; i++) { - // ii = large-FFT band number of this kernel sample - int ii = i + bpl->fq_offset_int + (int)bpl->sftsize / 2; - // this_ff = fractional frequency of this kernel sample - double this_ff = ii * plan.inv_fftsize_double; - // ki = kernel index - int ki = ii - bpl->fq_offset_int; - // When sftsize == fftsize, the support of the kernel can - // exceed sftsize, and in this case, it should be allowed - // to wrap so that it remains smooth. When sftsize < fftsize, - // sftsize is large enough for the support and no wrapping - // is needed or wanted. - if (bpl->kernel.size() == plan.fftsize && !mock) { - bpl->kernel[ki & (plan.fftsize - 1)] += - eval_kernel(¶ms, bp, this_ff); - if (plan.synthesis) - bpl->dual_kernel[ki & (plan.fftsize - 1)] += - eval_dual_kernel(¶ms, bp, this_ff); - } else { - if (ki >= 0 && ki < (int)bpl->kernel.size()) { - bpl->kernel[ki] += eval_kernel(¶ms, bp, this_ff); - if (plan.synthesis) - bpl->dual_kernel[ki] = eval_dual_kernel(¶ms, bp, this_ff); - } - } - } - } - - // Calculate complex exponentials for non-integer center - // frequency adjustment and phase convention adjustment - for (unsigned int obno = 0; obno < bv.size(); obno++) { - band_params *bp = bv[obno].get(); - band_plan *bpl = syn ? &bp->syn_plans[pno] : &bp->anl_plans[pno]; - for (unsigned int i = 0; i < bpl->sftsize; i++) { - double center = - (params.phase == coef_phase::global) ? bpl->center : 0; - double arg = tau * ((double)i / bpl->sftsize) * -(center - bpl->icenter); - C t(cos(arg), sin(arg)); - // Apply ifftshift of spectrum in time domain - bpl->shift_kernel[i] = (i & 1) ? -t : t; - if (plan.synthesis) - // Conjugate kernel does not have ifftshift - bpl->shift_kernel_conj[i] = conj(t); - } - } - } - - // Add the power of the kernel in "*bp" to "power" - void - accumulate_power(plan &plan, band_params *bp, band_plan *bpl, T *power) { - for (unsigned int i = 0; i < bpl->sftsize; i++) { - // ii = large-FFT bin number - unsigned int ii = i + bpl->fq_offset_int; - ii &= plan.fftsize - 1; - assert(ii >= 0 && ii < plan.fftsize); - T p = bpl->kernel[i] * bpl->dual_kernel[i]; - power[ii] += p; - if (params.lowpass_version == 2 || ! bp->dc) { - unsigned int ni = -ii; - ni &= plan.fftsize - 1; - power[ni] += p; - } - } - } - - // Create coefficient metadata based on a slice length - coefs_meta *make_meta(int slice_len) const { - coefs_meta *cmeta = new coefs_meta; - cmeta->n_octaves = n_octaves; - cmeta->n_bands_total = n_bands_total; - cmeta->bands_per_octave = params.bands_per_octave; - cmeta->slice_len = slice_len; - cmeta->zones.resize(zones.size()); - for (unsigned int zi = 0; zi < zones.size(); zi++) { - zone *z = zones[zi].get(); - typename zone_coefs_meta::band_vector bv(z->bandparams.size()); - for (unsigned int i = 0; i < z->bandparams.size(); i++) { - unsigned int step_log2 = z->bandparams[i]->step_log2; - bv[i].step_log2 = step_log2; - bv[i].slice_len_log2 = whichp2(slice_len) - step_log2; - bv[i].slice_len = 1 << bv[i].slice_len_log2; - } - cmeta->zones[zi].init(bv); - } - cmeta->octaves.resize(octaves.size()); - int tbno = 0; - for (unsigned int i = 0; i < n_octaves; i++) { - cmeta->octaves[i].z = &cmeta->zones[this->octaves[i].z->zno]; - cmeta->octaves[i].n_bands_above = tbno; - tbno += cmeta->octaves[i].z->bands.size(); - } - return cmeta; - } - - std::vector> octaves; // Per-octave parameters - std::vector>> zones; - unsigned int max_step_log2; - - std::vector> anl_plans; - std::vector> syn_plans; - - unsigned int n_octaves; - unsigned int n_bands_total; // Total number of frequency bands, including DC - - double top_band_log2ff; // log2 of fractional frequency of the highest-frequency band - int ffref_gbno; // Band number of the reference frequency - - // Width of the downsampling filter passband in terms of the - // downsampled sample rate (between 0.25 and 0.5) - double ds_passband; - double ds_ff; // Downsampling filter -6 dB transition frequency - double ds_ff_sd; // Downsampling filter standard deviation - double ds_time_support; // Downsampling filter time-domain kernel support, each side - - unsigned int fftsize_max; // Largest FFT size of any plan - unsigned int sftsize_max; // Largest SFT size of any plan - - ref cmeta_any; -}; - - -// Iterate over the slices of a row (band) having slice length -// 2^sh that contain coefficients with indices ranging from i0 -// (inclusive) to i1 (exclusive), and call the function f for -// each such slice (full or partial), with the arguments -// -// sli - slice index -// bvi - index of first coefficient to process within the slice -// len - number of coefficients to process within the slice - -template -void foreach_slice(unsigned int sh, coef_index_t i0, coef_index_t i1, F f) { - // Note that this can be called with i0 > i1 and needs to handle - // that case gracefully. - // Band size (power of two) - int bsize = 1 << sh; - // Adjust for t=0 being outside the filet - int fatsize = bsize >> 1; - i0 -= fatsize; - i1 -= fatsize; - coef_index_t i = i0; - while (i < i1) { - // Slice index - slice_index_t sli = i >> sh; - // Band vector index - unsigned int bvi = i & (bsize - 1); - unsigned int len = bsize - bvi; - unsigned int remain = (unsigned int)(i1 - i); - if (remain < len) - len = remain; - f(sli, bvi, len); - i += len; - } -} - -// As foreach_slice, but call the "process_existing_slice" method of -// the given "dest" object for each full or partial slice of -// coefficients, and/or the "process_missing_slice" method for each -// nonexistent slice. -// -// Template parameters: -// T is the spectrogram value type -// D is the dest object type -// C is the coefficient type - -template > -struct row_foreach_slice { - typedef C value_type; - row_foreach_slice(const coefs &msc, - int oct_, unsigned int obno_): - oct(oct_), obno(obno_), sc(msc.octaves[oct]), - sh(sc.meta->bands[obno].slice_len_log2) - { - assert(oct < (int)msc.octaves.size()); - } -public: - void operator()(coef_index_t i0, coef_index_t i1, D &dest) const { - foreach_slice(sh, i0, i1, - [this, &dest](slice_index_t sli, unsigned int bvi, unsigned int len) { - oct_coefs *c = get_existing_coefs(sc, sli); - if (c) { - dest.process_existing_slice(c->bands[obno] + bvi, len); - } else { - dest.process_missing_slice(len); - } - }); - } - int oct; - unsigned int obno; - const sliced_coefs ≻ - unsigned int sh; -}; - -// Helper class for row_source - -template -struct writer_dest { - writer_dest(OI output_): output(output_) { } - void process_existing_slice(C *bv, size_t len) { - // Can't use std::copy here because it takes the output - // iterator by value, and using the return value does not - // work, either. - for (size_t i = 0; i < len; i++) - *output++ = bv[i]; - } - void process_missing_slice(size_t len) { - for (size_t i = 0; i < len; i++) - *output++ = C(); - } - OI output; -}; - -// Retrieve a sequence of coefficients from a row (band) in the -// spectrogram, with indices ranging from i0 to i1. The indices can -// be negative, and can extend outside the available data, in which -// case zero is returned. The coefficients are written through the -// output iterator "output". -// Template arguments: -// T is the spectrogram value type -// OI is the output iterator type -// C is the coefficient value type -// Note defaults for template arguments defined in forward declaration. - -template -struct row_source { - row_source(const coefs &msc_, - int oct_, unsigned int obno_): - slicer(msc_, oct_, obno_) - { } - OI operator()(coef_index_t i0, coef_index_t i1, OI output) const { - writer_dest dest(output); - slicer(i0, i1, dest); - return dest.output; - } - row_foreach_slice, C> slicer; -}; - -// The opposite of row_source: store a sequence of coefficients into -// a row (band) in the spectrogram. This duplicates quite a lot of -// the row_source code above (without comments); the main part that's -// different is marked by the comments "Begin payload" and "End -// payload". Other differences: iterator is called II rather than OI, -// and the coefs are not const. - -template -struct row_dest { - typedef C value_type; - row_dest(coefs &msc, - int oct_, unsigned int obno_): - oct(oct_), obno(obno_), sc(msc.octaves[oct]), - sh(sc.meta->bands[obno].slice_len_log2) - { - assert(oct < (int)msc.octaves.size()); - } -public: - II operator()(coef_index_t i0, coef_index_t i1, II input) const { - assert(i0 <= i1); - int bsize = 1 << sh; - int fatsize = bsize >> 1; - i0 -= fatsize; - i1 -= fatsize; - coef_index_t i = i0; - while (i < i1) { - slice_index_t sli = i >> sh; - unsigned int bvi = i & (bsize - 1); - unsigned int len = bsize - bvi; - unsigned int remain = (unsigned int)(i1 - i); - if (remain < len) - len = remain; - int bvie = bvi + len; - // Begin payload - oct_coefs *c = &get_or_create_coefs(sc, sli); - C *bv = c->bands[obno]; - for (int j = bvi; j < bvie; j++) - bv[j] = *input++; - i += len; - // End payload - } - return input; - } - int oct; - unsigned int obno; - sliced_coefs ≻ - unsigned int sh; -}; - -// One more set of duplicated code, now for adding to coefficients - -template -struct row_add_dest { - typedef C value_type; - row_add_dest(coefs &msc, - int oct_, unsigned int obno_): - oct(oct_), obno(obno_), sc(msc.octaves[oct]), - sh(sc.meta->bands[obno].slice_len_log2) - { - assert(oct < (int)msc.octaves.size()); - } -public: - II operator()(coef_index_t i0, coef_index_t i1, II input) const { - assert(i0 <= i1); - int bsize = 1 << sh; - int fatsize = bsize >> 1; - i0 -= fatsize; - i1 -= fatsize; - coef_index_t i = i0; - while (i < i1) { - slice_index_t sli = i >> sh; - unsigned int bvi = i & (bsize - 1); - unsigned int len = bsize - bvi; - unsigned int remain = (unsigned int)(i1 - i); - if (remain < len) - len = remain; - int bvie = bvi + len; - // Begin payload - oct_coefs *c = &get_or_create_coefs(sc, sli); - C *bv = c->bands[obno]; - for (int j = bvi; j < bvie; j++) - bv[j] += *input++; - i += len; - // End payload - } - return input; - } - int oct; - unsigned int obno; - sliced_coefs ≻ - unsigned int sh; -}; - -// Helper for process() below. Here, the function f() operates on an -// array of consecutive coefficient samples rather than a single -// sample. - -template -void apply_to_slice(bool create, - F f, - int b0, // = INT_MIN - int b1, // = INT_MAX - sample_index_t st0, // = INT64_MIN - sample_index_t st1, // = INT64_MAX - coefs& coefs0, - coefs&... coefsi) -{ - b0 = std::max(b0, 0); - b1 = std::min(b1, (int)coefs0.meta->n_bands_total); - for (int band = b0; band < b1; band++) { - int oct; - unsigned int obno; - bool valid = gaborator::bno_split(*coefs0.meta, band, oct, obno, true); - assert(valid); - - int exp = coefs0.meta->octaves[oct].z->bands[obno].step_log2 + oct; - int time_step = 1 << exp; - - coef_index_t ci0 = (st0 + time_step - 1) >> exp; - coef_index_t ci1 = ((st1 - 1) >> exp) + 1; - if (! create) { - // Restrict to existing coefficient index range - coef_index_t cib0, cib1; - get_band_coef_bounds(coefs0, oct, obno, cib0, cib1); - ci0 = std::max(ci0, cib0); - ci1 = std::min(ci1, cib1); - } - unsigned int sh = coefs0.meta->octaves[oct].z->bands[obno].slice_len_log2; - sample_index_t st = shift_left(ci0, exp); - foreach_slice(sh, ci0, ci1, - [&](slice_index_t sli, unsigned int bvi, - unsigned int len) - { - oct_coefs *c = create ? - &get_or_create_coefs(coefs0.octaves[oct], sli) : - get_existing_coefs(coefs0.octaves[oct], sli); - if (c) { - // p0 points to coefficient from the first set - C0 *p0 = c->bands[obno] + bvi; - f(band, st, time_step, len, p0, - get_or_create_coefs(coefsi.octaves[oct], sli).bands[obno] + bvi...); - } - st += len * time_step; - }); - } -} - -// Common implementation of process() and fill() - -template -void apply_common(bool create, - F f, - int b0, // = INT_MIN - int b1, // = INT_MAX - sample_index_t st0, // = INT64_MIN - sample_index_t st1, // = INT64_MAX - coefs &coefs0, - coefs&... coefsi) -{ - apply_to_slice(create, - [&](int bno, int64_t st, int time_step, - unsigned len, C0 *p0, CI *...pi) - { - for (unsigned int i = 0; i < len; i++) { - f(bno, st, *p0++, *pi++...); - st += time_step; - } - }, b0, b1, st0, st1, coefs0, coefsi...); -} - -// Iterate over one or more coefficient sets in parallel and apply the -// function f, passing it a coefficient from each set as an argument. -// -// The application can be optionally limited to coefficients within -// the band range b0 to b1 and/or the sample time range st0 to st1. -// -// The first coefficient set ("channel 0") is treated specially; it -// determines which coefficients are iterated over (optionally further -// restricted by b0/b1/st0/st1). That is, the iteration is over the -// coefficients that already exist in channel 0, and no new -// coefficients will be allocated in channel 0. In the other -// channels, new coefficients will be created on demand when they are -// missing from that channel but present in channel 0. -// -// The coefficients may be of a different data type in each set. -// -// The arguments to f() are: -// -// int bno Band number -// int64_t st Sample time -// C &p0 Coefficient from first set -// C... &pi Coefficient from subsequent set - -template -void process(F f, - int b0, // = INT_MIN - int b1, // = INT_MAX - sample_index_t t0, // = INT64_MIN - sample_index_t t1, // = INT64_MAX - coefs &coefs0, - coefs&... coefsi) -{ - apply_common(false, f, b0, b1, t0, t1, coefs0, coefsi...); -} - -template -void fill(F f, - int b0, // = INT_MIN - int b1, // = INT_MAX - sample_index_t t0, // = INT64_MIN - sample_index_t t1, // = INT64_MAX - coefs &coefs0, - coefs&... coefsi) -{ - apply_common(true, f, b0, b1, t0, t1, coefs0, coefsi...); -} - -// Apply the function f to each existing coefficient in the -// coefficient set msc within the time range st0 to st1. The -// initial analyzer argument is ignored. -// -// The arguments to f() are: -// -// C &c Coefficient -// int b Band number -// int64_t t Time in samples -// -// This is for backwards compatibility; process() is now preferred. - -template -void apply(const analyzer &, coefs &msc, F f, - sample_index_t st0 = INT64_MIN, - sample_index_t st1 = INT64_MAX) -{ - process([&](int b, int64_t t, complex& c) { - f(c, b, t); - }, INT_MIN, INT_MAX, st0, st1, msc); -} - -template -void forget_before(const analyzer &, coefs &msc, - sample_index_t limit, bool clean_cut = false) -{ - typedef complex C; - unsigned int n_oct = (unsigned int) msc.octaves.size(); - for (unsigned int oct = 0; oct < n_oct; oct++) { - sliced_coefs &sc = msc.octaves[oct]; - // Convert limit from samples to slices, rounding down. - // This assumes all bands in the octave have the same - // time range, as they must. - // First convert samples to coefficients, rounding down - int obno = 0; // Any band would do, and band 0 always exists - zone_coefs_meta *zmeta = msc.meta->octaves[oct].z; - coef_index_t ci = limit >> band_scale_exp(*zmeta, oct, obno); - // Then convert coefficients to slices, accounting for - // fat and rounding down - unsigned int slice_len = zmeta->bands[obno].slice_len; - unsigned int slice_len_log2 = zmeta->bands[obno].slice_len_log2; - int fat = slice_len >> 1; - slice_index_t sli = (ci - fat) >> slice_len_log2; - sc.slices.erase_before(sli); - if (clean_cut) { - // Partially erase slice at boundary, if any - const ref> *t = sc.slices.get(sli); - if (! t) - continue; - if (! *t) - continue; - const oct_coefs &c = **t; - unsigned int n_bands = (unsigned int)c.bands.size(); - for (unsigned int obno = 0; obno < n_bands; obno++) { - C *band = c.bands[obno]; - unsigned int len = sc.meta->bands[obno].slice_len; - sample_index_t st = sample_time(*sc.meta, sli, 0, oct, obno); - int time_step = 1 << band_scale_exp(*sc.meta, oct, obno); - for (unsigned int i = 0; i < len; i++) { - if (st < limit) - band[i] = 0; - else - break; - st += time_step; - } - } - } - } -} - - -} // namespace - -#endif diff --git a/lib/gaborator/gaborator/gaussian.h b/lib/gaborator/gaborator/gaussian.h deleted file mode 100644 index 1ee76b3..0000000 --- a/lib/gaborator/gaborator/gaussian.h +++ /dev/null @@ -1,123 +0,0 @@ -// -// The Gaussian and related functions -// -// Copyright (C) 2015-2021 Andreas Gustafsson. This file is part of -// the Gaborator library source distribution. See the file LICENSE at -// the top level of the distribution for license information. -// - -#ifndef _GABORATOR_GAUSSIAN_H -#define _GABORATOR_GAUSSIAN_H - -#include - -namespace gaborator { - -// A rough approximation of erfc_inv(), the inverse complementary -// error function. This is good for arguments in the range 1e-8 to 1, -// to within a few percent. - -inline float erfc_inv(float x) { - return sqrtf(-logf(x)) - 0.3f; -} - -// Gaussian with peak = 1 - -inline double norm_gaussian(double sd, double x) { - return exp(-(x * x) / (2 * sd * sd)); -} - -// Gaussian with integral = 1 - -inline double gaussian(double sd, double x) { - double a = 1.0 / (sd * sqrt(2.0 * M_PI)); - return a * norm_gaussian(sd, x); -} - -// The convolution of a Heaviside step function with a Gaussian of -// standard deviation sd. Goes smoothly from 0 to 1, with the 0.5 -// point at x=0. - -static inline -double gaussian_edge(double sd, double x) { - double erf_arg = x / (sd * M_SQRT2); - if (erf_arg < -7) - return 0; // error < 5e-23 - if (erf_arg > 7) - return 1; // error < 5e-23 - return (erf(erf_arg) + 1) * 0.5; -} - -// Translate the time-domain standard deviation of a Gaussian -// (in samples) into the corresponding frequency-domain standard -// deviation (as a fractional frequency), or vice versa. - -static inline double sd_t2f(double st_sd) { - return 1.0 / (2.0 * M_PI * st_sd); -} - -static inline double sd_f2t(double ff_sd) { - return sd_t2f(ff_sd); -} - -// Given a Gaussian with standard deviation "sd" and a maximum error -// "max_error", calculate the support needed on each side to keep the -// area below the curve within max_error of the exact value. - -static inline -double gaussian_area_support(double sd, double max_error) { - return sd * M_SQRT2 * erfc_inv(max_error); -} - -// Inverse of the above: given a support and maximum error, calculate -// the standard deviation. - -static inline -double gaussian_area_support_inv(double support, double max_error) { - return support / (M_SQRT2 * erfc_inv(max_error)); -} - -// Given a gaussian with standard deviation "sd" and a maximum error -// "max_error", calculate the support needed on each side for the -// value to fall to a factor of "max_error" of the peak. - -static inline -double gaussian_value_support(double sd, double max_error) { - return sd * M_SQRT2 * sqrt(-log(max_error)); -} - -// Inverse of the above: given a support and maximum error, calculate -// the standard deviation. - -static inline -double gaussian_value_support_inv(double support, double max_error) { - return support / (M_SQRT2 * sqrt(-log(max_error))); -} - -// Choose which criterion to use - -#if 1 -static inline -double gaussian_support(double support, double max_error) { - return gaussian_area_support(support, max_error); -}; - -static inline -double gaussian_support_inv(double support, double max_error) { - return gaussian_area_support_inv(support, max_error); -}; -#else -static inline -double gaussian_support(double support, double max_error) { - return gaussian_value_support(support, max_error); -}; - -static inline -double gaussian_support_inv(double support, double max_error) { - return gaussian_value_support_inv(support, max_error); -}; -#endif - -} // namespace - -#endif diff --git a/lib/gaborator/gaborator/pod_vector.h b/lib/gaborator/gaborator/pod_vector.h deleted file mode 100644 index e578864..0000000 --- a/lib/gaborator/gaborator/pod_vector.h +++ /dev/null @@ -1,116 +0,0 @@ -// -// A vector class without default-initialization, for "plain old data" -// -// Copyright (C) 2016-2021 Andreas Gustafsson. This file is part of -// the Gaborator library source distribution. See the file LICENSE at -// the top level of the distribution for license information. -// - -#ifndef _GABORATOR_POD_VECTOR_H -#define _GABORATOR_POD_VECTOR_H - -#include // size_t - -#include // std::swap - -namespace gaborator { - -// A vector for storing Plain Old Data. This is similar to a -// std::vector, except that it does not zero-initialize elements, -// and that it guarantees that data() returns a non-NULL pointer -// even when the vector contains zero elements. - -template -struct pod_vector { - typedef T value_type; - typedef T *iterator; - pod_vector() { - b = e = 0; - } - explicit pod_vector(size_t size_) { - // Allocate raw uninitialized memory - b = static_cast(::operator new(size_ * sizeof(T))); - e = b + size_; - } - ~pod_vector() -#if __cplusplus >= 201103L - noexcept -#endif - { - _free(); - } - void swap(pod_vector &x) { - std::swap(b, x.b); - std::swap(e, x.e); - } - iterator begin() const { return b; } - iterator end() const { return e; } - T *data() { return b; } - const T *data() const { return b; } - T &operator[](size_t i) { return b[i]; } - const T &operator[](size_t i) const { return b[i]; } - size_t size() const { return e - b; } - void resize(size_t new_size) { - if (new_size == size()) - return; - T *n = static_cast(::operator new(new_size * sizeof(T))); - size_t ncopy = std::min(size(), new_size); - std::copy(b, b + ncopy, n); - _free(); - b = n; - e = n + new_size; - } - pod_vector(const pod_vector &a) - - { - b = new T[a.size()]; - e = b + a.size(); - std::copy(a.b, a.e, b); - //if (size()) fprintf(stderr, "pod_vector cc %d\n", (int)size()); - } - void clear() { - _free(); - b = e = 0; - } -#if __cplusplus >= 201103L - pod_vector(pod_vector&& x) noexcept: - b(x.b), e(x.e) - { - x.b = x.e = 0; - //if (size()) fprintf(stderr, "pod_vector mv %d\n", (int)size()); - } -#endif - pod_vector &operator=(const pod_vector &a) { - if (&a == this) - return *this; - _free(); - b = static_cast(::operator new(a.size() * sizeof(T))); - e = b + a.size(); - std::copy(a.b, a.e, b); - //if (size()) fprintf(stderr, "pod_vector = %d\n", (int)size()); - return *this; - } -#if __cplusplus >= 201103L - pod_vector &operator=(pod_vector &&x) noexcept { - if (&x == this) - return *this; - _free(); - b = x.b; - e = x.e; - x.b = x.e = 0; - return *this; - } -#endif -private: - void _free() { - // Free as raw uninitialized memory - ::operator delete(b); - } -private: - T *b; - T *e; -}; - -} // namespace - -#endif diff --git a/lib/gaborator/gaborator/pool.h b/lib/gaborator/gaborator/pool.h deleted file mode 100644 index 47e8322..0000000 --- a/lib/gaborator/gaborator/pool.h +++ /dev/null @@ -1,47 +0,0 @@ -// -// Pool of shared objects -// -// Copyright (C) 2015-2018 Andreas Gustafsson. This file is part of -// the Gaborator library source distribution. See the file LICENSE at -// the top level of the distribution for license information. -// - -#ifndef _GABORATOR_POOL_H -#define _GABORATOR_POOL_H - -#include - -namespace gaborator { - -// The "pool" class is for sharing FFT objects so that we don't -// create multiple FFTs of the same size. It could also be used -// to share objects of some other class T where we don't want to -// create multiple Ts with the same K. - -template -struct pool { - typedef std::map m_t; - ~pool() { - for (typename m_t::iterator it = m.begin(); it != m.end(); it++) { - delete (*it).second; - } - } - T *get(const K &k) { - std::pair r = m.insert(std::make_pair(k, (T *)0)); - if (r.second) { - // New element was inserted - assert((*(r.first)).second == 0); - (*r.first).second = new T(k); - } - return (*r.first).second; - } - m_t m; - static pool shared; -}; - -template -pool pool::shared; - -} // namespace - -#endif diff --git a/lib/gaborator/gaborator/ref.h b/lib/gaborator/gaborator/ref.h deleted file mode 100644 index fb4562f..0000000 --- a/lib/gaborator/gaborator/ref.h +++ /dev/null @@ -1,80 +0,0 @@ -// -// Intrusive reference counting smart pointer -// -// Copyright (C) 2016-2018 Andreas Gustafsson. This file is part of -// the Gaborator library source distribution. See the file LICENSE at -// the top level of the distribution for license information. -// - -#ifndef _GABORATOR_REF_H -#define _GABORATOR_REF_H - -namespace gaborator { - -template struct ref; - -struct refcounted { - refcounted() { refcount = 0; } - unsigned int refcount; -}; - -// Template functions for manual reference counting, without using the -// ref<> class. It would be tempting to make these methods of struct -// refcounted, but that won't work because it would lose the full -// object type and invoke operator delete on the base class. - -template -void incref(T &r) { - r.refcount++; -} - -template -void decref(T &r) { - r.refcount--; - if (r.refcount == 0) - delete &r; -} - -template -struct ref { - ref(): p(0) { } - ref(T *p_): p(p_) { - _incref(); - } - ref(const ref &o): p(o.p) { - _incref(); - } - ref &operator=(const ref &o) { reset(o.p); return *this; } - ~ref() { reset(); } - void reset() { - _decref(); - p = 0; - } - void reset(T *n) { - if (n == p) - return; - _decref(); - p = n; - _incref(); - } - T *get() const { return p; } - T *operator->() const { return p; } - T &operator*() const { return *p; } - operator bool() const { return p; } -private: - void _incref() { - if (! p) - return; - incref(*p); - } - void _decref() { - if (! p) - return; - decref(*p); - } - T *p; -}; - -} // namespace - -#endif diff --git a/lib/gaborator/gaborator/render.h b/lib/gaborator/gaborator/render.h deleted file mode 100644 index 2b8df98..0000000 --- a/lib/gaborator/gaborator/render.h +++ /dev/null @@ -1,506 +0,0 @@ -// -// Rendering of spectrogram images -// -// Copyright (C) 2015-2021 Andreas Gustafsson. This file is part of -// the Gaborator library source distribution. See the file LICENSE at -// the top level of the distribution for license information. -// - -#ifndef _GABORATOR_RENDER_H -#define _GABORATOR_RENDER_H - -#include "gaborator/gaborator.h" -#include "gaborator/resample2.h" - -namespace gaborator { - - -// Convert a floating-point linear brightness value in the range 0..1 -// into an 8-bit pixel value, with clamping and (rough) gamma -// correction. This nominally uses the sRGB gamma curve, but the -// current implementation cheats and uses a gamma of 2 because it can -// be calculated quickly using a square root. - -template -unsigned int float2pixel_8bit(T val) { - // Clamp before gamma correction so we don't take the square root - // of a negative number; those can arise from bicubic - // interpolation. While we're at it, let's also skip the gamma - // correction for small numbers that will round to zero anyway, - // and especially denormals which could rigger GCC bug target/83240. - static const T almost_zero = 1.0 / 65536; - if (val < almost_zero) - val = 0; - if (val > 1) - val = 1; - return (unsigned int)(sqrtf(val) * 255.0f); -} - -////////////////////////////////////////////////////////////////////////////// - -// Power-of-two rendering - - -// Magnitude - -template -struct complex_abs_fob { - T operator()(const complex &c) { - return complex_abs(c); - } - typedef T return_type; -}; - - -// T -> f() -> OI::value_type - -template -struct transform_output_iterator: public std::iterator { - typedef T value_type; - transform_output_iterator(F f_, OI output_): f(f_), output(output_) { } - transform_output_iterator& operator=(T v) { - *output++ = f(v); - return *this; - } - transform_output_iterator& operator*() { return *this; } - transform_output_iterator& operator++() { return *this; } - transform_output_iterator& operator++(int) { return *this; } - F f; - OI output; -}; - -// A source object for resample2() that provides the -// values of a row of spectrogram coordinates transformed -// by function normf (typically an absolute value function). -// -// T is the analyzer signal type -// C is the coefficient type -// OI is the output iterator type - -template -struct abs_row_source { - typedef transform_output_iterator abs_writer_t; - abs_row_source(const coefs &coefs_, - int oct_, unsigned int obno_, - NORMF normf_): - rs(coefs_, oct_, obno_), - normf(normf_) - { } - OI operator()(sample_index_t i0, sample_index_t i1, OI output) const { - abs_writer_t abswriter(normf, output); - abs_writer_t abswriter_end = rs(i0, i1, abswriter); - return abswriter_end.output; - } - row_source rs; - NORMF normf; -}; - -// Helper class for abs_row_source specialization below - -template -struct abs_writer_dest { - abs_writer_dest(OI output_): output(output_) { } - void process_existing_slice(C *bv, size_t len) { - complex_magnitude(bv, output, len); - output += len; - } - void process_missing_slice(size_t len) { - for (size_t i = 0; i < len; i++) - *output++ = 0; - } - OI output; -}; - -// Partial specialization of class abs_row_source for NORMF = complex_abs_fob, -// for vectorization. - -template -struct abs_row_source> { - // Note unused last arg - abs_row_source(const coefs &coefs_, - int oct_, unsigned int obno_, - complex_abs_fob): - slicer(coefs_, oct_, obno_) - { } - OI operator()(coef_index_t i0, coef_index_t i1, OI output) const { - abs_writer_dest dest(output); - slicer(i0, i1, dest); - return dest.output; - } - row_foreach_slice, C> slicer; -}; - -// Render a single line (single frequency band), with scaling in the -// horizontal (time) dimension, and filtering to avoid aliasing when -// minifying. - -template -OI render_line(const analyzer &anl, - const coefs &msc, - int gbno, - affine_transform xf, - sample_index_t i0, sample_index_t i1, - OI output, - NORMF normf) -{ - typedef typename NORMF::return_type RST; - - int oct; - unsigned int obno; // Band number within octave - bool clip = ! bno_split(*msc.meta, gbno, oct, obno, false); - if (clip) { - for (sample_index_t i = i0; i < i1; i++) - *output++ = (T)0; - return output; - } - abs_row_source - abs_rowsource(msc, oct, obno, normf); - - // Scale by the downsampling factor of the band - int scale_exp = band_scale_exp(*msc.octaves[oct].meta, oct, obno); - RESAMPLER x_resampler(zoom_p2(xf, -scale_exp)); - output = x_resampler(abs_rowsource, i0, i1, output); - return output; -} - -// Render a two-dimensional image with scaling in the horizontal -// direction only. In the vertical direction, there is always a -// one-to-one correspondence between bands and pixels. yi0 and yi1 -// already have the yorigin applied, so there is no yorigin argument. - -template -OI render_noyscale(const analyzer &anl, - const coefs &msc, - affine_transform x_xf, - int64_t xi0, int64_t xi1, - int64_t yi0, int64_t yi1, - OI output, - int64_t output_stride, - NORMF normf) -{ - assert(xi1 >= xi0); - int w = (int)(xi1 - xi0); - assert(w >= 0); - int gbno0 = (int)yi0; - int gbno1 = (int)yi1; - for (int gbno = gbno0; gbno < gbno1; gbno++) { - render_line - (anl, msc, gbno, x_xf, - xi0, xi1, - output, normf); - output += output_stride; - } - return output; -} - -// Source data from a column of a row-major two-dimensional array. -// data points to the beginning of a row-major array with an x -// range of x0..x1 and an y range from y0..y1, and operator() -// returns data from column x (where x is within the range x0..x1). - -template -struct transverse_source { - transverse_source(T *data_, - int64_t x0_, int64_t x1_, int64_t y0_, int64_t y1_, - int64_t x_): - data(data_), - x0(x0_), x1(x1_), y0(y0_), y1(y1_), - x(x_), - stride(x1 - x0) - { } - OI operator()(int64_t i0, int64_t i1, OI out) const { - assert(x >= x0); - assert(x <= x1); - assert(i1 >= i0); - assert(i0 >= y0); - assert(i1 <= y1); - T *p = data + (x - x0) + (i0 - y0) * stride; - while (i0 != i1) { - *out++ = *p; - p += stride; - ++i0; - } - return out; - } - T *data; - int64_t x0, x1, y0, y1, x; - size_t stride; -}; - -template -struct stride_iterator: public std::iterator - ::value_type> -{ - typedef typename std::iterator_traits::value_type T; - stride_iterator(I it_, size_t stride_): it(it_), stride(stride_) { } - T& operator*() { return *it; } - stride_iterator& operator++() { - it += stride; - return *this; - } - stride_iterator operator++(int) { - stride_iterator old = *this; - it += stride; - return old; - } - I it; - size_t stride; -}; - -struct updated_nop { - void operator()(int64_t x0, int64_t x1, int64_t y0, int64_t y1) { } -}; - -// Render a two-dimensional image with scaling in both the horizontal -// (time) and vertical (frequency) directions. The output may be -// written through "output" out of order, so "output" must be a random -// access iterator. - -// Note the default template argument for NORMF. This is needed -// because the compiler won't deduce the type of NORMF from the -// default function argument "NORMF normf = complex_abs_fob()" -// when the normf argument is omitted; it is considered a "non-deduced -// context", being "a template parameter used in the parameter type of -// a function parameter that has a default argument that is being used -// in the call for which argument deduction is being done". -// Unfortuantely, this work-around of providing a default template -// argument requires C++11. - -// This supports incremental rendering where only the parts of the -// output image affected by analyzing a given time range of samples -// are updated, the time range being from inc_i0 (inclusive) to inc_i1 -// (exclusive). The updated parts of the image consist of zero or -// more non-overlapping rectangles; to find out what those are, pass a -// function "update" which will be called will each rectangle in turn -// after it has been updated. - -// For non-incremental rendering, pass inc_i0 = INT64_MIN and inc_i1 = -// INT64_MAX. - -// OI is the output iterator type -// T is the coefficient type - -template , - class RESAMPLER = lanczos2_pow2_resampler, - class UPDATEDF = updated_nop> -void render_incremental( - const analyzer &anl, - const coefs &msc, - affine_transform x_xf, affine_transform y_xf, - int64_t xi0, int64_t xi1, - int64_t yi0, int64_t yi1, - int64_t inc_i0, int64_t inc_i1, - OI output, - int64_t output_stride, - NORMF normf = complex_abs_fob(), - UPDATEDF updated = updated_nop()) -{ - assert(xi1 >= xi0); - assert(yi1 >= yi0); - assert(inc_i1 >= inc_i0); - - // The data type to reasmple - typedef typename NORMF::return_type RST; - - // Vertical resampler - RESAMPLER y_resampler(y_xf); - - // Find the image bounds in the spectrogram coordinate system, - // including the interpolation margin. The Y bounds are in - // bands and are used both to determine what to render into the - // temporary image and for short-circuiting. The X bounds are in - // coefficient samples, and are only used for short-circuiting. - // The X bounds will be calculated later if we need them. - int64_t ysi0, ysi1; - y_resampler.support(yi0, yi1, ysi0, ysi1); - - // Calculate adjusted image X bounds based on the updated signal - // range for incremental rendering, and return an estimate of the - // numbers of pixels we can avoid rendering at this Y coordinate - // by using the adjusted X bounds. - - auto savings = [&](int64_t y, int64_t &adj_x0, int64_t &adj_x1) { - // Find the highest-index / lowest-frequency band used - // as a resampler input for output pixel y; it will have - // the widest analysis support in the x direction. - // Note that we pass y twice, and ignore the ysi0 result. - int64_t ysi0, ysi1; - y_resampler.support(y, y, ysi0, ysi1); - int64_t band = ysi1; - // Clamp the band to the valid range - band = std::max(band, (int64_t)anl.bandpass_bands_begin()); - band = std::min(band, (int64_t)anl.bandpass_bands_end() - 1); - - // Find the analysis support in the time (x) dimension, - // in signal samples - double support = anl.analysis_support(band); - // Convert from signal samples to coefficient samples - int scale_exp = anl.band_scale_exp((int)band); - - // Extend the updated coefficient range by the analysis - // support, and map it back to pixel space to find the - // affected pixel range, taking the resampler support - // into account. - RESAMPLER x_resampler(zoom_p2(x_xf, -scale_exp)); - int64_t ceil_support = ceil(support); - - // inv_support() calculates both sides of the support at once, - // but in the one-sided case, passing INT64_MIN/MAX may cause - // overflow and undefined behavior. Therefore, we pass a - // dummy value of zero instead, and make sure not to use the - // corresponding output value. This may cause the two inputs - // to inv_support() to be out of order, so it needs to accept - // that. - x_resampler.inv_support( - inc_i0 == INT64_MIN ? (int64_t)0 : (inc_i0 - ceil_support) >> scale_exp, - inc_i1 == INT64_MAX ? (int64_t)0 : (inc_i1 + ceil_support) >> scale_exp, - adj_x0, - adj_x1); - - if (inc_i0 == INT64_MIN) { - adj_x0 = xi0; - } else { - adj_x0 = std::max(xi0, adj_x0); - // Don't let width go negative - adj_x0 = std::min(adj_x0, xi1); - } - - if (inc_i1 == INT64_MAX) { - adj_x1 = xi1; - } else { - adj_x1 = std::min(xi1, adj_x1); - // Don't let width go negative - adj_x1 = std::max(adj_x1, adj_x0); - } - - assert(adj_x0 <= adj_x1); - - return (adj_x0 - xi0) + (xi1 - adj_x1); - }; - - if (!(inc_i0 == INT64_MIN && inc_i1 == INT64_MAX)) { - // Incremental rendering has been requested - int64_t adj_x0_top, adj_x1_top, adj_x0_bot, adj_x1_bot; - // See how much rendering we can save per line at the bottom, - // and calcualate adjusted bounds - int64_t bot_savings = savings(ysi1, adj_x0_bot, adj_x1_bot); - // See how much rendering we can save per line at the top - int64_t top_savings = savings(ysi0, adj_x0_top, adj_x1_top); - // Adjust bounds and output pointer to realize the bottom - // savings - if (adj_x0_bot == adj_x1_bot) - return; - output += adj_x0_bot - xi0; - xi0 = adj_x0_bot; - xi1 = adj_x1_bot; - - // If the savings at the top are significantly greater than - // at the bottom, it pays to subdivde the area to render, - // so that the top part can benefit from the greater savings - // there. - if (((top_savings - bot_savings) * (yi1 - yi0)) > 1000) { - // Subdivide vertically - int64_t ysplit = (yi1 + yi0) >> 1; - size_t output_offset = (ysplit - yi0) * output_stride; - render_incremental - (anl, msc, x_xf, y_xf, - xi0, xi1, - yi0, ysplit, - inc_i0, inc_i1, - output, output_stride, normf, updated); - render_incremental - (anl, msc, x_xf, y_xf, - xi0, xi1, - ysplit, yi1, - inc_i0, inc_i1, - output + output_offset, output_stride, normf, updated); - return; - } - } - - // Horizontal resampler, used only to calculate the support for - // short-circuiting. Since the resampling factor varies by band, - // the support also varies; use the largest resampling factor of - // any band to get the worst-case support. - int worstcase_band = anl.bandpass_bands_end() - 1; - RESAMPLER x_resampler(zoom_p2(x_xf, -anl.band_scale_exp(worstcase_band))); - int64_t xsi0, xsi1; - x_resampler.support(xi0, xi1, xsi0, xsi1); - - // Short-circuiting: if the image to be rendered falls entirely - // outside the data, just set it to zero instead of resampling down - // (potentially) high-resolution zeros to the display resolution. - // This makes a difference when zooming out by a large factor, for - // example such that the entire spectrogram falls within a single - // tile; that tile will necessarily be expensive to calculate, but - // the other tiles need not be, and mustn't be if we are to keep - // the total amount of work bounded by O(L) with respect to the - // signal length L regardless of zoom. - coef_index_t cxi0, cxi1; - get_band_coef_bounds(msc, worstcase_band, cxi0, cxi1); - if (ysi1 < 0 || // Entirely above - ysi0 >= anl.n_bands_total - 1 || // Entirely below - xsi1 < cxi0 || // Entirely to the left - xsi0 >= cxi1) // Entirely to the right - { - size_t n = (size_t)((yi1 - yi0) * (xi1 - xi0)); - for (size_t i = 0; i < n; i++) - output[i] = (T)0; - return; - } - - if (y_xf.a == 1 && y_xf.b == 0) { - // No Y resampling needed, render data resampled in the X - // direction only. - render_noyscale - (anl, msc, x_xf, xi0, xi1, yi0, yi1, - output, output_stride, normf); - } else { - // Construct a temporary image resampled in the - // X direction, but not yet in the Y direction. Include - // extra scanlines at the top and bottom for interpolation. - size_t n_pixels = (size_t)((ysi1 - ysi0) * (xi1 - xi0)); - pod_vector render_data(n_pixels); - - // Render data resampled in the X direction - RST *p = render_data.data(); - render_noyscale - (anl, msc, x_xf, xi0, xi1, - ysi0, ysi1, p, xi1 - xi0, normf); - - // Resample in the Y direction - for (int64_t xi = xi0; xi < xi1; xi++) { - transverse_source src(render_data.data(), - xi0, xi1, ysi0, ysi1, - xi); - stride_iterator dest(output + (xi - xi0), output_stride); - y_resampler(src, yi0, yi1, dest); - } - } - updated(xi0, xi1, yi0, yi1); -} - -template , - class RESAMPLER = lanczos2_pow2_resampler> -void render_p2scale(const analyzer &anl, - const coefs &msc, - int64_t xorigin, int64_t yorigin, - int64_t xi0, int64_t xi1, int xe, - int64_t yi0, int64_t yi1, int ye, - OI output, - NORMF normf = complex_abs_fob()) -{ - // Provide default inc_i0, inc_i1, and output_stride - render_incremental - (anl, msc, - affine_transform(ldexp(1, xe), xorigin), - affine_transform(ldexp(1, ye), yorigin), - xi0, xi1, yi0, yi1, - INT64_MIN, INT64_MAX, - output, xi1 - xi0, normf); -} - - -} // namespace - -#endif diff --git a/lib/gaborator/gaborator/resample2.h b/lib/gaborator/gaborator/resample2.h deleted file mode 100644 index c8b6e08..0000000 --- a/lib/gaborator/gaborator/resample2.h +++ /dev/null @@ -1,324 +0,0 @@ -// -// Fast resampling by powers of two -// -// Copyright (C) 2016-2021 Andreas Gustafsson. This file is part of -// the Gaborator library source distribution. See the file LICENSE at -// the top level of the distribution for license information. -// - -// Uses a two-lobe Lanczos kernel. Good enough for image data, not -// intended for audio. - -#ifndef _GABORATOR_RESAMPLE2_H -#define _GABORATOR_RESAMPLE2_H - -#include -#include -#include -#include // std::copy - -#include "gaborator/affine_transform.h" -#include "gaborator/pod_vector.h" - -namespace gaborator { - -// -// There are two ways to look at this. -// -// In one point of view, there is only one coordinate space, and -// coordinates are floating-point numbers. The various sub- and -// supersampled views differ in step sizes and the number of -// fractional coordinate bits, but any given coordinates refer to the -// same point in the image at any scale. Steps are powers of two, -// with integer exponents that may be negative. A step size > 1 -// implies downsampling (antialias lowpass filtering and subsampling), -// and a step size < 1 implies upsampling (aka interpolation). -// -// The coordinates are always integer multiples of the step size. -// -// e.g., -// x0 = 33.5, xstep = 0.5 -// x0 = 12, xstep = 4 -// -// In the other point of view, we introduce an integer exponent e and -// substitute x0 = i0 * 2^e and xstep = 2^e. Now instead of floating -// point coordinates, we use integer "indices". The above example -// now looks like his: -// -// i0 = 67, e = -1 -// i0 = 3, e = 2 -// -// This latter point of view is how the code actually works. -// - -// A power-of-two transform, as in y = 2^e x + origin - -struct p2_transform { - p2_transform(int e_, int64_t origin_): e(e_), origin(origin_) { } - // Convert a linear transform into a p2_transform - p2_transform(affine_transform xf) { - int exp; - double m = frexp(xf.a, &exp); - assert(m == 0.5); - e = exp - 1; - origin = xf.b; - assert(origin == xf.b); // No fraction - } - int e; - int64_t origin; -}; - -// Scale a transform by a power of two - -static inline p2_transform -zoom_p2(p2_transform xf, int e) { - return p2_transform(xf.e + e, xf.origin); -} - -static inline affine_transform -zoom_p2(affine_transform xf, int e) { - return affine_transform(ldexp(xf.a, e), xf.b); -} - -// Resample data from "source", generating a view between indices -// i0 and i1 of the scale determined by exponent e, and storing -// i1 - i0 samples starting at *out. -// -// The source must implement an operator() taking the arguments -// (int64_t i0, int64_t i1, T *out) and generating data for the base -// resolution (e=0). -// -// S is the type of the data source -// OI is the output iterator type - -template -OI resample2_p2xf(const S &source, p2_transform xf, - int64_t i0, int64_t i1, - bool interpolate, OI out) -{ - typedef typename std::iterator_traits::value_type T; - assert(i1 >= i0); - if (xf.e > 0) { - // Downsample - // Calculate super-octave coordinates - // margin is the support of the resampling kernel (on each side, - // not counting the center sample) - int margin = interpolate ? 1 : 0; - // When margin = 1, we use three samples, at 2i-1, 2i, 2i+1 - // and the corresponding half-open inverval is 2i-1...2i+1+1 - int64_t si0 = 2 * i0 - margin; - int64_t si1 = 2 * i1 + margin + 1; - // Get super-octave data - gaborator::pod_vector super_data(si1 - si0); - T *p = super_data.data(); - p = resample2_p2xf(source, p2_transform(xf.e - 1, xf.origin), - si0, si1, interpolate, p); - assert(p == super_data.data() + si1 - si0); - for (int64_t i = i0; i < i1; i++) { - int64_t si = 2 * i - si0; - T val; - if (!interpolate) { - // Point sampling - val = super_data[si]; - } else { - // Triangular kernel - T v1 = super_data[si - 1]; - T v0 = super_data[si]; - v1 += super_data[si + 1]; - val = - v0 * (T)0.5 + - v1 * (T)0.25; -#if 0 // Lanczos2 is overkill when downsampling. - } else { - // Two-lobe Lanczos kernel, needs margin = 2 - // Always aligned - T v3 = super_data[si - 3]; - // There is no v2 - T v1 = super_data[si - 1]; - T v0 = super_data[si]; - // There is still no v2 - v1 += super_data[si + 1]; - v3 += super_data[si + 3]; - val = - v0 * (T)0.49530706 + - v1 * (T)0.28388978 + - v3 * (T)-0.03154331; -#endif - } - *out++ = val; - } - } else if (xf.e < 0) { - // Upsample - if (! interpolate) { - // Return nearest neighbor. If the pixel lies - // exactly at the midpoint between the neighbors, - // return their average. - int sh = -xf.e; - int64_t si0 = i0 >> sh; - int64_t si1 = ((i1 - 1) >> sh) + 1 + 1; - gaborator::pod_vector source_data(si1 - si0); - source(xf.origin + si0, xf.origin + si1, source_data.data()); - for (int64_t i = i0; i < i1; i++) { - int64_t si = (i >> sh) - si0; - T val; - int rem = i & ((1 << sh) - 1); - int half = (1 << sh) >> 1; - if (rem < half) { - val = source_data[si]; - } else if (rem == half) { - val = (source_data[si] + source_data[si + 1]) * (T)0.5; - } else { // rem > half - val = source_data[si + 1]; - } - *out++ = val; - } - } else { - // Interpolate - // Calculate sub-octave coordinates - int margin = 2; - int64_t si0 = (i0 >> 1) - margin; - int64_t si1 = ((i1 - 1) >> 1) + margin + 1; - // Get sub-octave data - gaborator::pod_vector sub_data(si1 - si0); - T *p = sub_data.data(); - p = resample2_p2xf(source, p2_transform(xf.e + 1, xf.origin), - si0, si1, interpolate, p); - assert(p == sub_data.data() + si1 - si0); - for (int64_t i = i0; i < i1; i++) { - int64_t si = (i >> 1) - si0; - T val; - if (i & 1) { - T v1 = sub_data[si - 1]; - T v0 = sub_data[si]; - v0 += sub_data[si + 1]; - v1 += sub_data[si + 2]; - val = - v0 * (T)0.5625 + - v1 * (T)-0.0625; - } else { - val = sub_data[si]; - } - *out++ = val; - } - } - } else { - // e == 0 - out = source(xf.origin + i0, xf.origin + i1, out); - } - return out; -} - -// As above, but taking an affine_transform. - -template -OI resample2(const S &source, affine_transform lxf, - int64_t i0, int64_t i1, - bool interpolate, OI out) -{ - p2_transform xf(lxf); - typedef typename std::iterator_traits::value_type T; - gaborator::pod_vector data(i1 - i0); - T *p = data.data(); - p = resample2_p2xf(source, xf, i0, i1, interpolate, p); - return std::copy(data.data(), data.data() + (i1 - i0), out); -} - -// Calculate the range of source indices that will be accessed -// by a call to resample2(source, i0, i1, e) and return it -// through si0_ret and si1_ret. - -// XXX this should take an "interpolate" argument so we don't -// return an unnecessarily large support when interpolation is off - -inline void -resample2_support(affine_transform lxf, int64_t i0, int64_t i1, - int64_t &si0_ret, int64_t &si1_ret) -{ - p2_transform xf(lxf); - int margin = 2; - if (xf.e > 0) { - // Note code duplication wrt resample2_p2xf(). - // Also note tail recursion. - int64_t si0 = i0 * 2 - margin + 1; - int64_t si1 = i1 * 2 + margin; - resample2_support(zoom_p2(lxf, -1), - si0, si1, si0_ret, si1_ret); - } else if (xf.e < 0) { - int64_t si0 = (i0 >> 1) - margin; - int64_t si1 = ((i1 - 1) >> 1) + margin + 1; - resample2_support(zoom_p2(lxf, +1), - si0, si1, si0_ret, si1_ret); - } else { - si0_ret = xf.origin + i0; - si1_ret = xf.origin + i1; - } -} - -// Inverse of the above, more or less: calculate the range of -// destination indices that depend on a given range of source indices. - -inline void -resample2_inv_support(affine_transform lxf, int64_t si0, int64_t si1, - int64_t &i0_ret, int64_t &i1_ret) -{ - p2_transform xf(lxf); - // Conservative - int margin = 2; - if (xf.e > 0) { - int64_t di0, di1; - resample2_inv_support(zoom_p2(lxf, -1), - si0, si1, di0, di1); - i0_ret = di0 >> 1; - i1_ret = (di1 >> 1) + 1; - } else if (xf.e < 0) { - int64_t di0, di1; - resample2_inv_support(zoom_p2(lxf, +1), - si0, si1, di0, di1); - i0_ret = di0 * 2 - margin + 1; - i1_ret = di1 * 2 + margin; - } else { - i0_ret = si0 - xf.origin; - i1_ret = si1 - xf.origin; - } -} - -// Class wrappers for compatibility with other resamplers. - -// Lanczos2 power-of-two resampler - -struct lanczos2_pow2_resampler { - lanczos2_pow2_resampler(affine_transform xf_): xf(xf_) { } - template - OI operator()(const S &source, int64_t i0, int64_t i1, OI out) const { - return resample2(source, xf, i0, i1, true, out); - } - void support(int64_t i0, int64_t i1, int64_t &si0_ret, int64_t &si1_ret) const { - return resample2_support(xf, i0, i1, si0_ret, si1_ret); - } - void inv_support(int64_t si0, int64_t si1, int64_t &i0_ret, int64_t &i1_ret) { - return resample2_inv_support(xf, si0, si1, i0_ret, i1_ret); - } - affine_transform xf; -}; - -// Nearest-neighbor power-of-two resampler -// XXX simplify - -struct nn_pow2_resampler { - nn_pow2_resampler(affine_transform xf_): xf(xf_){ } - template - OI operator()(const S &source, int64_t i0, int64_t i1, OI out) const { - return resample2(source, xf, i0, i1, false, out); - } - void support(int64_t i0, int64_t i1, int64_t &si0_ret, int64_t &si1_ret) const { - return resample2_support(xf, i0, i1, si0_ret, si1_ret); - } - void inv_support(int64_t si0, int64_t si1, int64_t &i0_ret, int64_t &i1_ret) { - return resample2_inv_support(xf, si0, si1, i0_ret, i1_ret); - } - affine_transform xf; -}; - -} // namespace - -#endif diff --git a/lib/gaborator/gaborator/vector_math.h b/lib/gaborator/gaborator/vector_math.h deleted file mode 100644 index f1d9cca..0000000 --- a/lib/gaborator/gaborator/vector_math.h +++ /dev/null @@ -1,301 +0,0 @@ -// -// Vector math operations -// -// Copyright (C) 2016-2018 Andreas Gustafsson. This file is part of -// the Gaborator library source distribution. See the file LICENSE at -// the top level of the distribution for license information. -// - -#ifndef _GABORATOR_VECTOR_MATH_H -#define _GABORATOR_VECTOR_MATH_H - -#include - -#if GABORATOR_USE_SSE3_INTRINSICS -#include -#endif - -#include - -namespace gaborator { - -// The _naive versions are used when SSE is not available, and as -// references when testing the SSE versions - -// Naive or not, this is faster than the macOS std::complex -// multiplication operator, which checks for infinities even with -// -ffast-math. - -template -std::complex complex_mul_naive(std::complex a, - std::complex b) -{ - return std::complex(a.real() * b.real() - a.imag() * b.imag(), - a.real() * b.imag() + a.imag() * b.real()); -} - -#if GABORATOR_USE_SSE3_INTRINSICS -// Note: this is sometimes slower than the naive code. -static inline -std::complex complex_mul(std::complex a_, - std::complex b_) -{ - __v2df a = _mm_setr_pd(a_.real(), a_.imag()); - __v2df b = _mm_setr_pd(b_.real(), b_.imag()); - __v2df as = (__v2df) _mm_shuffle_pd(a, a, 0x1); - __v2df t0 = _mm_mul_pd(a, _mm_shuffle_pd(b, b, 0x0)); - __v2df t1 = _mm_mul_pd(as, _mm_shuffle_pd(b, b, 0x3)); - __v2df c = __builtin_ia32_addsubpd(t0, t1); // SSE3 - return std::complex(c[0], c[1]); -} -#else -static inline -std::complex complex_mul(std::complex a_, - std::complex b_) -{ - return complex_mul_naive(a_, b_); -} -#endif - -static inline -std::complex complex_mul(std::complex a_, - std::complex b_) -{ - return complex_mul_naive(a_, b_); -} - -template -static inline void -elementwise_product_naive(T *r, - U *a, - V *b, - size_t n) -{ - for (size_t i = 0; i < n; i++) - r[i] = complex_mul(a[i], b[i]); -} - -template -static inline void -elementwise_product_times_scalar_naive(T *r, - U *a, - V *b, - S s, - size_t n) -{ - for (size_t i = 0; i < n; i++) - r[i] = a[i] * b[i] * s; -} - -// I is the input complex data type, O is the output data type -template -static inline void -complex_magnitude_naive(I *inv, - O *outv, - size_t n) -{ - for (size_t i = 0; i < n; i++) - outv[i] = std::sqrt(inv[i].real() * inv[i].real() + inv[i].imag() * inv[i].imag()); -} - -#if GABORATOR_USE_SSE3_INTRINSICS - -#include - -// Perform two complex float multiplies in parallel - -static inline -__v4sf complex_mul_vec2(__v4sf aa, __v4sf bb) { - __v4sf aas =_mm_shuffle_ps(aa, aa, 0xb1); - __v4sf t0 = _mm_mul_ps(aa, _mm_moveldup_ps(bb)); - __v4sf t1 = _mm_mul_ps(aas, _mm_movehdup_ps(bb)); - return __builtin_ia32_addsubps(t0, t1); // SSE3 -} - -// Calculate the elementwise product of a complex float vector -// and another complex float vector. - -static inline void -elementwise_product(std::complex *cv, - const std::complex *av, - const std::complex *bv, - size_t n) -{ - assert((n & 1) == 0); - n >>= 1; - __v4sf *c = (__v4sf *) cv; - const __v4sf *a = (const __v4sf *) av; - const __v4sf *b = (const __v4sf *) bv; - while (n--) { - __v4sf aa = *a++; - __v4sf bb = *b++; - *c++ = complex_mul_vec2(aa, bb); - } -} - -// Calculate the elementwise product of a complex float vector -// and real float vector. -// -// The input "a" may be unaligned; the output "c" must be aligned. - -static inline void -elementwise_product(std::complex *cv, - const std::complex *av, - const float *bv, - size_t n) -{ - assert((n & 3) == 0); - n >>= 2; - __v4sf *c = (__v4sf *) cv; - const __v4sf *a = (const __v4sf *) av; - const __v4sf *b = (const __v4sf *) bv; - while (n--) { - __v4sf a0 = (__v4sf) _mm_loadu_si128((const __m128i *) a++); - __v4sf a1 = (__v4sf) _mm_loadu_si128((const __m128i *) a++); - __v4sf bb = *b++; - *c++ = _mm_mul_ps(a0, _mm_unpacklo_ps(bb, bb)); - *c++ = _mm_mul_ps(a1, _mm_unpackhi_ps(bb, bb)); - } -} - -static inline void -elementwise_product_times_scalar(std::complex *cv, - const std::complex *av, - const std::complex *bv, - std::complex d, - size_t n) -{ - assert((n & 1) == 0); - n >>= 1; - const __v4sf *a = (const __v4sf *) av; - const __v4sf *b = (const __v4sf *) bv; - const __v4sf dd = (__v4sf) { d.real(), d.imag(), d.real(), d.imag() }; - __v4sf *c = (__v4sf *) cv; - while (n--) { - __v4sf aa = *a++; - __v4sf bb = *b++; - *c++ = complex_mul_vec2(complex_mul_vec2(aa, bb), dd); - } -} - -// XXX arguments reversed wrt others - -static inline void -complex_magnitude(std::complex *inv, - float *outv, - size_t n) -{ - // Processes four complex values (32 bytes) at a time , - // outputting four scalar magnitudes (16 bytes) at a time. - while ((((uintptr_t) inv) & 0x1F) && n) { - std::complex v = *inv++; - *outv++ = std::sqrt(v.real() * v.real() + v.imag() * v.imag()); - n--; - } - const __v4sf *in = (const __v4sf *) inv; - __v4sf *out = (__v4sf *) outv; - while (n >= 4) { - __v4sf aa = *in++; // c0re c0im c1re c1im - __v4sf aa2 = _mm_mul_ps(aa, aa); // c0re^2 c0im^2 c1re^2 c1im^2 - __v4sf bb = *in++; // c2re c2im c3re c3im - __v4sf bb2 = _mm_mul_ps(bb, bb); // etc - // Gather the real parts: x0 x2 y0 y2 - // 10 00 10 00 = 0x88 - __v4sf re2 =_mm_shuffle_ps(aa2, bb2, 0x88); - __v4sf im2 =_mm_shuffle_ps(aa2, bb2, 0xdd); - __v4sf mag2 = _mm_add_ps(re2, im2); - __v4sf mag = __builtin_ia32_sqrtps(mag2); - // Unaligned store - _mm_storeu_si128((__m128i *)out, (__m128i)mag); - out++; - n -= 4; - } - inv = (std::complex *) in; - outv = (float *)out; - while (n) { - std::complex v = *inv++; - *outv++ = std::sqrt(v.real() * v.real() + v.imag() * v.imag()); - n--; - } -} - -// Double-precision version is unoptimized - -static inline void -elementwise_product(std::complex *c, - const std::complex *a, - const std::complex *b, - size_t n) -{ - elementwise_product_naive(c, a, b, n); -} - -static inline void -elementwise_product(std::complex *c, - const std::complex *a, - const double *b, - size_t n) -{ - elementwise_product_naive(c, a, b, n); -} - -template -static inline void -elementwise_product_times_scalar(T *r, - U *a, - V *b, - S s, - size_t n) -{ - elementwise_product_times_scalar_naive(r, a, b, s, n); -} - -template -static inline void -complex_magnitude(std::complex *inv, - O *outv, - size_t n) -{ - complex_magnitude_naive(inv, outv, n); -} - -#else // ! GABORATOR_USE_SSE3_INTRINSICS - -// Forward to the naive implementations. These are inline functions -// rather than #defines to avoid namespace pollution. - -template -static inline void -elementwise_product(T *r, - U *a, - V *b, - size_t n) -{ - elementwise_product_naive(r, a, b, n); -} - -template -static inline void -elementwise_product_times_scalar(T *r, - U *a, - V *b, - S s, - size_t n) -{ - elementwise_product_times_scalar_naive(r, a, b, s, n); -} - -template -static inline void -complex_magnitude(I *inv, - O *outv, - size_t n) -{ - complex_magnitude_naive(inv, outv, n); -} - -#endif // ! USE_SSE3_INTINSICS - -} // namespace - -#endif diff --git a/lib/gaborator/gaborator/version.h b/lib/gaborator/gaborator/version.h deleted file mode 100644 index 41ddc8a..0000000 --- a/lib/gaborator/gaborator/version.h +++ /dev/null @@ -1,15 +0,0 @@ -// -// Version number -// -// Copyright (C) 2015-2021 Andreas Gustafsson. This file is part of -// the Gaborator library source distribution. See the file LICENSE at -// the top level of the distribution for license information. -// - -#ifndef _GABORATOR_VERSION_H -#define _GABORATOR_VERSION_H - -#define GABORATOR_VERSION_MAJOR 1 -#define GABORATOR_VERSION_MINOR 7 - -#endif