diff options
| -rw-r--r-- | CMakeLists.txt | 139 | ||||
| -rw-r--r-- | LICENCE | 502 | ||||
| -rw-r--r-- | README.md | 73 | ||||
| -rw-r--r-- | cmake/Modules/FindCsound.cmake | 29 | ||||
| -rw-r--r-- | examples/transformhelpers.udo | 99 | ||||
| -rw-r--r-- | include/cosine.hpp | 29 | ||||
| -rw-r--r-- | include/haar.hpp | 29 | ||||
| -rw-r--r-- | include/walsh.hpp | 43 | ||||
| -rw-r--r-- | include/walsh2.hpp | 30 | ||||
| -rw-r--r-- | src/cosine.cpp | 131 | ||||
| -rw-r--r-- | src/haar.cpp | 418 | ||||
| -rw-r--r-- | src/opcodes.cpp | 220 | ||||
| -rw-r--r-- | src/walsh.cpp | 1081 | ||||
| -rw-r--r-- | src/walsh2.cpp | 58 | ||||
| -rw-r--r-- | src/windowing.cpp | 208 |
15 files changed, 3089 insertions, 0 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..42753d0 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,139 @@ +project("csound_transforms") + +cmake_minimum_required(VERSION 3.8) + +set(APIVERSION "6.0") + +# Release or Debug +set(CMAKE_BUILD_TYPE "Release") + +# force make to print the command lines +set(CMAKE_VERBOSE_MAKEFILE on) + +# path to Csound cmake module +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules/") + +# set compilation flags +set(CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS} "-std=c++11 -fpermissive -fPIC -w -DUSE_DOUBLE -DB64BIT") + +# options +option(USE_LIB64 "Set to on to set installation dir for libs to lib64" OFF) +option(USE_DOUBLE "Use doubles for audio calculations" ON) +option(CPP11 "c++11" ON) + +set(BUILDING_CSOUND_PLUGINS ON) + +# ---------------------------------------------- + +include(FindCsound) + +include(CheckCCompilerFlag) +include(CheckCXXCompilerFlag) + +# ----------------------------------------------- + +function(addflag flag flagname) + check_c_compiler_flag(${flag} ${flagname}) + if (${flagname}) + # message(STATUS "Setting C flag ${flag}") + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${flag}" PARENT_SCOPE) + endif() + check_cxx_compiler_flag(${flag} CXX_${flagname}) + if (CXX_${flagname}) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${flag}" PARENT_SCOPE) + endif() +endfunction(addflag) + + +# set optimization flags +if(CMAKE_COMPILER_IS_GNUCXX OR CMAKE_COMPILER_IS_CLANG) + add_definitions(-fvisibility=hidden) + if(NATIVE) + add_definitions(-march=native) + endif() + + include(CheckCCompilerFlag) + include(CheckCXXCompilerFlag) + + addflag(-msse HAS_SSE) + addflag(-msse2 HAS_SSE2) + addflag(-mfgpath=sse HAS_FPMATH_SSE) + +endif() + +if(MINGW) + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -mstackrealign") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -mstackrealign") +endif() + +addflag(-ftree-vectorize HAS_TREE_VECTORIZE) +addflag(-ffast-math HAS_FAST_MATH) +addflag(-fomit-frame-pointer HAS_OMIT_FRAME_POINTER) + + +# ------------------------------------------------------------------- + +set(CS_FRAMEWORK_DEST "~/Library/Frameworks") + + +if(USE_LIB64) + set(LIBRARY_INSTALL_DIR "lib64") + add_definitions("-DLIB64") +else() + set(LIBRARY_INSTALL_DIR "lib") +endif() + +message(STATUS "LIBRARY INSTALL DIR: ${LIBRARY_INSTALL_DIR}") + +# ------------------------------------------------------------------- + + +if(USE_DOUBLE) + message(STATUS ">>> using doubles") + + if(APPLE) + set(CSOUNDLIB "CsoundLib64") + set(PLUGIN_INSTALL_DIR "${CS_FRAMEWORK_DEST}/${CSOUNDLIB}.framework/Versions/${APIVERSION}/Resources/Opcodes64") + else() + set(CSOUNDLIB "csound64") + set(PLUGIN_INSTALL_DIR "${LIBRARY_INSTALL_DIR}/csound/plugins64-${APIVERSION}") + endif() +else() + message(STATUS ">>> not using doubles") + if(APPLE) + set(CSOUNDLIB "CsoundLib") + set(PLUGIN_INSTALL_DIR "${CS_FRAMEWORK_DEST}/${CSOUNDLIB}.framework/Versions/${APIVERSION}/Resources/Opcodes") + else() + set(CSOUNDLIB "csound") + set(PLUGIN_INSTALL_DIR "${LIBRARY_INSTALL_DIR}/csound/plugins-${APIVERSION}") + endif() +endif() + + +# ------------------------------------------------------------------- + +# Csound opcode build +find_package(Csound) + + +set(BUILD_PLUGINS_DIR ${CMAKE_CURRENT_BINARY_DIR}) + +if(NOT CSOUND_FOUND) + message(FATAL_ERROR "Csound installation not found") +endif() + + +set(CPPFILES src/opcodes.cpp src/cosine.cpp src/haar.cpp src/walsh.cpp src/walsh2.cpp) + + +add_library(transforms SHARED ${CPPFILES}) +include_directories(${CSOUND_INCLUDE_DIRS}) +include_directories(include) + + +set_target_properties(transforms PROPERTIES + RUNTIME_OUTPUT_DIRECTORY ${BUILD_PLUGINS_DIR} + LIBRARY_OUTPUT_DIRECTORY ${BUILD_PLUGINS_DIR}) + +install(TARGETS transforms LIBRARY DESTINATION "${PLUGIN_INSTALL_DIR}" ) + @@ -0,0 +1,502 @@ + GNU LESSER GENERAL PUBLIC LICENSE + Version 2.1, February 1999 + + Copyright (C) 1991, 1999 Free Software Foundation, Inc. + 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + +[This is the first released version of the Lesser GPL. It also counts + as the successor of the GNU Library Public License, version 2, hence + the version number 2.1.] + + Preamble + + The licenses for most software are designed to take away your +freedom to share and change it. By contrast, the GNU General Public +Licenses are intended to guarantee your freedom to share and change +free software--to make sure the software is free for all its users. + + This license, the Lesser General Public License, applies to some +specially designated software packages--typically libraries--of the +Free Software Foundation and other authors who decide to use it. You +can use it too, but we suggest you first think carefully about whether +this license or the ordinary General Public License is the better +strategy to use in any particular case, based on the explanations below. + + When we speak of free software, we are referring to freedom of use, +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 this service if you wish); that you receive source code or can get +it if you want it; that you can change the software and use pieces of +it in new free programs; and that you are informed that you can do +these things. + + To protect your rights, we need to make restrictions that forbid +distributors to deny you these rights or to ask you to surrender these +rights. These restrictions translate to certain responsibilities for +you if you distribute copies of the library or if you modify it. + + For example, if you distribute copies of the library, whether gratis +or for a fee, you must give the recipients all the rights that we gave +you. You must make sure that they, too, receive or can get the source +code. If you link other code with the library, you must provide +complete object files to the recipients, so that they can relink them +with the library after making changes to the library and recompiling +it. And you must show them these terms so they know their rights. + + We protect your rights with a two-step method: (1) we copyright the +library, and (2) we offer you this license, which gives you legal +permission to copy, distribute and/or modify the library. + + To protect each distributor, we want to make it very clear that +there is no warranty for the free library. Also, if the library is +modified by someone else and passed on, the recipients should know +that what they have is not the original version, so that the original +author's reputation will not be affected by problems that might be +introduced by others. +^L + Finally, software patents pose a constant threat to the existence of +any free program. We wish to make sure that a company cannot +effectively restrict the users of a free program by obtaining a +restrictive license from a patent holder. Therefore, we insist that +any patent license obtained for a version of the library must be +consistent with the full freedom of use specified in this license. + + Most GNU software, including some libraries, is covered by the +ordinary GNU General Public License. This license, the GNU Lesser +General Public License, applies to certain designated libraries, and +is quite different from the ordinary General Public License. We use +this license for certain libraries in order to permit linking those +libraries into non-free programs. + + When a program is linked with a library, whether statically or using +a shared library, the combination of the two is legally speaking a +combined work, a derivative of the original library. The ordinary +General Public License therefore permits such linking only if the +entire combination fits its criteria of freedom. The Lesser General +Public License permits more lax criteria for linking other code with +the library. + + We call this license the "Lesser" General Public License because it +does Less to protect the user's freedom than the ordinary General +Public License. It also provides other free software developers Less +of an advantage over competing non-free programs. These disadvantages +are the reason we use the ordinary General Public License for many +libraries. However, the Lesser license provides advantages in certain +special circumstances. + + For example, on rare occasions, there may be a special need to +encourage the widest possible use of a certain library, so that it becomes +a de-facto standard. To achieve this, non-free programs must be +allowed to use the library. A more frequent case is that a free +library does the same job as widely used non-free libraries. In this +case, there is little to gain by limiting the free library to free +software only, so we use the Lesser General Public License. + + In other cases, permission to use a particular library in non-free +programs enables a greater number of people to use a large body of +free software. For example, permission to use the GNU C Library in +non-free programs enables many more people to use the whole GNU +operating system, as well as its variant, the GNU/Linux operating +system. + + Although the Lesser General Public License is Less protective of the +users' freedom, it does ensure that the user of a program that is +linked with the Library has the freedom and the wherewithal to run +that program using a modified version of the Library. + + The precise terms and conditions for copying, distribution and +modification follow. Pay close attention to the difference between a +"work based on the library" and a "work that uses the library". The +former contains code derived from the library, whereas the latter must +be combined with the library in order to run. +^L + GNU LESSER GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License Agreement applies to any software library or other +program which contains a notice placed by the copyright holder or +other authorized party saying it may be distributed under the terms of +this Lesser General Public License (also called "this License"). +Each licensee is addressed as "you". + + A "library" means a collection of software functions and/or data +prepared so as to be conveniently linked with application programs +(which use some of those functions and data) to form executables. + + The "Library", below, refers to any such software library or work +which has been distributed under these terms. A "work based on the +Library" means either the Library or any derivative work under +copyright law: that is to say, a work containing the Library or a +portion of it, either verbatim or with modifications and/or translated +straightforwardly into another language. (Hereinafter, translation is +included without limitation in the term "modification".) + + "Source code" for a work means the preferred form of the work for +making modifications to it. For a library, complete source code means +all the source code for all modules it contains, plus any associated +interface definition files, plus the scripts used to control compilation +and installation of the library. + + Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running a program using the Library is not restricted, and output from +such a program is covered only if its contents constitute a work based +on the Library (independent of the use of the Library in a tool for +writing it). Whether that is true depends on what the Library does +and what the program that uses the Library does. + + 1. You may copy and distribute verbatim copies of the Library's +complete source code as you receive it, in any medium, provided that +you conspicuously and appropriately publish on each copy an +appropriate copyright notice and disclaimer of warranty; keep intact +all the notices that refer to this License and to the absence of any +warranty; and distribute a copy of this License along with the +Library. + + You may charge a fee for the physical act of transferring a copy, +and you may at your option offer warranty protection in exchange for a +fee. +^L + 2. You may modify your copy or copies of the Library or any portion +of it, thus forming a work based on the Library, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) The modified work must itself be a software library. + + b) You must cause the files modified to carry prominent notices + stating that you changed the files and the date of any change. + + c) You must cause the whole of the work to be licensed at no + charge to all third parties under the terms of this License. + + d) If a facility in the modified Library refers to a function or a + table of data to be supplied by an application program that uses + the facility, other than as an argument passed when the facility + is invoked, then you must make a good faith effort to ensure that, + in the event an application does not supply such function or + table, the facility still operates, and performs whatever part of + its purpose remains meaningful. + + (For example, a function in a library to compute square roots has + a purpose that is entirely well-defined independent of the + application. Therefore, Subsection 2d requires that any + application-supplied function or table used by this function must + be optional: if the application does not supply it, the square + root function must still compute square roots.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Library, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Library, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote +it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Library. + +In addition, mere aggregation of another work not based on the Library +with the Library (or with a work based on the Library) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may opt to apply the terms of the ordinary GNU General Public +License instead of this License to a given copy of the Library. To do +this, you must alter all the notices that refer to this License, so +that they refer to the ordinary GNU General Public License, version 2, +instead of to this License. (If a newer version than version 2 of the +ordinary GNU General Public License has appeared, then you can specify +that version instead if you wish.) Do not make any other change in +these notices. +^L + Once this change is made in a given copy, it is irreversible for +that copy, so the ordinary GNU General Public License applies to all +subsequent copies and derivative works made from that copy. + + This option is useful when you wish to copy part of the code of +the Library into a program that is not a library. + + 4. You may copy and distribute the Library (or a portion or +derivative of it, under Section 2) in object code or executable form +under the terms of Sections 1 and 2 above provided that you accompany +it with the complete corresponding machine-readable source code, which +must be distributed under the terms of Sections 1 and 2 above on a +medium customarily used for software interchange. + + If distribution of object code is made by offering access to copy +from a designated place, then offering equivalent access to copy the +source code from the same place satisfies the requirement to +distribute the source code, even though third parties are not +compelled to copy the source along with the object code. + + 5. A program that contains no derivative of any portion of the +Library, but is designed to work with the Library by being compiled or +linked with it, is called a "work that uses the Library". Such a +work, in isolation, is not a derivative work of the Library, and +therefore falls outside the scope of this License. + + However, linking a "work that uses the Library" with the Library +creates an executable that is a derivative of the Library (because it +contains portions of the Library), rather than a "work that uses the +library". The executable is therefore covered by this License. +Section 6 states terms for distribution of such executables. + + When a "work that uses the Library" uses material from a header file +that is part of the Library, the object code for the work may be a +derivative work of the Library even though the source code is not. +Whether this is true is especially significant if the work can be +linked without the Library, or if the work is itself a library. The +threshold for this to be true is not precisely defined by law. + + If such an object file uses only numerical parameters, data +structure layouts and accessors, and small macros and small inline +functions (ten lines or less in length), then the use of the object +file is unrestricted, regardless of whether it is legally a derivative +work. (Executables containing this object code plus portions of the +Library will still fall under Section 6.) + + Otherwise, if the work is a derivative of the Library, you may +distribute the object code for the work under the terms of Section 6. +Any executables containing that work also fall under Section 6, +whether or not they are linked directly with the Library itself. +^L + 6. As an exception to the Sections above, you may also combine or +link a "work that uses the Library" with the Library to produce a +work containing portions of the Library, and distribute that work +under terms of your choice, provided that the terms permit +modification of the work for the customer's own use and reverse +engineering for debugging such modifications. + + You must give prominent notice with each copy of the work that the +Library is used in it and that the Library and its use are covered by +this License. You must supply a copy of this License. If the work +during execution displays copyright notices, you must include the +copyright notice for the Library among them, as well as a reference +directing the user to the copy of this License. Also, you must do one +of these things: + + a) Accompany the work with the complete corresponding + machine-readable source code for the Library including whatever + changes were used in the work (which must be distributed under + Sections 1 and 2 above); and, if the work is an executable linked + with the Library, with the complete machine-readable "work that + uses the Library", as object code and/or source code, so that the + user can modify the Library and then relink to produce a modified + executable containing the modified Library. (It is understood + that the user who changes the contents of definitions files in the + Library will not necessarily be able to recompile the application + to use the modified definitions.) + + b) Use a suitable shared library mechanism for linking with the + Library. A suitable mechanism is one that (1) uses at run time a + copy of the library already present on the user's computer system, + rather than copying library functions into the executable, and (2) + will operate properly with a modified version of the library, if + the user installs one, as long as the modified version is + interface-compatible with the version that the work was made with. + + c) Accompany the work with a written offer, valid for at + least three years, to give the same user the materials + specified in Subsection 6a, above, for a charge no more + than the cost of performing this distribution. + + d) If distribution of the work is made by offering access to copy + from a designated place, offer equivalent access to copy the above + specified materials from the same place. + + e) Verify that the user has already received a copy of these + materials or that you have already sent this user a copy. + + For an executable, the required form of the "work that uses the +Library" must include any data and utility programs needed for +reproducing the executable from it. However, as a special exception, +the materials to be distributed need not include anything that is +normally distributed (in either source or binary form) with the major +components (compiler, kernel, and so on) of the operating system on +which the executable runs, unless that component itself accompanies +the executable. + + It may happen that this requirement contradicts the license +restrictions of other proprietary libraries that do not normally +accompany the operating system. Such a contradiction means you cannot +use both them and the Library together in an executable that you +distribute. +^L + 7. You may place library facilities that are a work based on the +Library side-by-side in a single library together with other library +facilities not covered by this License, and distribute such a combined +library, provided that the separate distribution of the work based on +the Library and of the other library facilities is otherwise +permitted, and provided that you do these two things: + + a) Accompany the combined library with a copy of the same work + based on the Library, uncombined with any other library + facilities. This must be distributed under the terms of the + Sections above. + + b) Give prominent notice with the combined library of the fact + that part of it is a work based on the Library, and explaining + where to find the accompanying uncombined form of the same work. + + 8. You may not copy, modify, sublicense, link with, or distribute +the Library except as expressly provided under this License. Any +attempt otherwise to copy, modify, sublicense, link with, or +distribute the Library is void, and will automatically terminate your +rights under this License. However, parties who have received copies, +or rights, from you under this License will not have their licenses +terminated so long as such parties remain in full compliance. + + 9. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Library or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Library (or any work based on the +Library), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Library or works based on it. + + 10. Each time you redistribute the Library (or any work based on the +Library), the recipient automatically receives a license from the +original licensor to copy, distribute, link with or modify the Library +subject to these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties with +this License. +^L + 11. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +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 +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Library at all. For example, if a patent +license would not permit royalty-free redistribution of the Library by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Library. + +If any portion of this section is held invalid or unenforceable under any +particular circumstance, the balance of the section is intended to apply, +and the section as a whole is intended to apply in other circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 12. If the distribution and/or use of the Library is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Library under this License may add +an explicit geographical distribution limitation excluding those countries, +so that distribution is permitted only in or among countries not thus +excluded. In such case, this License incorporates the limitation as if +written in the body of this License. + + 13. The Free Software Foundation may publish revised and/or new +versions of the Lesser 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 Library +specifies a version number of this License which applies to it and +"any later version", you have the option of following the terms and +conditions either of that version or of any later version published by +the Free Software Foundation. If the Library does not specify a +license version number, you may choose any version ever published by +the Free Software Foundation. +^L + 14. If you wish to incorporate parts of the Library into other free +programs whose distribution conditions are incompatible with these, +write to the author to ask for permission. For software which is +copyrighted by the Free Software Foundation, write to the Free +Software Foundation; we sometimes make exceptions for this. Our +decision will be guided by the two goals of preserving the free status +of all derivatives of our free software and of promoting the sharing +and reuse of software generally. + + NO WARRANTY + + 15. BECAUSE THE LIBRARY IS LICENSED FREE OF CHARGE, THERE IS NO +WARRANTY FOR THE LIBRARY, TO THE EXTENT PERMITTED BY APPLICABLE LAW. +EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR +OTHER PARTIES PROVIDE THE LIBRARY "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 +LIBRARY IS WITH YOU. SHOULD THE LIBRARY PROVE DEFECTIVE, YOU ASSUME +THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN +WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY +AND/OR REDISTRIBUTE THE LIBRARY 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 +LIBRARY (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 LIBRARY TO OPERATE WITH ANY OTHER SOFTWARE), EVEN IF +SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH +DAMAGES. + + END OF TERMS AND CONDITIONS +^L + How to Apply These Terms to Your New Libraries + + If you develop a new library, and you want it to be of the greatest +possible use to the public, we recommend making it free software that +everyone can redistribute and change. You can do so by permitting +redistribution under these terms (or, alternatively, under the terms of the +ordinary General Public License). + + To apply these terms, attach the following notices to the library. It is +safest to attach them to the start of each source file to most effectively +convey 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 library's name and a brief idea of what it does.> + Copyright (C) <year> <name of author> + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library 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 + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Also add information on how to contact you by electronic and paper mail. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the library, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the + library `Frob' (a library for tweaking knobs) written by James Random Hacker. + + <signature of Ty Coon>, 1 April 1990 + Ty Coon, President of Vice + +That's all there is to it! diff --git a/README.md b/README.md new file mode 100644 index 0000000..b91227f --- /dev/null +++ b/README.md @@ -0,0 +1,73 @@ +# Introduction + +These opcodes present variations on the following transforms: + Haar, + Walsh/Hadamard + Cosine + +At current the opcodes will function simply and similar to the fft and fftinv +opcodes in Csound - ie, the user must implement windowing and ksmps in respect +to the desired analysis. + +These are experimental and exact usage/manipulation of the functions is subject +to some testing so there may be further development to make the opcodes more +usable. + +Tested on Linux with Csound 6.14 + + +# Installation + + mkdir build && cd build + cmake .. + make + sudo make install + + +# Opcodes provided +Each of the following opcodes receives and emits a k-rate array. The size of the +input array should be a power of two. The output size will match the input. + + +### tfcosine +Cosine transform. Notably slow (probably not for realtime usage) + +### tfcosineinv +Inverse cosine transform. As Above + +### tfhaar1 +Haar transform, algorithm 1 + +### tfhaar1inv +Inverse Haar transform, algorithm 1 + +### tfhaar2 +Haar transform, algorithm 2 + +### tfhaar2inv +Inverse Haar transform, algorithm 1 + +### tfwalsh1 +Walsh/Hadamard transform, algorithm 1 + +### tfwalsh2 +Walsh/Hadamard transform, algorithm 2 + +### tfwalsh2inv +Inverse Walsh/Hadamard transform, algorithm 2 + + + + +# Credits + +Code, inspiration and relevant licensing has been derived from the following: + https://people.sc.fsu.edu/~jburkardt/f_src/walsh/walsh.html + https://people.sc.fsu.edu/~jburkardt/f_src/haar/haar.html + https://people.sc.fsu.edu/~jburkardt/f_src/cosine_transform/cosine_transform.html + https://github.com/mochow13/competitive-programming-library/blob/master/Math/Fast%20Walsh-Hadamard%20Transform.cpp + + +Todo / possibly forthcoming + Implement Daubechies wavelets, sine transform. + diff --git a/cmake/Modules/FindCsound.cmake b/cmake/Modules/FindCsound.cmake new file mode 100644 index 0000000..e55b269 --- /dev/null +++ b/cmake/Modules/FindCsound.cmake @@ -0,0 +1,29 @@ +# Try to find the Csound library. +# Once done this will define: +# CSOUND_FOUND - System has the Csound library +# CSOUND_INCLUDE_DIRS - The Csound include directories. +# CSOUND_LIBRARIES - The libraries needed to use the Csound library. + +if(APPLE) +find_path(CSOUND_INCLUDE_DIR csound.h HINTS /Library/Frameworks/CsoundLib64.framework/Headers +"$ENV{HOME}/Library/Frameworks/CsoundLib64.framework/Headers") +else() +find_path(CSOUND_INCLUDE_DIR csound.h PATH_SUFFIXES csound) +endif() + +if(APPLE) +find_library(CSOUND_LIBRARY NAMES CsoundLib64 HINTS /Library/Frameworks/CsoundLib64.framework/ +"$ENV{HOME}/Library/Frameworks/CsoundLib64.framework") +else() +find_library(CSOUND_LIBRARY NAMES csound64 csound) +endif() + +include(FindPackageHandleStandardArgs) +# handle the QUIETLY and REQUIRED arguments and set CSOUND_FOUND to TRUE +# if all listed variables are TRUE +find_package_handle_standard_args(CSOUND + CSOUND_LIBRARY CSOUND_INCLUDE_DIR) +mark_as_advanced(CSOUND_INCLUDE_DIR CSOUND_LIBRARY) + +set(CSOUND_INCLUDE_DIRS ${CSOUND_INCLUDE_DIR}) +set(CSOUND_LIBRARIES ${CSOUND_LIBRARY} ) diff --git a/examples/transformhelpers.udo b/examples/transformhelpers.udo new file mode 100644 index 0000000..d3168a1 --- /dev/null +++ b/examples/transformhelpers.udo @@ -0,0 +1,99 @@ + +#define DEFAULT_SIZE #1024# +#define DEFAULT_HOPSIZE #128# + +opcode tfin, k[], ajj + a1, isize, ihopsize xin + if (isize == -1) then + isize = $DEFAULT_SIZE + endif + if (ihopsize == -1) then + ihopsize = $DEFAULT_HOPSIZE + endif + + iolaps = isize/ihopsize + kcnt init 0 + krow init 0 + kTemp[] init isize + kOla[] init isize + kIn[] init isize + kOut[][] init iolaps, isize + + if (kcnt == ihopsize) then + kWin[] window kIn, krow*ihopsize + kSpec[] tfhaar1 kWin + krow = (krow+1) % iolaps + kcnt = 0 + endif + + kIn shiftin a1 + + xout kSpec + kcnt += ksmps +endop + + +opcode tfout, a, k[]jj + kSpec[], isize, ihopsize xin + if (isize == -1) then + isize = $DEFAULT_SIZE + endif + if (ihopsize == -1) then + ihopsize = $DEFAULT_HOPSIZE + endif + + iolaps = isize/ihopsize ; overlaps + kcnt init 0 ; counting vars + krow init 0 + kTemp[] init isize + kOla[] init isize ; overlap-add buffer + kIn[] init isize ; input buffer + kOut[][] init iolaps, isize ; output buffers + + if (kcnt == ihopsize) then + kRow[] tfhaar1inv kSpec + kWin[] window kRow, krow*ihopsize + kOut setrow kWin, krow + kOla = 0 + ki = 0 + until (ki == iolaps) do + kRow getrow kOut, ki + kOla = kOla + kRow + ki += 1 + od + + krow = (krow+1)%iolaps + kcnt = 0 + endif + a2 shiftout kOla + aout = (a2/iolaps) + xout aout + kcnt += ksmps +endop + + +opcode scramble, k[], k[] + kin[] xin + inum random lenarray(kin)*0.2, lenarray(kin) + kndx = 0 + while (kndx < inum) do + ksrc random 0, lenarray(kin) + kdest random 0, lenarray(kin) + kin[int(kdest)] = kin[int(ksrc)] + kndx += 1 + od + xout kin +endop + + +opcode invert, k[], k[] + kin[] xin + kout[] init lenarray(kin) + kwritex = 0 + imax = lenarray(kin) - 1 + while (kwritex < lenarray(kin)) do + kout[kwritex] = kin[imax - kwritex] + kwritex += 1 + od + xout kout +endop diff --git a/include/cosine.hpp b/include/cosine.hpp new file mode 100644 index 0000000..b9a4804 --- /dev/null +++ b/include/cosine.hpp @@ -0,0 +1,29 @@ +/* + cosine.hpp + Copyright (C) 2019 Richard Knight + + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser 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 + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with this program; if not, write to the Free Software Foundation, + Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + */ + +#ifndef COSINE_HPP +#define COSINE_HPP +#include <plugin.h> +void cosine_transform_data (csnd::Vector<MYFLT> d, csnd::Vector<MYFLT> c); +void cosine_transform_inverse (csnd::Vector<MYFLT> c, csnd::Vector<MYFLT> d); + +#endif /* COSINE_HPP */ + diff --git a/include/haar.hpp b/include/haar.hpp new file mode 100644 index 0000000..c18477d --- /dev/null +++ b/include/haar.hpp @@ -0,0 +1,29 @@ +/* + haar.hpp + Copyright (C) 2019 Richard Knight + + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser 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 + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with this program; if not, write to the Free Software Foundation, + Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + */ + +#ifndef HAAR_HPP +#define HAAR_HPP +#include <plugin.h> + +void haar_1d (csnd::Csound* csound, csnd::Vector<MYFLT> x); +void haar_1d_inverse (csnd::Csound* csound, csnd::Vector<MYFLT> x); +#endif /* HAAR_HPP */ + diff --git a/include/walsh.hpp b/include/walsh.hpp new file mode 100644 index 0000000..b3e3d4b --- /dev/null +++ b/include/walsh.hpp @@ -0,0 +1,43 @@ +/* + walsh.hpp + Copyright (C) 2019 Richard Knight + + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser 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 + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with this program; if not, write to the Free Software Foundation, + Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + */ + +#ifndef WALSH_HPP +#define WALSH_HPP +#include <plugin.h> +void ffwt (csnd::Csound* csound, csnd::Vector<MYFLT> x); +void fwt (csnd::Csound* csound, csnd::Vector<MYFLT> x); +void haar (csnd::Csound* csound, csnd::Vector<MYFLT> x); +void haarin (csnd::Csound* csound, csnd::Vector<MYFLT> x); +void hnorm (csnd::Vector<MYFLT> x); +int i4_log_2 ( int i ); +int i4_max ( int i1, int i2 ); +int i4_min ( int i1, int i2 ); +int i4_modp ( int i, int j ); +int i4_power ( int i, int j ); +int i4_wrap ( int ival, int ilo, int ihi ); +void r8vec_copy ( int n, double a1[], csnd::Vector<MYFLT> a2); +void r8vec_copy ( int n, csnd::Vector<MYFLT> a1, double a2[]); +double *r8vec_copy_new (csnd::Csound* csound, int n, double a1[] ); +void walsh (csnd::Csound* csound, csnd::Vector<MYFLT> x); + + +#endif /* WALSH_HPP */ + diff --git a/include/walsh2.hpp b/include/walsh2.hpp new file mode 100644 index 0000000..e7066ed --- /dev/null +++ b/include/walsh2.hpp @@ -0,0 +1,30 @@ +/* + walsh2.cpp + Copyright (C) 2019 Richard Knight + + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser 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 + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with this program; if not, write to the Free Software Foundation, + Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + */ +#ifndef WALSH2_HPP +#define WALSH2_HPP + +#include <plugin.h> + +void fwtanal(csnd::Vector<MYFLT> io); +void fwtsynth(csnd::Vector<MYFLT> io); + +#endif /* WALSH2_HPP */ + diff --git a/src/cosine.cpp b/src/cosine.cpp new file mode 100644 index 0000000..78cafd8 --- /dev/null +++ b/src/cosine.cpp @@ -0,0 +1,131 @@ +/* + cosine.cpp + Copyright (C) 2015 John Burkardt + Copyright (C) 2019 Richard Knight + + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser 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 + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with this program; if not, write to the Free Software Foundation, + Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + */ + +#include <cmath> +#include <cstdlib> +#include <iostream> +#include <ctime> +#include <plugin.h> + +using namespace std; + +#include "cosine.hpp" + +//****************************************************************************80 + +void cosine_transform_data (csnd::Vector<MYFLT> d, csnd::Vector<MYFLT> c) + +//****************************************************************************80 +// +// Purpose: +// +// COSINE_TRANSFORM_DATA does a cosine transform on a vector of data. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 27 August 2015 +// +// Author: +// +// John Burkardt +// +// Parameters: +// +// Input, integer N, the number of data points. +// +// Input, double D[N], the vector of data. +// +// Output, double COSINE_TRANSFORM_DATA[N], the transform coefficients. +// +{ + int n = d.len(); + double angle; + int i; + int j; + const double r8_pi = 3.141592653589793; + + + for ( i = 0; i < n; i++ ) + { + c[i] = 0.0; + for ( j = 0; j < n; j++ ) + { + angle = r8_pi * ( double ) ( i * ( 2 * j + 1 ) ) / ( double ) ( 2 * n ); + c[i] = c[i] + cos ( angle ) * d[j]; + } + c[i] = c[i] * sqrt ( 2.0 / ( double ) ( n ) ); + } +} +//****************************************************************************80 + +void cosine_transform_inverse (csnd::Vector<MYFLT> c, csnd::Vector<MYFLT> d) + +//****************************************************************************80 +// +// Purpose: +// +// COSINE_TRANSFORM_INVERSE does an inverse cosine transform. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 27 August 2015 +// +// Author: +// +// John Burkardt +// +// Parameters: +// +// Input, integer N, the number of data points. +// +// Input, double C[N], the vector of transform coefficients +// +// Output, double COSINE_TRANSFORM_INVERSE[N], the original data. +// +{ + int n = c.len(); + double angle; + int i; + int j; + double r8_pi = 3.141592653589793; + + + for ( i = 0; i < n; i++ ) + { + d[i] = c[0] / 2.0; + for ( j = 1; j < n; j++ ) + { + angle = r8_pi * ( double ) ( ( 2 * i + 1 ) * j ) / ( double ) ( 2 * n ); + d[i] = d[i] + cos ( angle ) * c[j]; + } + d[i] = d[i] * sqrt ( 2.0 / ( double ) ( n ) ); + } +} +//****************************************************************************80 diff --git a/src/haar.cpp b/src/haar.cpp new file mode 100644 index 0000000..36d39da --- /dev/null +++ b/src/haar.cpp @@ -0,0 +1,418 @@ +/* + haar.cpp + Copyright (C) 2014 John Burkardt + Copyright (C) 2019 Richard Knight + + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser 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 + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with this program; if not, write to the Free Software Foundation, + Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + */ + +#include <cstdlib> +#include <iostream> +#include <iomanip> +#include <cmath> +#include <ctime> +#include <plugin.h> + +using namespace std; + +#include "haar.hpp" + +//****************************************************************************80 + +void haar_1d (csnd::Csound* csound, csnd::Vector<MYFLT> x) + +//****************************************************************************80 +// +// Purpose: +// +// HAAR_1D computes the Haar transform of a vector. +// +// Discussion: +// +// For the classical Haar transform, N should be a power of 2. +// However, this is not required here. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 06 March 2014 +// +// Author: +// +// John Burkardt +// +// Parameters: +// +// Input, int N, the dimension of the vector. +// +// Input/output, double X[N], on input, the vector to be transformed. +// On output, the transformed vector. +// +{ + int n = x.len(); + int i; + int k; + double s; + double *y; + + s = sqrt ( 2.0 ); + + //y = new double[n]; + y = (double*) csound->malloc(sizeof(double) * n); +// +// Initialize. +// + for ( i = 0; i < n; i++ ) + { + y[i] = 0.0; + } +// +// Determine K, the largest power of 2 such that K <= N. +// + k = 1; + while ( k * 2 <= n ) + { + k = k * 2; + } + + while ( 1 < k ) + { + k = k / 2; + for ( i = 0; i < k; i++ ) + { + y[i] = ( x[2*i] + x[2*i+1] ) / s; + y[i+k] = ( x[2*i] - x[2*i+1] ) / s; + } + for ( i = 0; i < k * 2; i++ ) + { + x[i] = y[i]; + } + } + + //delete [] y; + csound->free(y); + + return; +} +//****************************************************************************80 + +void haar_1d_inverse (csnd::Csound* csound, csnd::Vector<MYFLT> x) + +//****************************************************************************80 +// +// Purpose: +// +// HAAR_1D_INVERSE computes the inverse Haar transform of a vector. +// +// Discussion: +// +// For the classical Haar transform, N should be a power of 2. +// However, this is not required here. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 06 March 2014 +// +// Author: +// +// John Burkardt +// +// Parameters: +// +// Input, int N, the dimension of the vector. +// +// Input/output, double X[N], on input, the vector to be transformed. +// On output, the transformed vector. +// +{ + + int n = x.len(); + int i; + int k; + double s; + double *y; + + s = sqrt ( 2.0 ); + + //y = new double[n]; + y = (double*) csound->malloc(sizeof(double) * n); +// +// Initialize. +// + for ( i = 0; i < n; i++ ) + { + y[i] = 0.0; + } + + k = 1; + while ( k * 2 <= n ) + { + for ( i = 0; i < k; i++ ) + { + y[2*i] = ( x[i] + x[i+k] ) / s; + y[2*i+1] = ( x[i] - x[i+k] ) / s; + } + for ( i = 0; i < k * 2; i++ ) + { + x[i] = y[i]; + } + k = k * 2; + } + + //delete [] y; + csound->free(y); + return; +} + + + +/* + + +//****************************************************************************80 + +void haar_2d ( int m, int n, double u[] ) + +//****************************************************************************80 +// +// Purpose: +// +// HAAR_2D computes the Haar transform of an array. +// +// Discussion: +// +// For the classical Haar transform, M and N should be a power of 2. +// However, this is not required here. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 06 March 2014 +// +// Author: +// +// John Burkardt +// +// Parameters: +// +// Input, int M, N, the dimensions of the array. +// +// Input/output, double U[M*N], the array to be transformed. +// +{ + int i; + int j; + int k; + double s; + double *v; + + s = sqrt ( 2.0 ); + + v = new double[m*n]; + + for ( j = 0; j < n; j++ ) + { + for ( i = 0; i < m; i++ ) + { + v[i+j*m] = u[i+j*m]; + } + } +// +// Determine K, the largest power of 2 such that K <= M. +// + k = 1; + while ( k * 2 <= m ) + { + k = k * 2; + } +// +// Transform all columns. +// + while ( 1 < k ) + { + k = k / 2; + + for ( j = 0; j < n; j++ ) + { + for ( i = 0; i < k; i++ ) + { + v[i +j*m] = ( u[2*i+j*m] + u[2*i+1+j*m] ) / s; + v[k+i+j*m] = ( u[2*i+j*m] - u[2*i+1+j*m] ) / s; + } + } + for ( j = 0; j < n; j++ ) + { + for ( i = 0; i < 2 * k; i++ ) + { + u[i+j*m] = v[i+j*m]; + } + } + } +// +// Determine K, the largest power of 2 such that K <= N. +// + k = 1; + while ( k * 2 <= n ) + { + k = k * 2; + } +// +// Transform all rows. +// + while ( 1 < k ) + { + k = k / 2; + + for ( j = 0; j < k; j++ ) + { + for ( i = 0; i < m; i++ ) + { + v[i+( j)*m] = ( u[i+2*j*m] + u[i+(2*j+1)*m] ) / s; + v[i+(k+j)*m] = ( u[i+2*j*m] - u[i+(2*j+1)*m] ) / s; + } + } + + for ( j = 0; j < 2 * k; j++ ) + { + for ( i = 0; i < m; i++ ) + { + u[i+j*m] = v[i+j*m]; + } + } + } + delete [] v; + + return; +} +//****************************************************************************80 + +void haar_2d_inverse ( int m, int n, double u[] ) + +//****************************************************************************80 +// +// Purpose: +// +// HAAR_2D_INVERSE inverts the Haar transform of an array. +// +// Discussion: +// +// For the classical Haar transform, M and N should be a power of 2. +// However, this is not required here. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 06 March 2014 +// +// Author: +// +// John Burkardt +// +// Parameters: +// +// Input, int M, N, the dimensions of the array. +// +// Input/output, double U[M*N], the array to be transformed. +// +{ + int i; + int j; + int k; + double s; + double *v; + + s = sqrt ( 2.0 ); + + v = new double[m*n]; + + for ( j = 0; j < n; j++ ) + { + for ( i = 0; i < m; i++ ) + { + v[i+j*m] = u[i+j*m]; + } + } +// +// Inverse transform of all rows. +// + k = 1; + + while ( k * 2 <= n ) + { + for ( j = 0; j < k; j++ ) + { + for ( i = 0; i < m; i++ ) + { + v[i+(2*j )*m] = ( u[i+j*m] + u[i+(k+j)*m] ) / s; + v[i+(2*j+1)*m] = ( u[i+j*m] - u[i+(k+j)*m] ) / s; + } + } + + for ( j = 0; j < 2 * k; j++ ) + { + for ( i = 0; i < m; i++ ) + { + u[i+j*m] = v[i+j*m]; + } + } + k = k * 2; + } +// +// Inverse transform of all columns. +// + k = 1; + + while ( k * 2 <= m ) + { + for ( j = 0; j < n; j++ ) + { + for ( i = 0; i < k; i++ ) + { + v[2*i +j*m] = ( u[i+j*m] + u[k+i+j*m] ) / s; + v[2*i+1+j*m] = ( u[i+j*m] - u[k+i+j*m] ) / s; + } + } + + for ( j = 0; j < n; j++ ) + { + for ( i = 0; i < 2 * k; i++ ) + { + u[i+j*m] = v[i+j*m]; + } + } + k = k * 2; + } + delete [] v; + + return; +} + + * + */
\ No newline at end of file diff --git a/src/opcodes.cpp b/src/opcodes.cpp new file mode 100644 index 0000000..ef4f2b1 --- /dev/null +++ b/src/opcodes.cpp @@ -0,0 +1,220 @@ +/* + opcodes.cpp + Copyright (C) 2019 Richard Knight + + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser 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 + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with this program; if not, write to the Free Software Foundation, + Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + */ + +#include <plugin.h> +#include "cosine.hpp" +#include "haar.hpp" +#include "walsh.hpp" +#include "walsh2.hpp" + + +struct tfcosine : csnd::Plugin<1, 1> { + static constexpr char const *otypes = "k[]"; + static constexpr char const *itypes = "k[]"; + + int init() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + out.init(csound, in.len()); + return OK; + } + + int kperf() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + cosine_transform_data(in, out); + return OK; + } +}; + +struct tfcosineinv : csnd::Plugin<1, 1> { + static constexpr char const *otypes = "k[]"; + static constexpr char const *itypes = "k[]"; + + int init() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + out.init(csound, in.len()); + return OK; + } + + int kperf() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + cosine_transform_inverse(in, out); + return OK; + } +}; + +struct tfhaar1 : csnd::Plugin<1, 1> { + static constexpr char const *otypes = "k[]"; + static constexpr char const *itypes = "k[]"; + + int init() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + out.init(csound, in.len()); + return OK; + } + + int kperf() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + std::copy(in.begin(), in.end(), out.begin()); + haar_1d(csound, out); + return OK; + } +}; + +struct tfhaar1inv : csnd::Plugin<1, 1> { + static constexpr char const *otypes = "k[]"; + static constexpr char const *itypes = "k[]"; + + int init() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + out.init(csound, in.len()); + return OK; + } + + int kperf() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + std::copy(in.begin(), in.end(), out.begin()); + haar_1d_inverse(csound, out); + return OK; + } +}; + +struct tfhaar2 : csnd::Plugin<1, 1> { + static constexpr char const *otypes = "k[]"; + static constexpr char const *itypes = "k[]"; + + int init() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + out.init(csound, in.len()); + return OK; + } + + int kperf() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + std::copy(in.begin(), in.end(), out.begin()); + haar(csound, out); + return OK; + } +}; + +struct tfhaar2inv : csnd::Plugin<1, 1> { + static constexpr char const *otypes = "k[]"; + static constexpr char const *itypes = "k[]"; + + int init() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + out.init(csound, in.len()); + return OK; + } + + int kperf() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + std::copy(in.begin(), in.end(), out.begin()); + haarin(csound, out); + return OK; + } +}; + +struct tfwalsh1 : csnd::Plugin<1, 1> { + static constexpr char const *otypes = "k[]"; + static constexpr char const *itypes = "k[]"; + + int init() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + out.init(csound, in.len()); + return OK; + } + + int kperf() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + std::copy(in.begin(), in.end(), out.begin()); + walsh(csound, out); + return OK; + } +}; + +struct tfwalsh2 : csnd::Plugin<1, 1> { + static constexpr char const *otypes = "k[]"; + static constexpr char const *itypes = "k[]"; + + int init() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + out.init(csound, in.len()); + return OK; + } + + int kperf() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + std::copy(in.begin(), in.end(), out.begin()); + fwtanal(out); + return OK; + } +}; + +struct tfwalsh2inv : csnd::Plugin<1, 1> { + static constexpr char const *otypes = "k[]"; + static constexpr char const *itypes = "k[]"; + + int init() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + out.init(csound, in.len()); + return OK; + } + + int kperf() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + std::copy(in.begin(), in.end(), out.begin()); + fwtsynth(out); + return OK; + } +}; + + +#include <modload.h> + +void csnd::on_load(csnd::Csound *csound) { + csnd::plugin<tfcosine>(csound, "tfcosine", csnd::thread::ik); + csnd::plugin<tfcosineinv>(csound, "tfcosineinv", csnd::thread::ik); + csnd::plugin<tfhaar1>(csound, "tfhaar1", csnd::thread::ik); + csnd::plugin<tfhaar1inv>(csound, "tfhaar1inv", csnd::thread::ik); + csnd::plugin<tfhaar2>(csound, "tfhaar2", csnd::thread::ik); + csnd::plugin<tfhaar2inv>(csound, "tfhaar2inv", csnd::thread::ik); + csnd::plugin<tfwalsh1>(csound, "tfwalsh1", csnd::thread::ik); + csnd::plugin<tfwalsh2>(csound, "tfwalsh2", csnd::thread::ik); + csnd::plugin<tfwalsh2inv>(csound, "tfwalsh2inv", csnd::thread::ik); +} diff --git a/src/walsh.cpp b/src/walsh.cpp new file mode 100644 index 0000000..0ae574a --- /dev/null +++ b/src/walsh.cpp @@ -0,0 +1,1081 @@ +/* + walsh.cpp + Copyright (C) 2004 John Burkardt + Copyright (C) 2011 Ken Beauchamp + Copyright (C) 2019 Richard Knight + + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser 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 + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with this program; if not, write to the Free Software Foundation, + Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + */ + +#include <cstdlib> +#include <iostream> +#include <iomanip> +#include <cmath> +#include <ctime> +#include <cstring> +#include <plugin.h> + +using namespace std; + +#include "walsh.hpp" + +//****************************************************************************80 + +void ffwt (csnd::Csound* csound, csnd::Vector<MYFLT> x) + +//****************************************************************************80 +// +// Purpose: +// +// FFWT performs an in-place fast Walsh transform. +// +// Discussion: +// +// This routine performs a fast Walsh transform on an input series X +// leaving the transformed results in X. +// X is dimensioned N, which must be a power of 2. +// The results of this Walsh transform are in sequency order. +// +// The output sequence could be normalized by dividing by N. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 16 March 2011 +// +// Author: +// +// Ken Beauchamp +// +// Reference: +// +// Ken Beauchamp, +// Walsh functions and their applications, +// Academic Press, 1975, +// ISBN: 0-12-084050-2, +// LC: QA404.5.B33. +// +// Parameters: +// +// Input, int N, the number of items in X. +// N must be a power of 2. +// +// Input/output, double X[N], the data to be transformed. +// +{ + double hold; + int i; + int ii; + int j; + int j2; + int js; + int k; + int l; + int m; + int mw; + int mw1; + int nw; + int nz; + int nz2; + int nzi; + int nzn; + int *two_power; + double z; + int n = x.len(); + + m = i4_log_2 ( n ); + + two_power = (int*) csound->malloc(sizeof(int) * m); // new int[m]; + + for ( i = 0; i < m; i++ ) + { + two_power[i] = i4_power ( 2, m - 1 - i ); + } + + for ( l = 0; l < m; l++ ) + { + nz = i4_power ( 2, l ); + nzi = 2 * nz; + nzn = n / nzi; + nz2 = nz / 2; + if ( nz2 == 0 ) + { + nz2 = 1; + } + + for ( i = 0; i < nzn; i++ ) + { + js = i * nzi; + z = 1.0; + for ( ii = 0; ii < 2; ii++ ) + { + for ( j = 0; j < nz2; j++ ) + { + js = js + 1; + j2 = js + nz; + hold = x[js-1] + z * x[j2-1]; + z = - z; + x[j2-1] = x[js-1] + z * x[j2-1]; + x[js-1] = hold; + z = - z; + } + if ( l == 0 ) + { + break; + } + z = - 1.0; + } + } + } +// +// Bit reversal section. +// + nw = 0; + for ( k = 0; k < n; k++ ) + { +// +// Choose correct index and switch elements if not already switched. +// + if ( k < nw ) + { + hold = x[nw]; + x[nw] = x[k]; + x[k] = hold; + } +// +// Bump up series by 1. +// + for ( i = 0; i < m; i++ ) + { + ii = i; + if ( nw < two_power[i] ) + { + break; + } + mw = nw / two_power[i]; + mw1 = mw / 2; + if ( mw <= 2 * mw1 ) + { + break; + } + nw = nw - two_power[i]; + } + nw = nw + two_power[ii]; + } + + //delete [] two_power; + csound->free(two_power); + + return; +} +//****************************************************************************80 + +void fwt (csnd::Csound* csound, csnd::Vector<MYFLT> x) + +//****************************************************************************80 +// +// Purpose: +// +// FWT performs a fast Walsh transform. +// +// Discussion: +// +// This routine performs a fast Walsh transform on an input series X +// leaving the transformed results in X. +// X is dimensioned N, which must be a power of 2. +// The results of this Walsh transform are in sequency order. +// +// The output sequence could be normalized by dividing by N. +// +// Note that the program text in the reference included the line +// y(jd) = abs ( x(j) - x(j2) ) +// which has been corrected to: +// y(jd) = x(j) - x(j2) +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 16 March 2011 +// +// Author: +// +// Ken Beauchamp +// +// Reference: +// +// Ken Beauchamp, +// Walsh functions and their applications, +// Academic Press, 1975, +// ISBN: 0-12-084050-2, +// LC: QA404.5.B33. +// +// Parameters: +// +// Input, int N, the number of items in X. +// N must be a power of 2. +// +// Input/output, double X[N], the data to be transformed. +// +{ + int i; + int j; + int j2; + int jd; + int js; + int l; + int m; + int n2; + int nx; + int ny; + int nz; + int nzi; + int nzn; + double *y; + int n = x.len(); + y = (double*) csound->malloc(sizeof(double) * n);//new double[n]; + + n2 = n / 2; + m = i4_log_2 ( n ); + + for ( l = 1; l <= m; l++ ) + { + ny = 0; + nz = i4_power ( 2, l - 1 ); + nzi = 2 * nz; + nzn = n / nzi; + for ( i = 1; i <= nzn; i++ ) + { + nx = ny + 1; + ny = ny + nz; + js = ( i - 1 ) * nzi; + jd = js + nzi + 1; + for ( j = nx; j <= ny; j++ ) + { + js = js + 1; + j2 = j + n2; + y[js-1] = x[j-1] + x[j2-1]; + jd = jd - 1; + y[jd-1] = x[j-1] - x[j2-1]; + } + } + r8vec_copy (n, y, x ); + } + csound->free(y); //delete [] y; + + return; +} +//****************************************************************************80 + +void haar (csnd::Csound* csound, csnd::Vector<MYFLT> x) + +//****************************************************************************80 +// +// Purpose: +// +// HAAR performs a Haar transform. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 16 March 2011 +// +// Author: +// +// Ken Beauchamp +// +// Reference: +// +// Ken Beauchamp, +// Walsh functions and their applications, +// Academic Press, 1975, +// ISBN: 0-12-084050-2, +// LC: QA404.5.B33. +// +// Parameters: +// +// Input, int N, the number of items in X. +// N must be a power of 2. +// +// Input/output, double X[N], the data to be transformed. +// +{ + int i; + int i1; + int j; + int jj; + int k; + int l; + int l2; + int l3; + double *y; + int n = x.len(); + + y = (double*) csound->malloc(sizeof(double) * n); //new double[n]; + + k = i4_log_2 ( n ); + + for ( i = 1; i <= k; i++ ) + { + l = k + 1 - i; + l2 = i4_power ( 2, l - 1 ); + + r8vec_copy (2 * l2, x, y ); + + for ( j = 1; j <= l2; j++ ) + { + l3 = l2 + j; + jj = 2 * j - 1; + x[j-1] = y[jj-1] + y[jj]; + x[l3-1] = y[jj-1] - y[jj]; + } + } + csound->free(y); //delete [] y; + + return; +} +//****************************************************************************80 + +void haarin (csnd::Csound* csound, csnd::Vector<MYFLT> x) + +//****************************************************************************80 +// +// Purpose: +// +// HAARIN inverts a Haar transform. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 16 March 2011 +// +// Author: +// +// Ken Beauchamp +// +// Reference: +// +// Ken Beauchamp, +// Walsh functions and their applications, +// Academic Press, 1975, +// ISBN: 0-12-084050-2, +// LC: QA404.5.B33. +// +// Parameters: +// +// Input, int N, the number of items in X. +// N must be a power of 2. +// +// Input/output, double X[N], the data to be transformed. +// +{ + int i; + int i1; + int j; + int jj; + int jj1; + int k; + int l; + int lj; + double *y; + int n = x.len(); + + y = (double*) csound->malloc(sizeof(double) * n); //new double[n]; + + k = i4_log_2 ( n ); + + for ( i = 1; i <= k; i++ ) + { + l = i4_power ( 2, i - 1 ); + r8vec_copy (2 * l, x, y ); + for ( j = 1; j <= l; j++ ) + { + lj = l + j; + jj = 2 * j; + jj1 = jj - 1; + x[jj-1] = y[j-1] - y[lj-1]; + x[jj1-1] = y[j-1] + y[lj-1]; + } + } + + csound->free(y); //delete [] y; + + return; +} +//****************************************************************************80 + +void hnorm (csnd::Vector<MYFLT> x) + +//****************************************************************************80 +// +// Purpose: +// +// HNORM computes normalization factors for a forward or inverse Haar transform. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 16 March 2011 +// +// Author: +// +// Ken Beauchamp +// +// Reference: +// +// Ken Beauchamp, +// Walsh functions and their applications, +// Academic Press, 1975, +// ISBN: 0-12-084050-2, +// LC: QA404.5.B33. +// +// Parameters: +// +// Input, int N, the number of items in X. +// N must be a power of 2. +// +// Input/output, double X[N], the data to be transformed. +// +{ + int i; + int ii; + int j; + int jmax; + int jmin; + int k; + double wlk; + int n = x.len(); + + k = i4_log_2 ( n ); + + x[0] = x[0] / pow ( 2.0, k ); + + if ( 1 <= k ) + { + x[1] = x[1] / pow ( 2.0, k ); + } + + for ( ii = 2; ii <= k; ii++ ) + { + i = ii - 1; + wlk = 1.0 / pow ( 2.0, k - i ); + jmin = i4_power ( 2, i ); + jmax = i4_power ( 2, ii ) - 1; + for ( j = jmin; j <= jmax; j++ ) + { + x[j] = x[j] * wlk; + } + } + return; +} +//****************************************************************************80 + +int i4_log_2 ( int i ) + +//****************************************************************************80 +// +// Purpose: +// +// I4_LOG_2 returns the integer part of the logarithm base 2 of an I4. +// +// Example: +// +// I I4_LOG_10 +// ----- -------- +// 0 0 +// 1 0 +// 2 1 +// 3 1 +// 4 2 +// 5 2 +// 7 2 +// 8 3 +// 9 3 +// 1000 9 +// 1024 10 +// +// Discussion: +// +// I4_LOG_2 ( I ) + 1 is the number of binary digits in I. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 04 January 2004 +// +// Author: +// +// John Burkardt +// +// Parameters: +// +// Input, int I, the number whose logarithm base 2 is desired. +// +// Output, int I4_LOG_2, the integer part of the logarithm base 2 of +// the absolute value of X. +// +{ + int i_abs; + int two_pow; + int value; + + if ( i == 0 ) + { + value = 0; + } + else + { + value = 0; + two_pow = 2; + + i_abs = abs ( i ); + + while ( two_pow <= i_abs ) + { + value = value + 1; + two_pow = two_pow * 2; + } + } + + return value; +} +//****************************************************************************80 + +int i4_max ( int i1, int i2 ) + +//****************************************************************************80 +// +// Purpose: +// +// I4_MAX returns the maximum of two I4's. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 13 October 1998 +// +// Author: +// +// John Burkardt +// +// Parameters: +// +// Input, int I1, I2, are two integers to be compared. +// +// Output, int I4_MAX, the larger of I1 and I2. +// +{ + int value; + + if ( i2 < i1 ) + { + value = i1; + } + else + { + value = i2; + } + return value; +} +//****************************************************************************80 + +int i4_min ( int i1, int i2 ) + +//****************************************************************************80 +// +// Purpose: +// +// I4_MIN returns the minimum of two I4's. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 13 October 1998 +// +// Author: +// +// John Burkardt +// +// Parameters: +// +// Input, int I1, I2, two integers to be compared. +// +// Output, int I4_MIN, the smaller of I1 and I2. +// +{ + int value; + + if ( i1 < i2 ) + { + value = i1; + } + else + { + value = i2; + } + return value; +} +//****************************************************************************80 + +int i4_modp ( int i, int j ) + +//****************************************************************************80 +// +// Purpose: +// +// I4_MODP returns the nonnegative remainder of I4 division. +// +// Discussion: +// +// If +// NREM = I4_MODP ( I, J ) +// NMULT = ( I - NREM ) / J +// then +// I = J * NMULT + NREM +// where NREM is always nonnegative. +// +// The MOD function computes a result with the same sign as the +// quantity being divided. Thus, suppose you had an angle A, +// and you wanted to ensure that it was between 0 and 360. +// Then mod(A,360) would do, if A was positive, but if A +// was negative, your result would be between -360 and 0. +// +// On the other hand, I4_MODP(A,360) is between 0 and 360, always. +// +// I J MOD I4_MODP I4_MODP Factorization +// +// 107 50 7 7 107 = 2 * 50 + 7 +// 107 -50 7 7 107 = -2 * -50 + 7 +// -107 50 -7 43 -107 = -3 * 50 + 43 +// -107 -50 -7 43 -107 = 3 * -50 + 43 +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 26 May 1999 +// +// Author: +// +// John Burkardt +// +// Parameters: +// +// Input, int I, the number to be divided. +// +// Input, int J, the number that divides I. +// +// Output, int I4_MODP, the nonnegative remainder when I is +// divided by J. +// +{ + int value; + + if ( j == 0 ) + { + cerr << "\n"; + cerr << "I4_MODP - Fatal error!\n"; + cerr << " I4_MODP ( I, J ) called with J = " << j << "\n"; + exit ( 1 ); + } + + value = i % j; + + if ( value < 0 ) + { + value = value + abs ( j ); + } + + return value; +} +//****************************************************************************80 + +int i4_power ( int i, int j ) + +//****************************************************************************80 +// +// Purpose: +// +// I4_POWER returns the value of I^J. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 01 April 2004 +// +// Author: +// +// John Burkardt +// +// Parameters: +// +// Input, int I, J, the base and the power. J should be nonnegative. +// +// Output, int I4_POWER, the value of I^J. +// +{ + int k; + int value; + + if ( j < 0 ) + { + if ( i == 1 ) + { + value = 1; + } + else if ( i == 0 ) + { + cerr << "\n"; + cerr << "I4_POWER - Fatal error!\n"; + cerr << " I^J requested, with I = 0 and J negative.\n"; + exit ( 1 ); + } + else + { + value = 0; + } + } + else if ( j == 0 ) + { + if ( i == 0 ) + { + cerr << "\n"; + cerr << "I4_POWER - Fatal error!\n"; + cerr << " I^J requested, with I = 0 and J = 0.\n"; + exit ( 1 ); + } + else + { + value = 1; + } + } + else if ( j == 1 ) + { + value = i; + } + else + { + value = 1; + for ( k = 1; k <= j; k++ ) + { + value = value * i; + } + } + return value; +} +//****************************************************************************80 + +int i4_wrap ( int ival, int ilo, int ihi ) + +//****************************************************************************80 +// +// Purpose: +// +// I4_WRAP forces an I4 to lie between given limits by wrapping. +// +// Example: +// +// ILO = 4, IHI = 8 +// +// I Value +// +// -2 8 +// -1 4 +// 0 5 +// 1 6 +// 2 7 +// 3 8 +// 4 4 +// 5 5 +// 6 6 +// 7 7 +// 8 8 +// 9 4 +// 10 5 +// 11 6 +// 12 7 +// 13 8 +// 14 4 +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 19 August 2003 +// +// Author: +// +// John Burkardt +// +// Parameters: +// +// Input, int IVAL, an integer value. +// +// Input, int ILO, IHI, the desired bounds for the integer value. +// +// Output, int I4_WRAP, a "wrapped" version of IVAL. +// +{ + int jhi; + int jlo; + int value; + int wide; + + jlo = i4_min ( ilo, ihi ); + jhi = i4_max ( ilo, ihi ); + + wide = jhi + 1 - jlo; + + if ( wide == 1 ) + { + value = jlo; + } + else + { + value = jlo + i4_modp ( ival - jlo, wide ); + } + + return value; +} +//****************************************************************************80 + +void r8vec_copy (int n, double a1[], csnd::Vector<MYFLT> a2) + +//****************************************************************************80 +// +// Purpose: +// +// R8VEC_COPY copies an R8VEC. +// +// Discussion: +// +// An R8VEC is a vector of R8's. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 03 July 2005 +// +// Author: +// +// John Burkardt +// +// Parameters: +// +// Input, int N, the number of entries in the vectors. +// +// Input, double A1[N], the vector to be copied. +// +// Output, double A2[N], the copy of A1. +// +{ + int i; + + for ( i = 0; i < n; i++ ) + { + a2[i] = a1[i]; + } + return; +} + +void r8vec_copy ( int n, csnd::Vector<MYFLT> a1, double a2[]) +{ + int i; + + for ( i = 0; i < n; i++ ) + { + a2[i] = a1[i]; + } + return; +} + +//****************************************************************************80 + +double *r8vec_copy_new (csnd::Csound* csound, int n, double a1[] ) + +//****************************************************************************80 +// +// Purpose: +// +// R8VEC_COPY_NEW copies an R8VEC to a "new" R8VEC. +// +// Discussion: +// +// An R8VEC is a vector of R8's. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 03 July 2008 +// +// Author: +// +// John Burkardt +// +// Parameters: +// +// Input, int N, the number of entries in the vectors. +// +// Input, double A1[N], the vector to be copied. +// +// Output, double R8VEC_COPY_NEW[N], the copy of A1. +// +{ + double *a2; + int i; + + a2 = (double*) csound->malloc(sizeof(double) * n); //new double[n]; + + for ( i = 0; i < n; i++ ) + { + a2[i] = a1[i]; + } + return a2; +} +//****************************************************************************80 + +//****************************************************************************80 + +void walsh (csnd::Csound* csound, csnd::Vector<MYFLT> x) + +//****************************************************************************80 +// +// Purpose: +// +// WALSH performs a fast Walsh transform. +// +// Discussion: +// +// This routine performs a fast Wash transform on an input series X +// leaving the transformed results in X. The array Y is working space. +// X and Y are dimensioned N, which must be a power of 2. +// The results of this Walsh transform are in sequency order. +// +// The output sequence could be normalized by dividing by N. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 16 March 2011 +// +// Author: +// +// Ken Beauchamp +// +// Reference: +// +// Ken Beauchamp, +// Walsh functions and their applications, +// Academic Press, 1975, +// ISBN: 0-12-084050-2, +// LC: QA404.5.B33. +// +// Parameters: +// +// Input, int N, the number of items in X. +// N must be a power of 2. +// +// Input/output, double X[N], the data to be transformed. +// +{ + double a; + int i; + int i1; + int is; + int j; + int j1; + int l; + int m; + int n1; + int n2; + double w; + double *y; + double z; + int n = x.len(); + + n2 = n / 2; + y = (double*) csound->malloc(sizeof(double) * n2); //new double[n2]; + m = i4_log_2 ( n ); + z = - 1.0; + + for ( j = 1; j <= m; j++ ) + { + n1 = i4_power ( 2, m - j + 1 ); + j1 = i4_power ( 2, j - 1 ); + for ( l = 1; l <= j1; l++ ) + { + is = ( l - 1 ) * n1 + 1; + i1 = 0; + w = z; + for ( i = is; i <= is + n1 - 1; i = i + 2 ) + { + a = x[i-1]; + x[is+i1-1] = a + x[i]; + i1 = i1 + 1; + y[i1-1] = ( x[i] - a ) * w; + w = w * z; + } + for ( i = 1; i <= n1 / 2; i++ ) + { + x[n1/2+is+i-2] = y[i-1]; + } + } + } + + csound->free(y); + + return; +} diff --git a/src/walsh2.cpp b/src/walsh2.cpp new file mode 100644 index 0000000..3e15592 --- /dev/null +++ b/src/walsh2.cpp @@ -0,0 +1,58 @@ +/* + walsh2.cpp + Copyright (C) 2018 Mottakin Chowdhury + Copyright (C) 2019 Richard Knight + + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser 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 + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with this program; if not, write to the Free Software Foundation, + Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + */ +#include <plugin.h> + +const int N = 1<<16; + + +void fwtanal(csnd::Vector<MYFLT> io) { + int n = io.len(); + for (int d = 1; d < n; d <<= 1) { + for (int i = 0, m = d<<1; i < n; i += m) { + for (int j = 0; j < d; j++) { /// Don't forget modulo if required + MYFLT x = io[i+j], y = io[i+j+d]; + io[i+j] = (x+y), io[i+j+d] = (x-y); // xor + // io[i+j] = x+y; // and + // io[i+j+d] = x+y; // or + } + } + } +} + + +void fwtsynth(csnd::Vector<MYFLT> io) { + int n = io.len(); + for (int d = 1; d < n; d <<= 1) { + for (int i = 0, m = d<<1; i < n; i += m) { + for (int j = 0; j < d; j++) { /// Don't forget modulo if required + MYFLT x = io[i+j]; + MYFLT y = io[i+j+d]; + /// Modular inverse if required here + ///////io[i+j] = (x+y)>>1, io[i+j+d] = (x-y)>>1; // xor + io[i+j] = (x+y)/2; + io[i+j+d] = (x-y)/2; + // io[i+j] = x-y; // and + // io[i+j+d] = y-x; // or + } + } + } +}
\ No newline at end of file diff --git a/src/windowing.cpp b/src/windowing.cpp new file mode 100644 index 0000000..bdc7fdf --- /dev/null +++ b/src/windowing.cpp @@ -0,0 +1,208 @@ +#include <plugin.h> + +typedef struct _fft { + OPDS h; + ARRAYDAT *out; + ARRAYDAT *in, *in2; + MYFLT *f; + MYFLT b; + int32_t n; + void *setup; + AUXCH mem; +} FFT; + + +static inline void tabinit(CSOUND *csound, ARRAYDAT *p, int size) +{ + size_t ss; + if (p->dimensions==0) { + p->dimensions = 1; + p->sizes = (int32_t*)csound->Malloc(csound, sizeof(int32_t)); + } + if (p->data == NULL) { + CS_VARIABLE* var = p->arrayType->createVariable(csound, NULL); + p->arrayMemberSize = var->memBlockSize; + ss = p->arrayMemberSize*size; + p->data = (MYFLT*)csound->Calloc(csound, ss); + p->allocated = ss; + } else if( (ss = p->arrayMemberSize*size) > p->allocated) { + p->data = (MYFLT*) csound->ReAlloc(csound, p->data, ss); + p->allocated = ss; + } + if (p->dimensions==1) p->sizes[0] = size; + //p->dimensions = 1; +} + +int32_t shiftin_init(CSOUND *csound, FFT *p) { + int32_t sizs = CS_KSMPS; + if(p->out->sizes[0] < sizs) + tabinit(csound, p->out, sizs); + p->n = 0; + return OK; +} + +int32_t shiftin_perf(CSOUND *csound, FFT *p) { + IGN(csound); + uint32_t siz = p->out->sizes[0], n = p->n; + MYFLT *in = ((MYFLT *) p->in); + if (n + CS_KSMPS < siz) { + memcpy(p->out->data+n,in,CS_KSMPS*sizeof(MYFLT)); + } + else { + int32_t num = siz - n; + memcpy(p->out->data+n,in,num*sizeof(MYFLT)); + memcpy(p->out->data,in+num,(CS_KSMPS-num)*sizeof(MYFLT)); + } + p->n = (n + CS_KSMPS)%siz; + return OK; +} + + +int32_t shiftout_init(CSOUND *csound, FFT *p) { + int32_t siz = p->in->sizes[0]; + p->n = ((int32_t)*((MYFLT *)p->in2) % siz); + if (UNLIKELY((uint32_t) siz < CS_KSMPS)) + return csound->InitError(csound, "%s", Str("input array too small\n")); + return OK; +} + +int32_t shiftout_perf(CSOUND *csound, FFT *p) { + IGN(csound); + uint32_t siz = p->in->sizes[0], n = p->n; + MYFLT *out = ((MYFLT *) p->out); + + if (n + CS_KSMPS < siz) { + memcpy(out,p->in->data+n,CS_KSMPS*sizeof(MYFLT)); + } + else { + int32_t num = siz - n; + memcpy(out,p->in->data+n,num*sizeof(MYFLT)); + memcpy(out+num,p->in->data,(CS_KSMPS-num)*sizeof(MYFLT)); + } + p->n = (n + CS_KSMPS)%siz; + return OK; +} + + +struct windowing : csnd::Plugin<1, 1> { + static constexpr char const *otypes = "k[]"; + static constexpr char const *itypes = "k[]"; + FWT<MYFLT> fwt; + + int init() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + out.init(csound, in.len()); + return OK; + } + + int kperf() { + csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0); + csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0); + std::copy(in.begin(), in.end(), out.begin()); + fwt.fwtanal(out); + return OK; + } +}; + +int32_t init_window(CSOUND *csound, FFT *p) { + int32_t N = p->in->sizes[0]; + int32_t i,type = (int32_t) *p->f; + MYFLT *w; + tabinit(csound, p->out, N); + if (p->mem.auxp == 0 || p->mem.size < N*sizeof(MYFLT)) + csound->AuxAlloc(csound, N*sizeof(MYFLT), &p->mem); + w = (MYFLT *) p->mem.auxp; + switch(type) { + case 0: + for (i=0; i<N; i++) w[i] = 0.54 - 0.46*cos(i*TWOPI/N); + break; + case 1: + default: + for (i=0; i<N; i++) w[i] = 0.5 - 0.5*cos(i*TWOPI/N); + } + return OK; +} + +int32_t perf_window(CSOUND *csound, FFT *p) { + IGN(csound); + int32_t i,end = p->out->sizes[0], off = *((MYFLT *)p->in2); + MYFLT *in, *out, *w; + in = p->in->data; + out = p->out->data; + w = (MYFLT *) p->mem.auxp; + /*while (off < 0) off += end; + for (i=0;i<end;i++) + out[(i+off)%end] = in[i]*w[i];*/ + if(off) off = end - off; + for(i=0;i<end;i++) + out[i] = in[i]*w[(i+off)%end]; + return OK; + +} + + + +static inline void tabensure2D(CSOUND *csound, ARRAYDAT *p, + int32_t rows, int32_t columns) +{ + if (p->data==NULL || p->dimensions == 0 || + (p->dimensions==2 && (p->sizes[0] < rows || p->sizes[1] < columns))) { + size_t ss; + if (p->data == NULL) { + CS_VARIABLE* var = p->arrayType->createVariable(csound, NULL); + p->arrayMemberSize = var->memBlockSize; + } + ss = p->arrayMemberSize*rows*columns; + if (p->data==NULL) { + p->data = (MYFLT*)csound->Calloc(csound, ss); + p->dimensions = 2; + p->sizes = (int32_t*)csound->Malloc(csound, sizeof(int32_t)*2); + } + else p->data = (MYFLT*) csound->ReAlloc(csound, p->data, ss); + p->sizes[0] = rows; p->sizes[1] = columns; + } +} +// setrow +int32_t set_rows_init(CSOUND *csound, FFT *p) { + int32_t sizs = p->in->sizes[0]; + int32_t row = *((MYFLT *)p->in2); + tabensure2D(csound, p->out, row+1, sizs); + return OK; +} + +int32_t set_rows_perf(CSOUND *csound, FFT *p) { + int32_t start = *((MYFLT *)p->in2); + if (UNLIKELY(start < 0 || start >= p->out->sizes[0])) + return csound->PerfError(csound, &(p->h), + Str("Error: index out of range\n")); + int32_t bytes = p->in->sizes[0]*sizeof(MYFLT); + start *= p->out->sizes[1]; + memcpy(p->out->data+start,p->in->data,bytes); + return OK; +} + + +// getrow +int32_t rows_init(CSOUND *csound, FFT *p) { + if (p->in->dimensions == 2) { + int32_t siz = p->in->sizes[1]; + tabinit(csound, p->out, siz); + return OK; + } + else + return csound->InitError(csound, "%s", + Str("in array not 2-dimensional\n")); +} + +int32_t rows_perf(CSOUND *csound, FFT *p) { + int32_t start = *((MYFLT *)p->in2); + if (LIKELY(start < p->in->sizes[0])) { + int32_t bytes = p->in->sizes[1]*sizeof(MYFLT); + start *= p->in->sizes[1]; + memcpy(p->out->data,p->in->data+start,bytes); + return OK; + } + else return csound->PerfError(csound, &(p->h), + Str("requested row is out of range\n")); +}
\ No newline at end of file |