Removed gaborator source, replace with submodule
All checks were successful
continuous-integration/drone/push Build is passing
All checks were successful
continuous-integration/drone/push Build is passing
This commit is contained in:
parent
21a76a98ed
commit
430ac5d89f
3
.gitmodules
vendored
3
.gitmodules
vendored
|
@ -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
|
||||
|
|
1
lib/gaborator
Submodule
1
lib/gaborator
Submodule
|
@ -0,0 +1 @@
|
|||
Subproject commit 80ef6d4c8703d165660f7ef8d321fce9b526adfc
|
|
@ -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
|
|
@ -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.
|
|
@ -1 +0,0 @@
|
|||
See doc/index.html for HTML documentation.
|
|
@ -1,661 +0,0 @@
|
|||
GNU AFFERO GENERAL PUBLIC LICENSE
|
||||
Version 3, 19 November 2007
|
||||
|
||||
Copyright (C) 2007 Free Software Foundation, Inc. <https://fsf.org/>
|
||||
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.
|
||||
|
||||
<one line to give the program's name and a brief idea of what it does.>
|
||||
Copyright (C) <year> <name of author>
|
||||
|
||||
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 <https://www.gnu.org/licenses/>.
|
||||
|
||||
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
|
||||
<https://www.gnu.org/licenses/>.
|
|
@ -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;
|
||||
}
|
Binary file not shown.
Before Width: | Height: | Size: 19 KiB |
|
@ -1,266 +0,0 @@
|
|||
<!DOCTYPE html>
|
||||
<!--
|
||||
Copyright (C) 2017-2020 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.
|
||||
-->
|
||||
<html>
|
||||
<head>
|
||||
<link rel="stylesheet" href="doc.css" />
|
||||
<title>Gaborator Example 2: Frequency-Domain Filtering</title>
|
||||
</head>
|
||||
<body>
|
||||
<h1>Example 2: Frequency-Domain Filtering</h1>
|
||||
|
||||
<h2>Introduction</h2>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<p>The specific filter implemented here is a 3 dB/octave lowpass
|
||||
filter. This is sometimes called a <i>pinking filter</i> 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.
|
||||
</p>
|
||||
|
||||
<p>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.
|
||||
</p>
|
||||
|
||||
<h2>Preamble</h2>
|
||||
|
||||
<pre>
|
||||
#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);
|
||||
}
|
||||
</pre>
|
||||
|
||||
<h2>Reading the Audio</h2>
|
||||
|
||||
<p>The code for reading the input audio file is identical to
|
||||
that in <a href="render.html">Example 1</a>:</p>
|
||||
|
||||
<pre>
|
||||
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);
|
||||
</pre>
|
||||
|
||||
<h2>Spectrum Analysis Parameters</h2>
|
||||
|
||||
<p>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.</p>
|
||||
<pre>
|
||||
gaborator::parameters params(100, 20.0 / fs);
|
||||
gaborator::analyzer<float> analyzer(params);
|
||||
</pre>
|
||||
|
||||
<h2>Precalculating Gains</h2>
|
||||
|
||||
<p>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 <code>band_gains</code>
|
||||
of one gain value per band, including one for the
|
||||
special lowpass band that contains the frequencies from 0 to 20 Hz.</p>
|
||||
|
||||
<pre>
|
||||
std::vector<float> band_gains(analyzer.bands_end());
|
||||
</pre>
|
||||
|
||||
<p>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
|
||||
<code>analyzer</code> method <code>band_ff()</code>, which
|
||||
returns the center frequency of the band in units of the
|
||||
sampling frequency. The gain is normalized to unity at 20 Hz.
|
||||
</p>
|
||||
<pre>
|
||||
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);
|
||||
}
|
||||
</pre>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<pre>
|
||||
band_gains[analyzer.band_lowpass()] = band_gains[analyzer.bandpass_bands_end() - 1];
|
||||
</pre>
|
||||
|
||||
<h2>De-interleaving</h2>
|
||||
|
||||
<p>To handle stereo and other multi-channel audio files,
|
||||
we will loop over the channels and filter each channel separately.
|
||||
Since <i>libsndfile</i> produces interleaved samples, we first
|
||||
de-interleave the current channel into a temporary vector called
|
||||
<code>channel</code>:</p>
|
||||
<pre>
|
||||
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];
|
||||
</pre>
|
||||
<h2>Spectrum Analysis</h2>
|
||||
<p>Now we can spectrum analyze the current channel, producing
|
||||
a set of coefficients:</p>
|
||||
<pre>
|
||||
gaborator::coefs<float> coefs(analyzer);
|
||||
analyzer.analyze(channel.data(), 0, channel.size(), coefs);
|
||||
</pre>
|
||||
|
||||
<h2>Filtering</h2>
|
||||
<p>
|
||||
The filtering is done using the function
|
||||
<code>process()</code>, 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 <code>int64_t</code> argument is the time in units
|
||||
of samples; this could be use to implement a time-varying filter if
|
||||
desired.</p>
|
||||
<p>
|
||||
The second and third argument to <code>process()</code> specify a
|
||||
range of frequency bands to process; here we pass <code>INT_MIN,
|
||||
INT_MAX</code> to process all of them. Similarly, the fourth and
|
||||
fifth argument specify a time range to process, and we pass
|
||||
<code>INT64_MIN, INT64_MAX</code> to process all the coefficients
|
||||
in <code>coefs</code> regardless of time.
|
||||
</p>
|
||||
<pre>
|
||||
process([&](int band, int64_t, std::complex<float> &coef) {
|
||||
coef *= band_gains[band];
|
||||
},
|
||||
INT_MIN, INT_MAX,
|
||||
INT64_MIN, INT64_MAX,
|
||||
coefs);
|
||||
</pre>
|
||||
|
||||
<h2>Resynthesis</h2>
|
||||
<p>We can now resynthesize audio from the filtered coefficients by
|
||||
calling <code>synthesize()</code>. This is a mirror image of the call to
|
||||
<code>analyze()</code>: now the coefficients are the input, and
|
||||
the buffer of samples is the output. The <code>channel</code>
|
||||
vector that originally contained the input samples for the channel
|
||||
is now reused to hold the output samples.</p>
|
||||
<pre>
|
||||
analyzer.synthesize(coefs, 0, channel.size(), channel.data());
|
||||
</pre>
|
||||
|
||||
<h2>Re-interleaving</h2>
|
||||
<p>The <code>audio</code> vector that contained the
|
||||
original interleaved audio is reused for the interleaved
|
||||
filtered audio. This concludes the loop over the channels.
|
||||
</p>
|
||||
<pre>
|
||||
for (sf_count_t i = 0; i < n_frames; i++)
|
||||
audio[i * sfinfo.channels + ch] = channel[i];
|
||||
}
|
||||
</pre>
|
||||
|
||||
<h2>Writing the Audio</h2>
|
||||
<p>The filtered audio is written using <i>libsndfile</i>,
|
||||
using code that closely mirrors that for reading.
|
||||
Note that we use <code>SFC_SET_CLIPPING</code>
|
||||
to make sure that any samples too loud for the file format
|
||||
will saturate; by default, <i>libsndfile</i> makes them
|
||||
wrap around, which sounds really bad.</p>
|
||||
<a name="writing_audio_code">
|
||||
<pre>
|
||||
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);
|
||||
</pre>
|
||||
</a>
|
||||
|
||||
<h2>Postamble</h2>
|
||||
<p>
|
||||
We need a couple more lines of boilerplate to make the example a
|
||||
complete program:
|
||||
</p>
|
||||
<pre>
|
||||
return 0;
|
||||
}
|
||||
</pre>
|
||||
|
||||
<h2>Compiling</h2>
|
||||
<p>Like <a href="render.html#compiling">Example 1</a>, this example
|
||||
can be built using a one-line build command:
|
||||
</p>
|
||||
<pre class="build Darwin Linux NetBSD FreeBSD">
|
||||
c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` filter.cc `pkg-config --libs sndfile` -o filter
|
||||
</pre>
|
||||
<p>Or using the vDSP FFT on macOS:</p>
|
||||
<pre class="build Darwin">
|
||||
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
|
||||
</pre>
|
||||
<p>Or using PFFFT (see <a href="render.html#compiling">Example 1</a> for how to download and build PFFFT):</p>
|
||||
<pre class="build">
|
||||
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
|
||||
</pre>
|
||||
|
||||
<h2>Running</h2>
|
||||
<p>Running the following shell commands will download an example
|
||||
audio file containing five seconds of white noise and filter it,
|
||||
producing pink noise.</p>
|
||||
<pre class="run">
|
||||
wget http://download.gaborator.com/audio/white_noise.wav
|
||||
./filter white_noise.wav pink_noise.wav
|
||||
</pre>
|
||||
|
||||
<h2>Frequency response</h2>
|
||||
<p>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:</p>
|
||||
<img src="filter-response.png" alt="Frequency response plot" data-autogen="no">
|
||||
|
||||
<div class="nav"><span class="prev"><a href="render.html">Previous: Example 1: Rendering a Spectrogram Image</a></span><span class="next"><a href="stream.html">Next: Example 3: Streaming</a></span></div>
|
||||
|
||||
</body>
|
||||
</html>
|
Binary file not shown.
Before Width: | Height: | Size: 122 KiB |
Binary file not shown.
Before Width: | Height: | Size: 107 KiB |
Binary file not shown.
Before Width: | Height: | Size: 14 KiB |
|
@ -1,78 +0,0 @@
|
|||
<!DOCTYPE html>
|
||||
<!--
|
||||
Copyright (C) 2018-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.
|
||||
-->
|
||||
<html>
|
||||
<head>
|
||||
<link rel="stylesheet" href="doc.css" />
|
||||
<link rel="icon" type="image/png" href="favicon64.png" />
|
||||
<title>The Gaborator</title>
|
||||
</head>
|
||||
<body>
|
||||
<h1>The Gaborator</h1>
|
||||
<p>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.</p>
|
||||
|
||||
<p>The Gaborator implements the invertible constant-Q transform of
|
||||
Velasco, Holighaus, Dörfler, and Grill, described in the papers
|
||||
<i><a href="http://www.univie.ac.at/nonstatgab/pdf_files/dohogrve11_amsart.pdf">
|
||||
Constructing an invertible constant-Q transform with nonstationary Gabor frames, 2011</a></i>
|
||||
and <i><a href="http://www.univie.ac.at/nonstatgab/pdf_files/dogrhove12_amsart.pdf">
|
||||
A Framework for invertible, real-time constant-Q transforms, 2012</a></i>,
|
||||
using Gaussian bandpass filters and an efficient multi-rate architecture.
|
||||
</p>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<p>The Gaborator is open source under the GNU Affero General Public
|
||||
License, version 3, and is also available for commercial licensing.
|
||||
See the file <a href="../LICENSE">LICENSE</a> for details.</p>
|
||||
|
||||
<h2>Example Code</h2>
|
||||
|
||||
<p>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 <code>.cc</code> file in
|
||||
the <code>examples/</code> directory.</p>
|
||||
<ul>
|
||||
<li><a href="render.html">Example 1: Rendering a Spectrogram Image</a></li>
|
||||
<li><a href="filter.html">Example 2: Frequency-Domain Filtering</a></li>
|
||||
<li><a href="stream.html">Example 3: Streaming</a></li>
|
||||
<li><a href="snr.html">Example 4: Measuring the Signal-to-Noise Ratio</a></li>
|
||||
<li><a href="synth.html">Example 5: Synthesis from Scratch</a></li>
|
||||
</ul>
|
||||
|
||||
<h2>API Reference</h2>
|
||||
<p>The following documents define the library API.
|
||||
</p>
|
||||
<ul>
|
||||
<li><a href="ref/intro.html">API Introduction</a></li>
|
||||
<li><a href="ref/gaborator_h.html">Spectrum analysis and synthesis: <code>gaborator.h</code></a></li>
|
||||
<li><a href="ref/render_h.html">Spectrogram rendering: <code>render.h</code></a></li>
|
||||
</ul>
|
||||
|
||||
<h2>How it Works</h2>
|
||||
<p>The following document outlines the operation of the library.</p>
|
||||
<ul>
|
||||
<lI><a href="overview.html">Overview of Operation</a>
|
||||
</ul>
|
||||
|
||||
<h2>FAQ</h2>
|
||||
<ul>
|
||||
<li><a href="realtime.html">Is it real-time?</a></li>
|
||||
</ul>
|
||||
|
||||
<h2>Contact</h2>
|
||||
<p>Email questions and bug reports to the author at info@gaborator.com.</p>
|
||||
|
||||
</body>
|
||||
</html>
|
|
@ -1,135 +0,0 @@
|
|||
<!DOCTYPE html>
|
||||
<!--
|
||||
Copyright (C) 2018-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.
|
||||
-->
|
||||
<html>
|
||||
<head>
|
||||
<link rel="stylesheet" href="doc.css" type="text/css" />
|
||||
<title>Overview of Operation</title>
|
||||
</head>
|
||||
<body>
|
||||
<h1>Overview of Operation</h1>
|
||||
|
||||
<p>The Gaborator performs three main functions:</p>
|
||||
<ul>
|
||||
<li>spectrum <i>analysis</i>, which turns a signal into a set
|
||||
of <i>spectrogram coefficients</i>
|
||||
<li><i>resynthesis</i> (aka <i>reconstruction</i>), which turns a
|
||||
set of coefficients back into a signal, and
|
||||
<li><i>rendering</i>, which
|
||||
turns a set of coefficients into a rectangular array of
|
||||
amplitude values that can be turned into pixels to display
|
||||
a spectrogram.
|
||||
</ul>
|
||||
|
||||
<p>The following sections give a high-level overview of each
|
||||
of these functions.</p>
|
||||
|
||||
<h2>Analysis</h2>
|
||||
|
||||
<p>The first step of the analysis is to run the signal through
|
||||
an <i>analysis filter bank</i>, to split it into a number of
|
||||
overlapping frequency <i>bands</i>.</p>
|
||||
|
||||
<p>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 <i>constant-Q</i> 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 <i>minimum frequency</i> or f<sub>min</sub>.
|
||||
</p>
|
||||
|
||||
<p>Although frequencies below f<sub>min</sub> 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.</P>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<p>The following plot shows the frequency responses of the analysis
|
||||
filters at 12 bands per octave and f<sub>min</sub> = 0.03. A more
|
||||
typical f<sub>min</sub> 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.</p>
|
||||
|
||||
<img src="gen/allkernels_v1_bpo12_ffmin0.03_ffref0.5_anl_wob.png" alt="Analysis filters">
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<p>The center frequencies of the analysis filters and the points in
|
||||
time at which they are sampled form a two-dimensional,
|
||||
multi-resolution <i>time-frequency grid</i>, where high frequencies
|
||||
are sampled sparsely in frequency but densely in time, and low
|
||||
frequencies are sampled densely in frequency but sparsely in time.</p>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<img src="gen/grid_v1_bpo12_ffmin0.03_ffref0.5_wob.png" alt="Sampling grid">
|
||||
|
||||
<h2>Resynthesis</h2>
|
||||
|
||||
<p>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 <i>reconstruction filter bank</i>
|
||||
that is a <i>dual</i> of the analysis filter bank. The following
|
||||
plot shows the frequency responses of the reconstruction filters
|
||||
corresponding to the analysis filters shown earlier.</p>
|
||||
|
||||
<img src="gen/allkernels_v1_bpo12_ffmin0.03_ffref0.5_syn_wob.png" alt="Reconstruction filters">
|
||||
|
||||
<p>Although the bandpass filters may look similar to the Gaussian
|
||||
filters of the analysis filter bank, their shapes are actually subtly
|
||||
different.</p>
|
||||
|
||||
<h2>Spectrogram Rendering</h2>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<p>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).</p>
|
||||
|
||||
<div class="nav"><span class="prev"><a href="ref/render_h.html">Previous: Spectrogram rendering: <code>render.h</code></a></span><span class="next"><a href="realtime.html">Next: Is it real-time?</a></span></div>
|
||||
|
||||
</body>
|
||||
</html>
|
|
@ -1,126 +0,0 @@
|
|||
<!DOCTYPE html>
|
||||
<!--
|
||||
Copyright (C) 2018-2020 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.
|
||||
-->
|
||||
<html>
|
||||
<head>
|
||||
<link rel="stylesheet" href="doc.css" />
|
||||
<title>Is it real-time?</title>
|
||||
</head>
|
||||
<body>
|
||||
<h1>Is it real-time?</h1>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<h2>Can it processes a recording in less time than its duration?</h2>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<h2>Does it have bounded latency?
|
||||
Can it start producing output before consuming the entire input?
|
||||
Will it stream?</h2>
|
||||
<p>Yes. See the <a href="stream.html">streaming example</a>.
|
||||
|
||||
<h2>Does it have low latency?</h2>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<p>The latency only affects <i>causal</i> 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.</p>
|
||||
|
||||
<p>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 <a href="https://waxingwave.com/spectrolite/">Spectrolite</a>
|
||||
iOS app.</p>
|
||||
|
||||
<h2>Does it support small blocks sizes?</h2>
|
||||
|
||||
<p>Yes, but there is a significant performance penalty.
|
||||
The Gaborator works most efficiently when the signal is processed
|
||||
in large blocks, preferably 2<sup>17</sup> samples or more,
|
||||
corresponding to several seconds of signal at typical audio sample
|
||||
rates.</p>
|
||||
|
||||
<p>A real-time application aiming for low latency will want to
|
||||
use smaller blocks, for examples 2<sup>5</sup> to 2<sup>10</sup>
|
||||
samples, and processing these will be significantly slower.
|
||||
For example, as of version 1.4, analyzing a signal in blocks of
|
||||
2<sup>10</sup> samples takes roughly five times as much CPU as
|
||||
analyzing it in blocks of 2<sup>20</sup> samples.</p>
|
||||
|
||||
<p>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 2<sup>4</sup> = 16 samples.</p>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<h2>Can it process a signal stream of any length?</h2>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<h2>Does it avoid dynamic memory allocation in the audio processing path?</h2>
|
||||
|
||||
<p>Currently, no — it dynamically allocates both the coefficient data
|
||||
structures and various temporary buffers.</p>
|
||||
|
||||
<div class="nav"><span class="prev"><a href="overview.html">Previous: Overview of Operation</a></span></div>
|
||||
|
||||
</body>
|
||||
</html>
|
|
@ -1,462 +0,0 @@
|
|||
<!DOCTYPE html>
|
||||
<!--
|
||||
Copyright (C) 2018-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.
|
||||
-->
|
||||
<html>
|
||||
<head>
|
||||
<link rel="stylesheet" href="../doc.css" type="text/css" />
|
||||
<title>Gaborator reference: gaborator.h</title>
|
||||
</head>
|
||||
<body>
|
||||
<h1>Gaborator reference: <code>gaborator.h</code></h1>
|
||||
|
||||
<h2>Spectrum Analysis Parameters</h2>
|
||||
|
||||
<p>A <code>parameters</code> object holds a set of parameters that
|
||||
determine the frequency range and resolution of the spectrum
|
||||
analysis.</p>
|
||||
|
||||
<pre>
|
||||
class parameters {
|
||||
</pre>
|
||||
|
||||
<div class="class_def">
|
||||
<h3>Constructor</h3>
|
||||
<pre>
|
||||
parameters(unsigned int bands_per_octave,
|
||||
double ff_min,
|
||||
double ff_ref = 1.0);
|
||||
</pre>
|
||||
<dl>
|
||||
<dt><code>bands_per_octave</code></dt>
|
||||
<dd>The number of frequency bands per octave.
|
||||
Values from 4 to 384 (inclusive) are supported.
|
||||
</dd>
|
||||
<dt><code>ff_min</code></dt>
|
||||
<dd>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 <code>ff_min</code> falls between the two lowest
|
||||
frequency bandpass filters.
|
||||
Values from 0.001 to 0.13 are supported.</dd>
|
||||
<dt><code>ff_ref</code></dt>
|
||||
<dd>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 <code>ff_ref</code>. If <code>ff_ref</code>
|
||||
falls outside the frequency range of the bandpass filter bank, this
|
||||
works as if the range were extended to include
|
||||
<code>ff_ref</code>. Must be positive. A typical value
|
||||
when analyzing music is <code>440.0 / fs</code>, where
|
||||
<code>fs</code> is the sample rate in Hz.
|
||||
</dd>
|
||||
</dl>
|
||||
<h3>Comparison</h3>
|
||||
<p>
|
||||
Comparison operators are provided for compatibility with
|
||||
standard container classes. The ordering is arbitrary but consistent.
|
||||
</p>
|
||||
<pre>
|
||||
bool operator<(const parameters &rhs) const;
|
||||
bool operator==(const parameters &rhs) const;
|
||||
</pre>
|
||||
|
||||
</div>
|
||||
<pre>
|
||||
};
|
||||
</pre>
|
||||
|
||||
<h2>Spectrogram Coefficients</h2>
|
||||
|
||||
<pre class="forward_decl">
|
||||
template<class T> class analyzer;
|
||||
</pre>
|
||||
|
||||
<p>
|
||||
A <code>coefs</code> 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 <code>T</code>
|
||||
must match that of the <code>analyzer</code> (usually <code>float</code>).
|
||||
The template argument <code>C</code> is the data type used to store each
|
||||
coefficient value; there is usually no need to specify it explicitly as
|
||||
it will default to <code>std::complex<T></code>.
|
||||
</p>
|
||||
|
||||
<pre>
|
||||
template<class T, class C = std::complex<T>>
|
||||
class coefs {
|
||||
</pre>
|
||||
<div class="class_def">
|
||||
<h3>Constructor</h3>
|
||||
<pre>
|
||||
coefs(analyzer<T> &a);
|
||||
</pre>
|
||||
<p>
|
||||
Construct an empty set of coefficients for use with the spectrum
|
||||
analyzer <code>a</code>. This represents a signal that is zero
|
||||
at all points in time.
|
||||
</p>
|
||||
|
||||
</div>
|
||||
<pre>
|
||||
};
|
||||
</pre>
|
||||
|
||||
<h2>Spectrum Analyzer</h2>
|
||||
|
||||
<p>
|
||||
The <code>analyzer</code> object performs spectrum analysis and/or resynthesis
|
||||
according to the given parameters. The template argument <code>T</code> is
|
||||
the floating-point type to use for the calculations. This is typically <code>float</code>;
|
||||
alternatively, <code>double</code> can be used for increased accuracy at the
|
||||
expense of speed and memory consumption.</p>
|
||||
|
||||
<pre>template<class T>
|
||||
class analyzer {</pre>
|
||||
<div class="class_def">
|
||||
|
||||
<h3>Constructor</h3>
|
||||
|
||||
<pre>
|
||||
analyzer(const parameters &params);
|
||||
</pre>
|
||||
<dl>
|
||||
<dt><code>params</code></dt>
|
||||
<dd>The spectrum analysis parameters.
|
||||
</dl>
|
||||
|
||||
<h3>Analysis and synthesis</h3>
|
||||
|
||||
<pre>
|
||||
void
|
||||
analyze(const T *signal,
|
||||
int64_t t0,
|
||||
int64_t t1,
|
||||
coefs<T> &coefs) const;
|
||||
</pre>
|
||||
<p>Spectrum analyze the samples at <code>*signal</code> and add the
|
||||
resulting coefficients to <code>coefs</code>.
|
||||
<dl>
|
||||
<dt><code>signal</code></dt>
|
||||
<dd>The signal samples to analyze, beginning with the sample from time <code>t0</code>
|
||||
and ending with the last sample before time <code>t1</code>, for a total of
|
||||
<code>t1 - t0</code> samples.
|
||||
<dt><code>t0</code></dt>
|
||||
<dd>The point in time when the sample at <code>signal[0]</code> 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 ±10<sup>8</sup> samples, so using
|
||||
large time values should be avoided when they are not necessary because
|
||||
of the length of the track.
|
||||
</dd>
|
||||
<dt><code>t1</code></dt>
|
||||
<dd>The point in time of the sample one past the
|
||||
end of the array of samples at <code>signal</code>,
|
||||
in samples.
|
||||
</dd>
|
||||
<dt><code>coefs</code></dt><dd>The coefficient object that the results of the
|
||||
spectrum analysis are added to.
|
||||
</dl>
|
||||
<p>If the <code>coefs</code> 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 <code>analyze()</code> with non-overlapping ranges that together
|
||||
cover the entire signal. For efficiency, the blocks should
|
||||
be large, as in
|
||||
<code>analyze(first_131072_samples, 0, 131072, coefs)</code>,
|
||||
<code>analyze(next_131072_samples, 131072, 262144, coefs)</code>,
|
||||
etc.
|
||||
</p>
|
||||
|
||||
<pre>
|
||||
void
|
||||
synthesize(const coefs<T> &coefs,
|
||||
uint64_t t0,
|
||||
uint64_t t1,
|
||||
T *signal) const;
|
||||
</pre>
|
||||
<p>Synthesize signal samples from the coefficients <code>coef</code> and store
|
||||
them at <code>*signal</code>.
|
||||
</p>
|
||||
<dl>
|
||||
<dt><code>coefs</code></dt><dd>The coefficients to synthesize the signal from.</dd>
|
||||
<dt><code>t0</code></dt>
|
||||
<dd>The point in time of the first sample to synthesize,
|
||||
in samples, using the same time scale as in <code>analyze()</code>.</dd>
|
||||
<dt><code>t1</code></dt>
|
||||
<dd>The point in time of the sample one past the last one to synthesize.</dd>
|
||||
<dt><code>signal</code></dt>
|
||||
<dd>The synthesized signal samples will be written here,
|
||||
beginning with the sample from time <code>t0</code> and
|
||||
and ending with the last sample before time <code>t1</code>,
|
||||
for a total of <code>t1 - t0</code> samples.</dd>
|
||||
</dl>
|
||||
<p>The time range <code>t0</code>...<code>t1</code> may extend outside
|
||||
the range analyzed using <code>analyze()</code>, in which case the
|
||||
signal is assumed to be zero in the un-analyzed range.</p>
|
||||
|
||||
<p>A signal may be synthesized in blocks by making multiple calls to
|
||||
<code>analyze()</code> with different sample ranges. For efficiency,
|
||||
the blocks should be large, and each <code>t0</code> should
|
||||
be multiple of a large power of two.<p>
|
||||
|
||||
<h3>Frequency Band Numbering</h3>
|
||||
|
||||
<p>The frequency bands of the analysis filter bank are numbered by
|
||||
nonnegative integers that increase towards lower (sic) frequencies.
|
||||
There is a number of <i>bandpass bands</i> corresponding to the
|
||||
logarithmically spaced bandpass analysis filters, from near 0.5
|
||||
(half the sample rate) to
|
||||
near f<sub>min</sub>, and a single <i>lowpass band</i> containing the
|
||||
residual signal from frequencies below f<sub>min</sub>.
|
||||
The numbering can be examined using the following methods:
|
||||
</p>
|
||||
|
||||
<pre>
|
||||
int bandpass_bands_begin() const;
|
||||
</pre>
|
||||
<p>
|
||||
Return the smallest valid bandpass band number, corresponding to the
|
||||
highest-frequency bandpass filter.</p>
|
||||
<pre>
|
||||
int bandpass_bands_end() const;
|
||||
</pre>
|
||||
<p>
|
||||
Return the bandpass band number one past the highest valid bandpass
|
||||
band number, corresponding to one past the lowest-frequency bandpass
|
||||
filter.
|
||||
</p>
|
||||
<pre>
|
||||
int band_lowpass() const;
|
||||
</pre>
|
||||
<p>
|
||||
Return the band number of the lowpass band.
|
||||
</p>
|
||||
<pre>
|
||||
int band_ref() const;
|
||||
</pre>
|
||||
<p>
|
||||
Return the band number corresponding to the reference frequency
|
||||
<code>ff_ref</code>. If <code>ff_ref</code> falls within
|
||||
the frequency range of the bandpass filter bank, this will
|
||||
be a valid bandpass band number, otherwise it will not.
|
||||
</p>
|
||||
<pre>
|
||||
double band_ff(int band) const;
|
||||
</pre>
|
||||
<p>
|
||||
Return the center frequency of band number <i>band</i>, in units of the
|
||||
sampling frequency.
|
||||
</p>
|
||||
|
||||
<h3>Support</h3>
|
||||
<pre>
|
||||
double analysis_support() const;
|
||||
</pre>
|
||||
<p>Returns the one-sided worst-case time domain <i>support</i> of any of the
|
||||
analysis filters. When calling <code>analyze()</code> with a sample at time <i>t</i>,
|
||||
only spectrogram coefficients within the time range <i>t ± support</i>
|
||||
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.</p>
|
||||
|
||||
<pre>
|
||||
double synthesis_support() const;
|
||||
</pre>
|
||||
<p>Returns the one-sided worst-case time domain <i>support</i> of any of the
|
||||
reconstruction filters. When calling <code>synthesize()</code> to
|
||||
synthesize a sample at time <i>t</i>, the sample will only be
|
||||
significantly affected by spectrogram coefficients in the time
|
||||
range <i>t ± support</i>. Coefficients outside the range may
|
||||
be used in the synthesis, but substituting zeroes for the actual
|
||||
coefficient values will not significantly reduce accuracy.</p>
|
||||
|
||||
</div>
|
||||
<pre>
|
||||
};
|
||||
</pre>
|
||||
|
||||
<h2>Functions</h2>
|
||||
|
||||
<h3>Iterating Over Existing Coefficients</h3>
|
||||
|
||||
<pre>
|
||||
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);
|
||||
</pre>
|
||||
|
||||
<p>
|
||||
Process one or more coefficient sets <code>coefs0</code>... by applying
|
||||
the function <code>f</code> to each coefficient present in <code>coefs0</code>,
|
||||
in an indeterminate order.</p>
|
||||
</p>
|
||||
<p>This can be optionally limited to coefficients whose
|
||||
band number <i>b</i> and sample time <i>t</i> satisfy
|
||||
<code>b0</code> ≤ <i>b</i> < <code>b1</code> and
|
||||
<code>t0</code> ≤ <i>t</i> < <code>t1</code>.
|
||||
To process every coefficient present
|
||||
in <code>coefs0</code>, pass <code>INT_MIN, INT_MAX, INT64_MIN, INT64_MAX</code>
|
||||
for the arguments <code>b0</code>, <code>b1</code>, <code>t0</code>,
|
||||
and <code>t1</code>, respectively.
|
||||
</p>
|
||||
<p>The function <code>f</code> should have the call signature</p>
|
||||
<dd>
|
||||
<pre>
|
||||
template<class T>
|
||||
void f(int b, int64_t t, std::complex<T> &c0, std::complex<T> &ci...);
|
||||
</pre>
|
||||
<p>where</p>
|
||||
<dl>
|
||||
<dt><code>b</code></dt>
|
||||
<dd>The band number of the frequency band the coefficients
|
||||
<code>c0</code> and <code>ci...</code> pertain to.
|
||||
This may be either a bandpass band or the lowpass band.</dd>
|
||||
<dt><code>t</code></dt>
|
||||
<dd>The point in time the coefficients <code>c0</code> and
|
||||
<code>ci...</code> pertain to, in samples</dd>
|
||||
<dt><code>c0</code></dt>
|
||||
<dd>A reference to a complex coefficient from <code>coefs0</code></dd>
|
||||
<dt><code>ci...</code></dt>
|
||||
<dd>Optional references to complex coefficients from the additional
|
||||
coefficient sets <code>coefsi...</code>.</dd>
|
||||
</dl>
|
||||
</dd>
|
||||
</dl>
|
||||
|
||||
<p>The function <code>f</code> may read and/or modify each of the
|
||||
coefficients passed through <code>c0</code> and each
|
||||
<code>ci...</code>.</p>
|
||||
|
||||
<p>The first coefficient set <code>c0</code> is a special case when
|
||||
it comes to the treatment of missing values. Coefficients missing
|
||||
from <code>c0</code> will not be iterated over at all, but when a
|
||||
coefficient <i>is</i> iterated over and is missing from one of the additional
|
||||
coefficient sets <code>ci...</code>, it will be automatically created
|
||||
and initialized to zero in that additional coefficient set.</p>
|
||||
|
||||
<p><i>Note: The template parameters <code>C0</code>
|
||||
and <code>CI</code>... exist to support the processing of coefficient
|
||||
sets containing data of types other
|
||||
than <code>std::complex<T></code>, which is not currently part of the
|
||||
documented API. In typical use, there is no need to specify them when
|
||||
calling <code>apply()</code> because the template parameter list
|
||||
can be deduced, but if they are expicitly specified, they should all
|
||||
be <code>std::complex<T></code>.
|
||||
</i></p>
|
||||
|
||||
<h3>Creating New Coefficients</h3>
|
||||
|
||||
<pre>
|
||||
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);
|
||||
</pre>
|
||||
<p>
|
||||
Fill a region of the time-frequency plane with coefficients
|
||||
and apply the function <code>f</code> to each.
|
||||
</p>
|
||||
<p>This works like <code>process()</code> except that it is not limited
|
||||
to processing coefficients that already exist in <code>coefs0</code>;
|
||||
instead, any missing coefficients in <code>coefs0</code> as well as
|
||||
any of the <code>coefsi</code>... are created and initialized to zero
|
||||
before <code>f</code> is called.</p>
|
||||
|
||||
<p>The <code>t0</code> and <code>t1</code> 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.</p>
|
||||
|
||||
<h3>Forgetting Coefficients</h3>
|
||||
<pre>
|
||||
template <class T>
|
||||
void forget_before(const analyzer<T> &a,
|
||||
coefs<T> &c,
|
||||
int64_t limit);
|
||||
</pre>
|
||||
<p>Allow the coefficients for points in time before <code>limit</code>
|
||||
(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
|
||||
<code>limit</code> are forgotten, only that ones for
|
||||
<code>limit</code> or later are not, and that the amount of memory
|
||||
consumed by any remaining coefficients before <code>limit</code> is
|
||||
bounded.</p>
|
||||
|
||||
<h3>Legacy API For Iterating Over Existing Coefficients</h3>
|
||||
|
||||
<p>Prior to version 1.5, the only way to iterate over
|
||||
coefficients was the <code>apply()</code> function.
|
||||
It is similar to <code>process()</code>, except that it
|
||||
</p>
|
||||
<ul>
|
||||
<li>requires an additional <code>analyzer</code> argument,
|
||||
<li>takes arguments in a different order,
|
||||
<li>applies a function <code>f</code> taking arguments in a different order,
|
||||
<li>does not support restricting the processing to a range of band numbers,
|
||||
<li>only supports iterating over a single coefficient set, and
|
||||
<li>provides default values for t0 and t1.
|
||||
</ul>
|
||||
<p>In new code, <code>process()</code> is preferred.</p>
|
||||
|
||||
<pre>
|
||||
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);
|
||||
</pre>
|
||||
<p>
|
||||
Apply the function <code>f</code> to each coefficient in the coefficient
|
||||
set <code>c</code> for points in time <i>t</i> that satisfy
|
||||
<code>t0</code> ≤ <i>t</i> < <code>t1</code>.
|
||||
If the <code>t0</code> and <code>t1</code> arguments are omitted, <code>f</code>
|
||||
is applied to every coefficient.
|
||||
</p>
|
||||
<dl>
|
||||
<dt><code>a</code></dt>
|
||||
<dd>The spectrum analyzer that produced the coefficients <code>c</code></dd>
|
||||
<dt><code>c</code></dt>
|
||||
<dd>A set of spectrogram coefficients</dd>
|
||||
<dt><code>f</code></dt>
|
||||
<dd>A function to apply to each coefficient in <code>c</code>,
|
||||
with the call signature
|
||||
<pre>
|
||||
template<class T>
|
||||
void f(std::complex<T> &coef, int band, int64_t t);
|
||||
</pre>
|
||||
<dl>
|
||||
<dt><code>coef</code></dt>
|
||||
<dd>A reference to a single complex coefficient. This may be read and/or modified.</dd>
|
||||
<dt><code>band</code></dt>
|
||||
<dd>The band number of the frequency band the coefficient <code>coef0</code> pertains to.
|
||||
This may be either a bandpass band or the lowpass band.</dd>
|
||||
<dt><code>t</code></dt>
|
||||
<dd>The point in time the coefficient <code>c0</code> pertains to, in samples</dd>
|
||||
<dt><code>t0</code></dt><dd>When not <code>INT64_MIN</code>, only apply <code>f</code> to the coefficients for time ≥ <code>t0</code></dd>
|
||||
<dt><code>t1</code></dt><dd>When not <code>INT64_MAX</code>, only apply <code>f</code> to the coefficients for time < <code>t1</code></dd>
|
||||
</dl>
|
||||
</dd>
|
||||
</dl>
|
||||
|
||||
<div class="nav"><span class="prev"><a href="intro.html">Previous: API Introduction</a></span><span class="next"><a href="render_h.html">Next: Spectrogram rendering: <code>render.h</code></a></span></div>
|
||||
|
||||
</body>
|
||||
</html>
|
|
@ -1,41 +0,0 @@
|
|||
<!DOCTYPE html>
|
||||
<!--
|
||||
Copyright (C) 2018-2020 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.
|
||||
-->
|
||||
<html>
|
||||
<head>
|
||||
<link rel="stylesheet" href="../doc.css" type="text/css" />
|
||||
<title>Gaborator reference: API Introdution</title>
|
||||
</head>
|
||||
<body>
|
||||
<h1>Gaborator reference: API Introduction</h1>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<p>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 <code>struct</code>
|
||||
rather than <code>class</code>, 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.
|
||||
</p>
|
||||
|
||||
<p>All definitions are in the namespace <code>gaborator</code>.
|
||||
Applications need to either prefix class names
|
||||
with <code>gaborator::</code>, or use <code>using namespace
|
||||
gaborator;</code>.
|
||||
</p>
|
||||
|
||||
<div class="nav"><span class="prev"><a href="../synth.html">Previous: Example 5: Synthesis from Scratch</a></span><span class="next"><a href="gaborator_h.html">Next: Spectrum analysis and synthesis: <code>gaborator.h</code></a></span></div>
|
||||
|
||||
</body>
|
||||
</html>
|
|
@ -1,74 +0,0 @@
|
|||
<!DOCTYPE html>
|
||||
<!--
|
||||
Copyright (C) 2018-2020 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.
|
||||
-->
|
||||
<html>
|
||||
<head>
|
||||
<link rel="stylesheet" href="../doc.css" type="text/css" />
|
||||
<title>Gaborator reference: render.h</title>
|
||||
</head>
|
||||
<body>
|
||||
<h1>Gaborator reference: <code>render.h</code></h1>
|
||||
|
||||
<h3>Spectrogram Rendering with Power-of-Two Scaling</h3>
|
||||
<pre>
|
||||
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);
|
||||
</pre>
|
||||
<p>Render a rectangular array of pixel values representing signal
|
||||
amplitudes in time-frequency space, optionally scaling up or
|
||||
down by powers of two.
|
||||
</p>
|
||||
<dl>
|
||||
<dt><code>a</code></dt>
|
||||
<dd>The spectrum analyzer that produced the coefficients <code>c</code></dd>
|
||||
<dt><code>c</code></dt>
|
||||
<dd>A set of spectrogram coefficients to render</dd>
|
||||
<dt><code>xorigin</code></dt>
|
||||
<dd>The point in time corresponding to pixel X coordinate 0, in samples</dd>
|
||||
<dt><code>yorigin</code></dt>
|
||||
<dd>The band number of the frequency band corresponding to pixel Y coordinate 0</dd>
|
||||
<dt><code>xi0</code></dt>
|
||||
<dd>The X coordinate of the first pixel to render</dd>
|
||||
<dt><code>xi1</code></dt>
|
||||
<dd>The X coordinate one past the last pixel to render</dd>
|
||||
<dt><code>xe</code></dt>
|
||||
<dd>The horizontal scaling exponent. One horizontal pixel corresponds to 2<sup>xe</sup> signal samples.</dd>
|
||||
<dt><code>yi0</code></dt>
|
||||
<dd>The Y coordinate of the first pixel to render</dd>
|
||||
<dt><code>yi1</code></dt>
|
||||
<dd>The Y coordinate one past the last pixel to render</dd>
|
||||
<dt><code>ye</code></dt>
|
||||
<dd>The vertical scaling exponent. One vertical pixel corresponds to 2<sup>ye</sup> frequency bands.</dd>
|
||||
<dt><code>output</code></dt>
|
||||
<dd>A random access iterator through which the output
|
||||
pixel amplitude values will be written. This is
|
||||
typically a <code>float *</code>. A total of
|
||||
<code>(xi1 - xi0) * (yi1 - yi0))</code> values will be written.
|
||||
</dd>
|
||||
</dl>
|
||||
|
||||
<h3>Utility Functions</h3>
|
||||
<pre>
|
||||
template <class T>
|
||||
unsigned int float2pixel_8bit(T amp);
|
||||
</pre>
|
||||
<p>Convert a normalized amplitude value to a 8-bit greyscale pixel value.</p>
|
||||
<dl>
|
||||
<dt><code>amp</code></dt>
|
||||
<dd>A floating point value representing a signal amplitude, nominally ranging from 0 to 1</dd>
|
||||
</dl>
|
||||
<p>Returns an pixel value ranging from 0 to 255 (inclusive), using an
|
||||
approximation of the sRGB gamma.</p>
|
||||
|
||||
<div class="nav"><span class="prev"><a href="gaborator_h.html">Previous: Spectrum analysis and synthesis: <code>gaborator.h</code></a></span><span class="next"><a href="../overview.html">Next: Overview of Operation</a></span></div>
|
||||
|
||||
</body>
|
||||
</html>
|
|
@ -1,395 +0,0 @@
|
|||
<!DOCTYPE html>
|
||||
<!--
|
||||
Copyright (C) 2017-2020 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.
|
||||
-->
|
||||
<html>
|
||||
<head>
|
||||
<link rel="stylesheet" href="doc.css" />
|
||||
<title>Gaborator Example 1: Rendering a Spectrogram Image</title>
|
||||
</head>
|
||||
<body>
|
||||
<h1>Example 1: Rendering a Spectrogram Image</h1>
|
||||
|
||||
<h2>Introduction</h2>
|
||||
|
||||
<p>This example shows how to generate a greyscale constant-Q
|
||||
spectrogram image from an audio file using the Gaborator library.
|
||||
</p>
|
||||
|
||||
<h2>Preamble</h2>
|
||||
|
||||
<p>We start off with some boilerplate #includes.</p>
|
||||
|
||||
<pre>
|
||||
#include <memory.h>
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <sndfile.h>
|
||||
</pre>
|
||||
|
||||
<p>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
|
||||
<code>gaborator/gaborator.h</code>, and the code for rendering
|
||||
images from the spectrogram coefficients is in
|
||||
<code>gaborator/render.h</code>.</p>
|
||||
|
||||
<pre>
|
||||
#include <gaborator/gaborator.h>
|
||||
#include <gaborator/render.h>
|
||||
</pre>
|
||||
|
||||
<p>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:</p>
|
||||
|
||||
<pre>
|
||||
int main(int argc, char **argv) {
|
||||
if (argc < 3) {
|
||||
std::cerr << "usage: render input.wav output.pgm\n";
|
||||
exit(1);
|
||||
}
|
||||
</pre>
|
||||
|
||||
<h2>Reading the Audio</h2>
|
||||
|
||||
<p>The audio file is read using the <i>libsndfile</i> library
|
||||
and stored in a <code>std::vector<float></code>.
|
||||
Note that although <i>libsndfile</i> is used in this example,
|
||||
the Gaborator library itself does not depend on or
|
||||
use <i>libsndfile</i>.</p>
|
||||
<pre>
|
||||
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);
|
||||
</pre>
|
||||
<p>In case the audio file is a stereo or multi-channel one,
|
||||
mix down the channels to mono, into a new <code>std::vector<float></code>:
|
||||
<pre>
|
||||
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;
|
||||
}
|
||||
</pre>
|
||||
|
||||
<h2>The Spectrum Analysis Parameters</h2>
|
||||
|
||||
<p>Next, we need to choose some parameters for the spectrum analysis:
|
||||
the frequency resolution, the frequency range, and optionally a
|
||||
reference frequency.</p>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<p>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: <code>20.0 / fs</code>.</p>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<p>If desired, one of the frequency bands can be exactly aligned with
|
||||
a <i>reference frequency</i>. When analyzing music signals, this is
|
||||
typically 440 Hz, the standard tuning of the note <i>A<sub>4</sub></i>.
|
||||
Like the minimum frequency, it is given in
|
||||
units of the sample rate, so we pass <code>440.0 / fs</code>.</p>
|
||||
|
||||
<p>The parameters are held in an object of type
|
||||
<code>gaborator::parameters</code>:
|
||||
<pre>
|
||||
gaborator::parameters params(48, 20.0 / fs, 440.0 / fs);
|
||||
</pre>
|
||||
|
||||
<h2>The Spectrum Analyzer</h2>
|
||||
|
||||
<p>Next, we create an object of type <code>gaborator::analyzer</code>;
|
||||
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 <code>float</code>.
|
||||
Constructing the <code>gaborator::analyzer</code> 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 <code>gaborator::analyzer</code>
|
||||
and reusing it is preferable to destroying and recreating it.</p>
|
||||
<pre>
|
||||
gaborator::analyzer<float> analyzer(params);
|
||||
</pre>
|
||||
|
||||
<h2>The Spectrogram Coefficients</h2>
|
||||
|
||||
<p>The result of the spectrum analysis will be a set of <i>spectrogram
|
||||
coefficients</i>. To store them, we will use a <code>gaborator::coefs</code>
|
||||
object. Like the <code>analyzer</code>, 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:</p>
|
||||
<pre>
|
||||
gaborator::coefs<float> coefs(analyzer);
|
||||
</pre>
|
||||
|
||||
<h2>Running the Analysis</h2>
|
||||
|
||||
<p>Now we are ready to do the actual spectrum analysis,
|
||||
by calling the <code>analyze</code> method of the spectrum
|
||||
analyzer object.
|
||||
The first argument to <code>analyze</code> is a <code>float</code> pointer
|
||||
pointing to the first element in the array of samples to analyze.
|
||||
The second and third arguments are of type
|
||||
<code>int64_t</code> 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 <code>mono.size()</code>. The fourth argument is a
|
||||
reference to the set of coefficients that the results of the spectrum
|
||||
analysis will be stored in.
|
||||
</p>
|
||||
<pre>
|
||||
analyzer.analyze(mono.data(), 0, mono.size(), coefs);
|
||||
</pre>
|
||||
|
||||
<h2>Rendering an Image</h2>
|
||||
|
||||
<p>Now there is a set of spectrogram coefficients in <code>coefs</code>.
|
||||
To render them as an image, we will use the function
|
||||
<code>gaborator::render_p2scale()</code>.
|
||||
</p>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<p>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:</p>
|
||||
|
||||
<pre>
|
||||
int64_t x_origin = 0;
|
||||
int64_t y_origin = analyzer.bandpass_bands_begin();
|
||||
</pre>
|
||||
|
||||
<p><code>render_p2scale()</code> 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.
|
||||
|
||||
<p>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
|
||||
2<sup>10</sup> = 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.</p>
|
||||
<pre>
|
||||
int x_scale_exp = 10;
|
||||
</pre>
|
||||
|
||||
<p>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:</p>
|
||||
<pre>
|
||||
while ((n_frames >> x_scale_exp) > 1000)
|
||||
x_scale_exp++;
|
||||
</pre>
|
||||
|
||||
<p>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 2<sup>0</sup>.</p>
|
||||
<pre>
|
||||
int y_scale_exp = 0;
|
||||
</pre>
|
||||
|
||||
<p>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).
|
||||
</p>
|
||||
|
||||
<pre>
|
||||
int64_t x0 = 0;
|
||||
int64_t y0 = 0;
|
||||
</pre>
|
||||
|
||||
<p>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 <code>n_frames</code>,
|
||||
and we get the number of bands as the difference of the end points of
|
||||
the range of band numbers:
|
||||
<code>analyzer.bandpass_bands_end() - analyzer.bandpass_bands_begin()</code>.
|
||||
The scale factor is taken into account by right shifting by the
|
||||
scale exponent.
|
||||
</p>
|
||||
|
||||
<pre>
|
||||
int64_t x1 = n_frames >> x_scale_exp;
|
||||
int64_t y1 = (analyzer.bandpass_bands_end() - analyzer.bandpass_bands_begin()) >> y_scale_exp;
|
||||
</pre>
|
||||
|
||||
<p>The right shift by <code>y_scale_exp</code> above doesn't actually
|
||||
do anything because <code>y_scale_exp</code> is zero, but it would be
|
||||
needed if, for example, you were to change <code>y_scale_exp</code> to
|
||||
1 to get a spectrogram scaled to half the height. You could also make a
|
||||
double-height spectrogram by setting <code>y_scale_exp</code> to -1,
|
||||
but then you also need to change the
|
||||
<code>>> y_scale_exp</code> to
|
||||
<code><< -y_scale_exp</code> since you can't shift by
|
||||
a negative number.
|
||||
</p>
|
||||
|
||||
<p>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
|
||||
<code>(y1 - y0)</code> rows of <code>(x1 - x0)</code> columns
|
||||
each, with the row indices increasing towards lower
|
||||
frequencies and column indices increasing towards later
|
||||
sampling times.
|
||||
</p>
|
||||
<pre>
|
||||
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());
|
||||
</pre>
|
||||
|
||||
<h2>Writing the Image File</h2>
|
||||
|
||||
<p>To keep the code simple and to avoid additional library
|
||||
dependencies, the image is stored in
|
||||
<code>pgm</code> (Portable GreyMap) format, which is simple
|
||||
enough to be generated with just a few lines of inline code.
|
||||
Each amplitude value in <code>amplitudes</code> is converted into an 8-bit
|
||||
gamma corrected pixel value by calling <code>gaborator::float2pixel_8bit()</code>.
|
||||
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.</p>
|
||||
<pre>
|
||||
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();
|
||||
</pre>
|
||||
|
||||
<h2>Postamble</h2>
|
||||
<p>
|
||||
To make the example code a complete program,
|
||||
we just need to finish <code>main()</code>:
|
||||
</p>
|
||||
<pre>
|
||||
return 0;
|
||||
}
|
||||
</pre>
|
||||
|
||||
<a name="compiling"><h2>Compiling</h2></a>
|
||||
<p>
|
||||
If you are using macOS, Linux, NetBSD, or a similar system, you can build
|
||||
the example by running the following command in the <code>examples</code>
|
||||
subdirectory.
|
||||
You need to have <i>libsndfile</i> is installed and supported by
|
||||
<code>pkg-config</code>.
|
||||
</p>
|
||||
<pre class="build Darwin Linux NetBSD FreeBSD">
|
||||
c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` render.cc `pkg-config --libs sndfile` -o render
|
||||
</pre>
|
||||
|
||||
<h2>Compiling for Speed</h2>
|
||||
<p>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
|
||||
<code>GABORATOR_USE_VDSP</code> and linking with the <code>Accelerate</code>
|
||||
framework:
|
||||
</p>
|
||||
<pre class="build Darwin">
|
||||
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
|
||||
</pre>
|
||||
|
||||
<p>On Linux and NetBSD, you can use the PFFFT (Pretty Fast FFT) library.
|
||||
You can get the latest version from
|
||||
<a href="https://bitbucket.org/jpommier/pffft">https://bitbucket.org/jpommier/pffft</a>,
|
||||
or the exact version that was used for testing from gaborator.com:
|
||||
</p>
|
||||
<!-- ftp https://bitbucket.org/jpommier/pffft/get/29e4f76ac53b.zip -->
|
||||
<pre class="build Linux NetBSD FreeBSD">
|
||||
wget http://download.gaborator.com/mirror/pffft/29e4f76ac53b.zip
|
||||
unzip 29e4f76ac53b.zip
|
||||
mv jpommier-pffft-29e4f76ac53b pffft
|
||||
</pre>
|
||||
<p>Then, compile it:</p>
|
||||
<pre class="build Linux NetBSD FreeBSD">
|
||||
cc -c -O3 -ffast-math pffft/pffft.c -o pffft/pffft.o
|
||||
</pre>
|
||||
<p>(If you are building for ARM, you will need to add <code>-mfpu=neon</code> to
|
||||
both the above compilation command and the ones below.)</p>
|
||||
<p>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:</p>
|
||||
<pre class="build Linux NetBSD FreeBSD">
|
||||
cc -c -O3 -ffast-math -DFFTPACK_DOUBLE_PRECISION pffft/fftpack.c -o pffft/fftpack.o
|
||||
</pre>
|
||||
<p>Then build the example and link it with both PFFFT and FFTPACK:</p>
|
||||
<pre class="build Linux NetBSD FreeBSD">
|
||||
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
|
||||
</pre>
|
||||
|
||||
<h2>Running</h2>
|
||||
<p>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 <code>.pgm</code> image, and then convert
|
||||
the <code>.pgm</code> image into a <code>JPEG</code> image:
|
||||
<pre class="run">
|
||||
wget http://download.gaborator.com/audio/guitar.wav
|
||||
./render guitar.wav guitar.pgm
|
||||
cjpeg <guitar.pgm >guitar.jpg
|
||||
</pre>
|
||||
|
||||
<h2>Example Output</h2>
|
||||
<p>The JPEG file produced by the above will look like this:</p>
|
||||
<img src="spectrogram.jpg" alt="Spectrogram" data-autogen="no">
|
||||
|
||||
<div class="nav"><span class="next"><a href="filter.html">Next: Example 2: Frequency-Domain Filtering</a></span></div>
|
||||
|
||||
</body>
|
||||
</html>
|
|
@ -1,121 +0,0 @@
|
|||
<!DOCTYPE html>
|
||||
<!--
|
||||
Copyright (C) 2019-2020 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.
|
||||
-->
|
||||
<html>
|
||||
<head>
|
||||
<link rel="stylesheet" href="doc.css" />
|
||||
<title>Gaborator Example 4: Measuring the Signal-to-Noise Ratio</title>
|
||||
</head>
|
||||
<body>
|
||||
<h1>Example 4: Measuring the Signal-to-Noise Ratio</h1>
|
||||
|
||||
<h2>Introduction</h2>
|
||||
|
||||
<p>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.
|
||||
</p>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<h2>Preamble</h2>
|
||||
|
||||
<pre>
|
||||
#include <iostream>
|
||||
#include <iomanip>
|
||||
#include <random>
|
||||
#include <gaborator/gaborator.h>
|
||||
</pre>
|
||||
|
||||
<h2>Amplitude Measurement</h2>
|
||||
<p>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 <code>rms()</code>.
|
||||
</p>
|
||||
<pre>
|
||||
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);
|
||||
}
|
||||
</pre>
|
||||
|
||||
<h2>Main Program</h2>
|
||||
|
||||
<p>For the test signal, we use a million samples of white noise with a
|
||||
uniform amplitude distribution between -1 and +1.</p>
|
||||
<pre>
|
||||
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);
|
||||
</pre>
|
||||
<p>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):</p>
|
||||
<pre>
|
||||
gaborator::parameters params(48, 5e-4);
|
||||
gaborator::analyzer<float> analyzer(params);
|
||||
</pre>
|
||||
<p>...and run the spectrum analyzis:</p>
|
||||
<pre>
|
||||
gaborator::coefs<float> coefs(analyzer);
|
||||
analyzer.analyze(signal_in.data(), 0, len, coefs);
|
||||
</pre>
|
||||
<p>...resynthesize the signal into <code>signal_out</code>:
|
||||
<pre>
|
||||
std::vector<float> signal_out(len);
|
||||
analyzer.synthesize(coefs, 0, len, signal_out.data());
|
||||
</pre>
|
||||
<p>...measure the resynthesis error:</p>
|
||||
<pre>
|
||||
std::vector<float> error(len);
|
||||
for (size_t i = 0; i < len; i++)
|
||||
error[i] = signal_out[i] - signal_in[i];
|
||||
</pre>
|
||||
<p>...calculate the signal-to-noise ratio:</p>
|
||||
<pre>
|
||||
double snr = rms(signal_in) / rms(error);
|
||||
</pre>
|
||||
<p>...and print it in decibels:</p>
|
||||
<pre>
|
||||
std::cout << std::fixed << std::setprecision(1) << 20 * log10(snr) << " dB\n";
|
||||
}
|
||||
</pre>
|
||||
<h2>Compiling</h2>
|
||||
<p>Like <a href="render.html#compiling">Example 1</a>, this example
|
||||
can be built using a one-line build command:
|
||||
</p>
|
||||
<pre class="build Darwin Linux NetBSD FreeBSD">
|
||||
c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` snr.cc `pkg-config --libs sndfile` -o snr
|
||||
</pre>
|
||||
<p>Or using the vDSP FFT on macOS:</p>
|
||||
<pre class="build Darwin">
|
||||
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
|
||||
</pre>
|
||||
<p>Or using PFFFT (see <a href="render.html#compiling">Example 1</a> for how to download and build PFFFT):</p>
|
||||
<pre class="build">
|
||||
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
|
||||
</pre>
|
||||
|
||||
<h2>Running</h2>
|
||||
<p>The program is run with no arguments:</p>
|
||||
<pre class="run">
|
||||
./snr
|
||||
</pre>
|
||||
<p>This will print the SNR which should be more than 100 dB if the library is working correctly.</p>
|
||||
|
||||
<div class="nav"><span class="prev"><a href="stream.html">Previous: Example 3: Streaming</a></span><span class="next"><a href="synth.html">Next: Example 5: Synthesis from Scratch</a></span></div>
|
||||
|
||||
</body>
|
||||
</html>
|
Binary file not shown.
Before Width: | Height: | Size: 22 KiB |
|
@ -1,279 +0,0 @@
|
|||
<!DOCTYPE html>
|
||||
<!--
|
||||
Copyright (C) 2018-2020 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.
|
||||
-->
|
||||
<html>
|
||||
<head>
|
||||
<link rel="stylesheet" href="doc.css" />
|
||||
<title>Gaborator Example 3: Streaming</title>
|
||||
</head>
|
||||
<body>
|
||||
<h1>Example 3: Streaming</h1>
|
||||
|
||||
<h2>Introduction</h2>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<p>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 <code>coef =
|
||||
-coef</code> a placeholder for your own streaming filter or effect
|
||||
code.</p>
|
||||
|
||||
<h2>Preamble</h2>
|
||||
|
||||
<pre>
|
||||
#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);
|
||||
}
|
||||
</pre>
|
||||
|
||||
<h2>Opening the Streams</h2>
|
||||
|
||||
<p>We again use <i>libsndfile</i> to read the input and write the output.
|
||||
To keep it simple, this example only handles mono files.</p>
|
||||
<pre>
|
||||
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);
|
||||
}
|
||||
</pre>
|
||||
<p>The next couple of lines work around a design flaw in
|
||||
<i>libsndfile</i>. By default, when reading a 16-bit
|
||||
audio file as floating point data and then writing them
|
||||
as another 16-bit audio file, <i>libsndfile</i> 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.</p>
|
||||
<pre>
|
||||
sf_command(sf_in, SFC_SET_NORM_FLOAT, NULL, SF_FALSE);
|
||||
sf_command(sf_out, SFC_SET_NORM_FLOAT, NULL, SF_FALSE);
|
||||
</pre>
|
||||
|
||||
<h2>Spectrum Analysis Parameters</h2>
|
||||
|
||||
<p>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.</p>
|
||||
<pre>
|
||||
gaborator::parameters params(12, 200.0 / fs, 440.0 / fs);
|
||||
gaborator::analyzer<float> analyzer(params);
|
||||
</pre>
|
||||
|
||||
<h2>Calculating Latency</h2>
|
||||
<p>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 <i>t</i> is a weighted average of samples
|
||||
from both before and after <i>t</i>, 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 <i>support</i>, and the worst-case (widest) support of
|
||||
any analysis filter can be found by calling the function
|
||||
<code>gaborator::analyzer::analysis_support()</code>. This returns
|
||||
the <i>one-sided</i> support, the width of the impulse
|
||||
response <i>to each side</i> of its center, as a floating point number.
|
||||
To be on the safe side, let's round this up to the next integer:</p>
|
||||
<pre>
|
||||
size_t analysis_support = ceil(analyzer.analysis_support());
|
||||
</pre>
|
||||
<p>Similarly, when resynthesizing audio from coefficients, calculating
|
||||
a sample at time <i>t</i> 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
|
||||
<code>gaborator::analyzer::synthesis_support()</code>:
|
||||
</p>
|
||||
<pre>
|
||||
size_t synthesis_support = ceil(analyzer.synthesis_support());
|
||||
</pre>
|
||||
|
||||
<p>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 <i>t</i> can only do so when it has received the input samples up
|
||||
to time <i>t + analysis_support</i>, and therefore has a minimum latency of
|
||||
<i>analysis_support</i>. 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 <i>t</i> once
|
||||
it has received the input samples up to
|
||||
<i>t + analysis_support + synthesis_support</i>,
|
||||
for a minimum latency of <i>analysis_support + synthesis_support</i>.
|
||||
Let's print this total latency to standard output:
|
||||
</p>
|
||||
<pre>
|
||||
std::cerr << "latency: " << analysis_support + synthesis_support << " samples\n";
|
||||
</pre>
|
||||
|
||||
<p>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.
|
||||
</p>
|
||||
|
||||
<h2>Streaming</h2>
|
||||
<p>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 <code>t_in</code> keeps track of time, indicating
|
||||
the sampling time of the first sample of the current input block,
|
||||
in units of samples.
|
||||
</p>
|
||||
<pre>
|
||||
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);
|
||||
</pre>
|
||||
<p>Now we can spectrum analyze the current block of samples. Note how
|
||||
the time range,
|
||||
<code>t_in</code>...<code>t_in + blocksize</code>,
|
||||
is explicitly passed to <code>analyze()</code>.
|
||||
</p>
|
||||
<pre>
|
||||
analyzer.analyze(buf.data(), t_in, t_in + blocksize, coefs);
|
||||
</pre>
|
||||
<p>The call to <code>analyze()</code> updates the coefficients
|
||||
for the time range from <code>t_in - analysis_support</code> to
|
||||
<code>t_in + blocksize + analysis_support</code>. The oldest
|
||||
<code>blocksize</code> samples of this time range,
|
||||
that is, from <code>t_in - analysis_support</code> to
|
||||
<code>t_in - analysis_support + blocksize</code>, 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 <code>prorcess()</code>
|
||||
applied a function to all the coefficients, here it is applied only to
|
||||
the coefficients within a limited time range.
|
||||
</p>
|
||||
<pre>
|
||||
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);
|
||||
</pre>
|
||||
|
||||
<p>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 <code>t</code> will depend on the coefficients of the
|
||||
time range <code>t - synthesis_support</code>...<code>t +
|
||||
synthesis_support</code>. To ensure that the resynthesis uses only
|
||||
coefficients that have already been processed by
|
||||
the <code>process()</code> call above, the most recent block of samples
|
||||
that can safely be resynthesized ranges from <code>t_out = t_in -
|
||||
analysis_support - synthesis_support</code> to <code>t_out +
|
||||
blocksize</code>.</p>
|
||||
<pre>
|
||||
int64_t t_out = t_in - analysis_support - synthesis_support;
|
||||
analyzer.synthesize(coefs, t_out, t_out + blocksize, buf.data());
|
||||
</pre>
|
||||
<p>The synthesized audio can now be written to the output file:</p>
|
||||
<pre>
|
||||
sf_count_t n_written = sf_writef_float(sf_out, buf.data(), blocksize);
|
||||
if (n_written != blocksize) {
|
||||
std::cerr << "write error\n";
|
||||
exit(1);
|
||||
}
|
||||
</pre>
|
||||
<p>Coefficients older than <code>t_out + blocksize - synthesis_support</code>
|
||||
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:
|
||||
</p>
|
||||
<pre>
|
||||
forget_before(analyzer, coefs, t_out + blocksize - synthesis_support);
|
||||
</pre>
|
||||
<p>This concludes the block-by-block processing loop.</p>
|
||||
<pre>
|
||||
t_in += blocksize;
|
||||
}
|
||||
</pre>
|
||||
|
||||
<h2>Postamble</h2>
|
||||
<pre>
|
||||
sf_close(sf_in);
|
||||
sf_close(sf_out);
|
||||
return 0;
|
||||
}
|
||||
</pre>
|
||||
|
||||
<h2>Compiling</h2>
|
||||
<p>Like the previous ones, this example can also be built using a one-line build command:
|
||||
</p>
|
||||
<pre class="build Darwin Linux NetBSD FreeBSD">
|
||||
c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` stream.cc `pkg-config --libs sndfile` -o stream
|
||||
</pre>
|
||||
<p>Or using the vDSP FFT on macOS:</p>
|
||||
<pre class="build Darwin">
|
||||
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
|
||||
</pre>
|
||||
<p>Or using PFFFT (see <a href="render.html#compiling">Example 1</a> for how to download and build PFFFT):</p>
|
||||
<pre class="build">
|
||||
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
|
||||
</pre>
|
||||
|
||||
<h2>Running</h2>
|
||||
<p>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.</p>
|
||||
<pre class="run">
|
||||
wget http://download.gaborator.com/audio/impulse.wav
|
||||
./stream impulse.wav impulse_streamed.wav
|
||||
</pre>
|
||||
|
||||
<p>The file <code>impulse_streamed.wav</code> will be identical to
|
||||
<code>impulse.wav</code> except that the impulse will be of
|
||||
opposite polarity, and delayed by the latency of
|
||||
<code>analysis_support + synthesis_support</code> samples.</p>
|
||||
|
||||
<div class="nav"><span class="prev"><a href="filter.html">Previous: Example 2: Frequency-Domain Filtering</a></span><span class="next"><a href="snr.html">Next: Example 4: Measuring the Signal-to-Noise Ratio</a></span></div>
|
||||
|
||||
</body>
|
||||
</html>
|
|
@ -1,213 +0,0 @@
|
|||
<!DOCTYPE html>
|
||||
<!--
|
||||
Copyright (C) 2020 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.
|
||||
-->
|
||||
<html>
|
||||
<head>
|
||||
<link rel="stylesheet" href="doc.css" />
|
||||
<title>Gaborator Example 5: Synthesis from Scratch</title>
|
||||
</head>
|
||||
<body>
|
||||
<h1>Example 5: Synthesis from Scratch</h1>
|
||||
|
||||
<h2>Introduction</h2>
|
||||
|
||||
<p>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.
|
||||
</p>
|
||||
|
||||
<h2>Preamble</h2>
|
||||
|
||||
<p>This example program takes a single command line argument, the name
|
||||
of the output file.</p>
|
||||
<pre>
|
||||
#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);
|
||||
}
|
||||
</pre>
|
||||
|
||||
<h2>Synthesis Parameters</h2>
|
||||
|
||||
<p>Although this example does not perform any analysis, we nonetheless
|
||||
need to create an <code>analyzer</code> 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.</p>
|
||||
|
||||
<p>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.</p>
|
||||
|
||||
<pre>
|
||||
double fs = 44100;
|
||||
gaborator::parameters params(12, 20.0 / fs, 8.18 / fs);
|
||||
gaborator::analyzer<float> analyzer(params);
|
||||
</pre>
|
||||
|
||||
<h2>Melody Parameters</h2>
|
||||
|
||||
<p>
|
||||
We will use the A minor pentatonic scale, which contains the
|
||||
following notes (using the MIDI note numbering):</p>
|
||||
<pre>
|
||||
static int pentatonic[] = { 57, 60, 62, 64, 67 };
|
||||
</pre>
|
||||
|
||||
<p>
|
||||
The melody will consist of 64 notes, at a tempo of 120 beats per
|
||||
minute:
|
||||
</p>
|
||||
<pre>
|
||||
int n_notes = 64;
|
||||
double tempo = 120.0;
|
||||
double beat_duration = 60.0 / tempo;
|
||||
</pre>
|
||||
|
||||
<p>
|
||||
The variable <code>volume</code> determines the amplitude of
|
||||
each note, and has been chosen such that there will be no clipping
|
||||
of the final output.
|
||||
</p>
|
||||
<pre>
|
||||
float volume = 0.2;
|
||||
</pre>
|
||||
|
||||
<h2>Composition</h2>
|
||||
|
||||
<p>We start with an empty coefficient set:</p>
|
||||
<pre>
|
||||
gaborator::coefs<float> coefs(analyzer);
|
||||
</pre>
|
||||
|
||||
<p>Each note is chosen randomly from the pentatonic scale and added
|
||||
to the coefficient set by calling the function <code>fill()</code>.
|
||||
The <code>fill()</code> function is similar to the <code>process()</code>
|
||||
function used in previous examples, except that it can be used to
|
||||
create new coefficients rather than just modifying existing ones.</p>
|
||||
|
||||
<p>Each note is created by calling <code>fill()</code> 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 <code>+=</code> operator
|
||||
rather than overwriting them.</p>
|
||||
|
||||
<p>Note that band numbers increase towards lower frequencies but MIDI
|
||||
note numbers increase towards higher frequencies, hence the minus sign
|
||||
in front of <code>midi_note</code>.
|
||||
</p>
|
||||
|
||||
<pre>
|
||||
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);
|
||||
}
|
||||
</pre>
|
||||
|
||||
<h2>Synthesis</h2>
|
||||
|
||||
<p>We can now synthesize audio from the coefficients by
|
||||
calling <code>synthesize()</code>. 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.
|
||||
</p>
|
||||
<pre>
|
||||
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());
|
||||
</pre>
|
||||
|
||||
<h2>Writing the Audio</h2>
|
||||
|
||||
<p>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
|
||||
<code>sfinfo</code> structure:</p>
|
||||
<pre>
|
||||
SF_INFO sfinfo;
|
||||
memset(&sfinfo, 0, sizeof(sfinfo));
|
||||
sfinfo.samplerate = fs;
|
||||
sfinfo.channels = 1;
|
||||
sfinfo.format = SF_FORMAT_WAV | SF_FORMAT_PCM_16;
|
||||
</pre>
|
||||
|
||||
<p>The rest is identical to
|
||||
<a href="filter.html#writing_audio_code">Example 2</a>:
|
||||
</p>
|
||||
<pre>
|
||||
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;
|
||||
}
|
||||
</pre>
|
||||
|
||||
<h2>Compiling</h2>
|
||||
<p>Like <a href="render.html#compiling">Example 1</a>, this example
|
||||
can be built using a one-line build command:
|
||||
</p>
|
||||
<pre class="build Darwin Linux NetBSD FreeBSD">
|
||||
c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` synth.cc `pkg-config --libs sndfile` -o synth
|
||||
</pre>
|
||||
<p>Or using the vDSP FFT on macOS:</p>
|
||||
<pre class="build Darwin">
|
||||
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
|
||||
</pre>
|
||||
<p>Or using PFFFT (see <a href="render.html">Example 1</a> for how to download and build PFFFT):</p>
|
||||
<pre class="build">
|
||||
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
|
||||
</pre>
|
||||
|
||||
<h2>Running</h2>
|
||||
<p>The example program can be run using the command</p>
|
||||
<pre class="run">
|
||||
./synth melody.wav
|
||||
</pre>
|
||||
<p>The resulting audio will be in <code>melody.wav</code>.</p>
|
||||
|
||||
<div class="nav"><span class="prev"><a href="snr.html">Previous: Example 4: Measuring the Signal-to-Noise Ratio</a></span><span class="next"><a href="ref/intro.html">Next: API Introduction</a></span></div>
|
||||
|
||||
</body>
|
||||
</html>
|
|
@ -1,69 +0,0 @@
|
|||
// See ../doc/filter.html for commentary
|
||||
|
||||
#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);
|
||||
}
|
||||
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);
|
||||
gaborator::parameters params(100, 20.0 / fs);
|
||||
gaborator::analyzer<float> analyzer(params);
|
||||
std::vector<float> 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<float> channel(n_frames);
|
||||
for (sf_count_t i = 0; i < n_frames; i++)
|
||||
channel[i] = audio[i * sfinfo.channels + ch];
|
||||
gaborator::coefs<float> coefs(analyzer);
|
||||
analyzer.analyze(channel.data(), 0, channel.size(), coefs);
|
||||
process([&](int band, int64_t, std::complex<float> &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;
|
||||
}
|
|
@ -1,69 +0,0 @@
|
|||
// See ../doc/render.html for commentary
|
||||
|
||||
#include <memory.h>
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <sndfile.h>
|
||||
#include <gaborator/gaborator.h>
|
||||
#include <gaborator/render.h>
|
||||
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<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);
|
||||
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;
|
||||
}
|
||||
gaborator::parameters params(48, 20.0 / fs, 440.0 / fs);
|
||||
gaborator::analyzer<float> analyzer(params);
|
||||
gaborator::coefs<float> 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<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());
|
||||
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;
|
||||
}
|
|
@ -1,32 +0,0 @@
|
|||
// See ../doc/snr.html for commentary
|
||||
|
||||
#include <iostream>
|
||||
#include <iomanip>
|
||||
#include <random>
|
||||
#include <gaborator/gaborator.h>
|
||||
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);
|
||||
}
|
||||
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);
|
||||
gaborator::parameters params(48, 5e-4);
|
||||
gaborator::analyzer<float> analyzer(params);
|
||||
gaborator::coefs<float> coefs(analyzer);
|
||||
analyzer.analyze(signal_in.data(), 0, len, coefs);
|
||||
std::vector<float> signal_out(len);
|
||||
analyzer.synthesize(coefs, 0, len, signal_out.data());
|
||||
std::vector<float> 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";
|
||||
}
|
|
@ -1,72 +0,0 @@
|
|||
// See ../doc/stream.html for commentary
|
||||
|
||||
#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);
|
||||
}
|
||||
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<float> 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<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);
|
||||
analyzer.analyze(buf.data(), t_in, t_in + blocksize, coefs);
|
||||
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);
|
||||
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;
|
||||
}
|
|
@ -1,62 +0,0 @@
|
|||
// See ../doc/synth.html for commentary
|
||||
|
||||
#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);
|
||||
}
|
||||
double fs = 44100;
|
||||
gaborator::parameters params(12, 20.0 / fs, 8.18 / fs);
|
||||
gaborator::analyzer<float> 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<float> 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<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);
|
||||
}
|
||||
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());
|
||||
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;
|
||||
}
|
|
@ -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
|
|
@ -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
|
|
@ -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 <algorithm>
|
||||
#include <complex>
|
||||
#include <iterator>
|
||||
#include <vector>
|
||||
|
||||
#include <memory.h>
|
||||
|
||||
namespace gaborator {
|
||||
|
||||
template <class I>
|
||||
struct fft {
|
||||
typedef typename std::iterator_traits<I>::value_type C; // complex
|
||||
typedef typename C::value_type T; // float/double
|
||||
typedef typename std::vector<C> 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 <class CI>
|
||||
struct rfft {
|
||||
typedef typename std::iterator_traits<CI>::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<CI> cf;
|
||||
};
|
||||
|
||||
} // Namespace
|
||||
|
||||
#endif
|
|
@ -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 <assert.h>
|
||||
#include <complex>
|
||||
|
||||
#include <iterator>
|
||||
#include <vector>
|
||||
|
||||
|
||||
#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<std::complex<float> *> {
|
||||
typedef std::complex<float> *I;
|
||||
typedef const std::complex<float> *CONST_I;
|
||||
typedef std::iterator_traits<I>::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<std::complex<float>>::iterator
|
||||
|
||||
template <>
|
||||
struct fft<std::vector<std::complex<float>>::iterator>:
|
||||
public fft<std::complex<float> *>
|
||||
{
|
||||
typedef fft<std::complex<float> *> base;
|
||||
typedef std::vector<std::complex<float>>::iterator I;
|
||||
fft(unsigned int n_): fft<std::complex<float> *>(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<std::complex<double> *> {
|
||||
typedef std::complex<double> *I;
|
||||
typedef const std::complex<double> *CONST_I;
|
||||
typedef std::iterator_traits<I>::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<double> wsave;
|
||||
};
|
||||
|
||||
// Real FFT
|
||||
|
||||
template <>
|
||||
struct rfft<std::complex<float> *> {
|
||||
typedef std::complex<float> *CI; // Complex iterator
|
||||
typedef const std::complex<float> *CONST_CI;
|
||||
typedef typename std::iterator_traits<CI>::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<CI>(in)[0] = C(tmp.real(), in[n >> 1].real());
|
||||
pffft_transform_ordered(setup, (const float *) in, out, NULL, PFFFT_BACKWARD);
|
||||
const_cast<CI>(in)[0] = tmp;
|
||||
}
|
||||
|
||||
private:
|
||||
// Size of the transform
|
||||
unsigned int n;
|
||||
PFFFT_Setup *setup;
|
||||
};
|
||||
|
||||
#undef GABORATOR_PFFFT_CHECK_ALIGN
|
||||
|
||||
} // namespace
|
||||
|
||||
#endif
|
|
@ -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 <assert.h>
|
||||
|
||||
#include <iterator>
|
||||
#include <vector>
|
||||
|
||||
#include <memory.h>
|
||||
|
||||
#include <mach/mach.h>
|
||||
#include <mach/task.h>
|
||||
#include <mach/task_info.h>
|
||||
#include <mach/vm_map.h>
|
||||
|
||||
#include <Accelerate/Accelerate.h>
|
||||
|
||||
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<std::complex<float> *> {
|
||||
typedef std::complex<float> *I;
|
||||
typedef typename std::iterator_traits<I>::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<std::complex<float> *> {
|
||||
typedef std::complex<float> *CI; // Complex iterator
|
||||
typedef const std::complex<float> *CONST_CI;
|
||||
typedef typename std::iterator_traits<CI>::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<CI>(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<CI>(in)[0] = tmp;
|
||||
}
|
||||
|
||||
private:
|
||||
// Size of the transform
|
||||
unsigned int n;
|
||||
unsigned int log2n;
|
||||
FFTSetup setup;
|
||||
};
|
||||
|
||||
// Support transforming std::vector<std::complex<float>>::iterator
|
||||
|
||||
template <>
|
||||
struct fft<typename std::vector<std::complex<float>>::iterator>:
|
||||
public fft<std::complex<float> *>
|
||||
{
|
||||
typedef fft<std::complex<float> *> base;
|
||||
typedef typename std::vector<std::complex<float>>::iterator I;
|
||||
fft(unsigned int n_): fft<std::complex<float> *>(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
|
File diff suppressed because it is too large
Load diff
|
@ -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 <math.h>
|
||||
|
||||
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
|
|
@ -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 <stdlib.h> // size_t
|
||||
|
||||
#include <algorithm> // 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 <class T>
|
||||
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<T *>(::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<T *>(::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<T *>(::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
|
|
@ -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 <map>
|
||||
|
||||
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 <class T, class K>
|
||||
struct pool {
|
||||
typedef std::map<K, T *> 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<typename m_t::iterator, bool> 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 <class T, class K>
|
||||
pool<T, K> pool<T, K>::shared;
|
||||
|
||||
} // namespace
|
||||
|
||||
#endif
|
|
@ -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 <class T> 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 <class T>
|
||||
void incref(T &r) {
|
||||
r.refcount++;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
void decref(T &r) {
|
||||
r.refcount--;
|
||||
if (r.refcount == 0)
|
||||
delete &r;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
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
|
|
@ -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 <class T>
|
||||
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<class T>
|
||||
struct complex_abs_fob {
|
||||
T operator()(const complex<T> &c) {
|
||||
return complex_abs(c);
|
||||
}
|
||||
typedef T return_type;
|
||||
};
|
||||
|
||||
|
||||
// T -> f() -> OI::value_type
|
||||
|
||||
template <class F, class OI, class T>
|
||||
struct transform_output_iterator: public std::iterator<std::output_iterator_tag, T> {
|
||||
typedef T value_type;
|
||||
transform_output_iterator(F f_, OI output_): f(f_), output(output_) { }
|
||||
transform_output_iterator<F, OI, T>& operator=(T v) {
|
||||
*output++ = f(v);
|
||||
return *this;
|
||||
}
|
||||
transform_output_iterator<F, OI, T>& operator*() { return *this; }
|
||||
transform_output_iterator<F, OI, T>& operator++() { return *this; }
|
||||
transform_output_iterator<F, OI, T>& 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 <class T, class C, class OI, class NORMF>
|
||||
struct abs_row_source {
|
||||
typedef transform_output_iterator<NORMF, OI, C> abs_writer_t;
|
||||
abs_row_source(const coefs<T, C> &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<T, abs_writer_t, C> rs;
|
||||
NORMF normf;
|
||||
};
|
||||
|
||||
// Helper class for abs_row_source specialization below
|
||||
|
||||
template <class C, class OI>
|
||||
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 <class T, class C, class OI>
|
||||
struct abs_row_source<T, C, OI, struct complex_abs_fob<T>> {
|
||||
// Note unused last arg
|
||||
abs_row_source(const coefs<T, C> &coefs_,
|
||||
int oct_, unsigned int obno_,
|
||||
complex_abs_fob<T>):
|
||||
slicer(coefs_, oct_, obno_)
|
||||
{ }
|
||||
OI operator()(coef_index_t i0, coef_index_t i1, OI output) const {
|
||||
abs_writer_dest<C, OI> dest(output);
|
||||
slicer(i0, i1, dest);
|
||||
return dest.output;
|
||||
}
|
||||
row_foreach_slice<T, abs_writer_dest<C, OI>, C> slicer;
|
||||
};
|
||||
|
||||
// Render a single line (single frequency band), with scaling in the
|
||||
// horizontal (time) dimension, and filtering to avoid aliasing when
|
||||
// minifying.
|
||||
|
||||
template <class OI, class T, class C, class NORMF, class RESAMPLER>
|
||||
OI render_line(const analyzer<T> &anl,
|
||||
const coefs<T, C> &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<T, C, RST *, NORMF>
|
||||
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 <class OI, class T, class C, class NORMF, class RESAMPLER>
|
||||
OI render_noyscale(const analyzer<T> &anl,
|
||||
const coefs<T, C> &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<OI, T, C, NORMF, RESAMPLER>
|
||||
(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 <class T, class OI>
|
||||
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 <class I>
|
||||
struct stride_iterator: public std::iterator
|
||||
<std::forward_iterator_tag, typename std::iterator_traits<I>::value_type>
|
||||
{
|
||||
typedef typename std::iterator_traits<I>::value_type T;
|
||||
stride_iterator(I it_, size_t stride_): it(it_), stride(stride_) { }
|
||||
T& operator*() { return *it; }
|
||||
stride_iterator<I>& 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<T>()"
|
||||
// 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 OI, class T, class C, class NORMF = complex_abs_fob<T>,
|
||||
class RESAMPLER = lanczos2_pow2_resampler,
|
||||
class UPDATEDF = updated_nop>
|
||||
void render_incremental(
|
||||
const analyzer<T> &anl,
|
||||
const coefs<T, C> &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<T>(),
|
||||
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<OI, T, C, NORMF, RESAMPLER, UPDATEDF>
|
||||
(anl, msc, x_xf, y_xf,
|
||||
xi0, xi1,
|
||||
yi0, ysplit,
|
||||
inc_i0, inc_i1,
|
||||
output, output_stride, normf, updated);
|
||||
render_incremental<OI, T, C, NORMF, RESAMPLER, UPDATEDF>
|
||||
(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<OI, T, C, NORMF, RESAMPLER>
|
||||
(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<RST> render_data(n_pixels);
|
||||
|
||||
// Render data resampled in the X direction
|
||||
RST *p = render_data.data();
|
||||
render_noyscale<OI, T, C, NORMF, RESAMPLER>
|
||||
(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<RST, OI> src(render_data.data(),
|
||||
xi0, xi1, ysi0, ysi1,
|
||||
xi);
|
||||
stride_iterator<OI> dest(output + (xi - xi0), output_stride);
|
||||
y_resampler(src, yi0, yi1, dest);
|
||||
}
|
||||
}
|
||||
updated(xi0, xi1, yi0, yi1);
|
||||
}
|
||||
|
||||
template <class OI, class T, class C, class NORMF = complex_abs_fob<T>,
|
||||
class RESAMPLER = lanczos2_pow2_resampler>
|
||||
void render_p2scale(const analyzer<T> &anl,
|
||||
const coefs<T, C> &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<T>())
|
||||
{
|
||||
// Provide default inc_i0, inc_i1, and output_stride
|
||||
render_incremental<OI, T, C, NORMF, RESAMPLER>
|
||||
(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
|
|
@ -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 <assert.h>
|
||||
#include <inttypes.h>
|
||||
#include <math.h>
|
||||
#include <algorithm> // 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 <class S, class OI>
|
||||
OI resample2_p2xf(const S &source, p2_transform xf,
|
||||
int64_t i0, int64_t i1,
|
||||
bool interpolate, OI out)
|
||||
{
|
||||
typedef typename std::iterator_traits<OI>::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<T> 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<T> 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<T> 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 <class S, class OI>
|
||||
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<OI>::value_type T;
|
||||
gaborator::pod_vector<T> 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 <class S, class OI>
|
||||
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 <class S, class OI>
|
||||
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
|
|
@ -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 <assert.h>
|
||||
|
||||
#if GABORATOR_USE_SSE3_INTRINSICS
|
||||
#include <pmmintrin.h>
|
||||
#endif
|
||||
|
||||
#include <complex>
|
||||
|
||||
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 <class T>
|
||||
std::complex<T> complex_mul_naive(std::complex<T> a,
|
||||
std::complex<T> b)
|
||||
{
|
||||
return std::complex<T>(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<double> complex_mul(std::complex<double> a_,
|
||||
std::complex<double> 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<double>(c[0], c[1]);
|
||||
}
|
||||
#else
|
||||
static inline
|
||||
std::complex<double> complex_mul(std::complex<double> a_,
|
||||
std::complex<double> b_)
|
||||
{
|
||||
return complex_mul_naive(a_, b_);
|
||||
}
|
||||
#endif
|
||||
|
||||
static inline
|
||||
std::complex<float> complex_mul(std::complex<float> a_,
|
||||
std::complex<float> b_)
|
||||
{
|
||||
return complex_mul_naive(a_, b_);
|
||||
}
|
||||
|
||||
template <class T, class U, class V>
|
||||
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 <class T, class U, class V, class S>
|
||||
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 <class I, class O>
|
||||
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 <pmmintrin.h>
|
||||
|
||||
// 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<float> *cv,
|
||||
const std::complex<float> *av,
|
||||
const std::complex<float> *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<float> *cv,
|
||||
const std::complex<float> *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<float> *cv,
|
||||
const std::complex<float> *av,
|
||||
const std::complex<float> *bv,
|
||||
std::complex<float> 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<float> *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<float> 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<float> *) in;
|
||||
outv = (float *)out;
|
||||
while (n) {
|
||||
std::complex<float> 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<double> *c,
|
||||
const std::complex<double> *a,
|
||||
const std::complex<double> *b,
|
||||
size_t n)
|
||||
{
|
||||
elementwise_product_naive(c, a, b, n);
|
||||
}
|
||||
|
||||
static inline void
|
||||
elementwise_product(std::complex<double> *c,
|
||||
const std::complex<double> *a,
|
||||
const double *b,
|
||||
size_t n)
|
||||
{
|
||||
elementwise_product_naive(c, a, b, n);
|
||||
}
|
||||
|
||||
template <class T, class U, class V, class S>
|
||||
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 <class O>
|
||||
static inline void
|
||||
complex_magnitude(std::complex<double> *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 <class T, class U, class V>
|
||||
static inline void
|
||||
elementwise_product(T *r,
|
||||
U *a,
|
||||
V *b,
|
||||
size_t n)
|
||||
{
|
||||
elementwise_product_naive(r, a, b, n);
|
||||
}
|
||||
|
||||
template <class T, class U, class V, class S>
|
||||
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 <class I, class O>
|
||||
static inline void
|
||||
complex_magnitude(I *inv,
|
||||
O *outv,
|
||||
size_t n)
|
||||
{
|
||||
complex_magnitude_naive(inv, outv, n);
|
||||
}
|
||||
|
||||
#endif // ! USE_SSE3_INTINSICS
|
||||
|
||||
} // namespace
|
||||
|
||||
#endif
|
|
@ -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
|
Loading…
Reference in a new issue