Mercurial > hg > dcase2013_sc_rnh
changeset 0:1b2ac32f7152
import code sent by Gerard Roma
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/COPYING Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,661 @@ + GNU AFFERO GENERAL PUBLIC LICENSE + Version 3, 19 November 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. <http://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 <http://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 +<http://www.gnu.org/licenses/>.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/README.md Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,23 @@ +The files in this folder represent our submission for the Scene Classificatoin task of the IEEE D-CASE AASP Challenge (http://c4dm.eecs.qmul.ac.uk/sceneseventschallenge/) + +The code has been tested mainly on Matlab2012 on OSX +Required libraries: rastamat and libsvm + +The implemented approach is described in: +G. Roma, W. Nogueira, P.Herrera, _Recurrence Quantification Analysis Features for Environmental Sound Recognition_. Proceedings of IEEE Workshop on Applications of Signal Processing to Audio and Acoustics. New Paltz, USA 2013. + +The main files are: + +* classify_scenes.m +* analyze_files.m +* RQA.m + +The rest of matlab files can be used to test the code. Some of the files are taken from: https://soundsoftware.ac.uk/projects/aasp-d-case-metrics + +Two submissions were sent to the challenge, one uses hardcoded SVM parameters, the other does grid search: + +classify_scenes(tmp_path, train_file,test_file, output_file, 0) % use hardcoded parameters for SVM +classify_scenes(tmp_path, train_file,test_file, output_file, 1) % use grid search + +The temp_path is used to store features. +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/RQA.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,87 @@ +% Copyright 2013 MUSIC TECHNOLOGY GROUP, UNIVERSITAT POMPEU FABRA +% +% Written by Gerard Roma <gerard.roma@upf.edu> +% +% 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 <http://www.gnu.org/licenses/>. + +function [ftr] = RQA(D, lmin,vmin) + N = max(size(D)); + pts = sum(D(:)); + + det = 0; + diags=[]; + for i = 1:N-1 + line = [0;diag(D,i);0]; + starts=find(diff(line)==1); % line starts + ends=find(diff(line)==-1); % line ends + diags = [diags; ends-starts]; + end + + verts = []; + intervals = []; + for i = 1:N + line = [0;D(:,i);0]; + starts=find(diff(line)==1); % line starts + ends=find(diff(line)==-1); % line ends + verts = [verts;ends-starts]; + end + + rr = pts/N^2; % recurrence rate + det = sum(diags(diags>=lmin))/pts; % determinism + adll = mean(diags(diags>=lmin)); % average diagonal line length + sddll = std(diags(diags>=lmin)); % standard deviation + + dratio = N^2 * sum(diags(diags>=lmin))/(pts^2); + + lam = sum(verts(verts>=vmin))/pts; % laminarity + tt = mean(verts(verts>=vmin)); % trappnig time + sdvll = std(verts(verts>=lmin)); % standard deviation + if(isnan(tt)) + tt = 0; + end + + vratio = N^2 * sum(verts(verts>=vmin))/(pts^2); % ratio + + lmax = max(1,max(diags)); % maximum diagonal line length + vmax = max(1,max(verts)); % maximum vertical line length + + if(isempty(lmax)) + lmax = 1; + end + + + if(isempty(vmax)) + vmax = 1; + end + + ddiv = 1/lmax; % divergence + vdiv = 1/vmax; + + diagsH = hist(diags,max(diags)); + diagsH = diagsH/sum(diagsH); + z = find(diagsH==0); + temp = diagsH.*log(diagsH); + temp(z)=0; % remove NaNs from zero prob + dentropy = sum(temp); % shannon entropy (diagonal) + + + vertsH = hist(verts,max(verts)); + vertsH = vertsH/sum(vertsH); + z = find(vertsH==0); + temp = vertsH.*log(vertsH); + temp(z)=0; % remove NaNs from zero prob + ventropy = sum(temp); % shannon entropy (vertical) + + ftr = [rr,det,adll,lam,tt,dratio,ddiv,dentropy,vratio,ventropy,vdiv]; + end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/analyze_files.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,61 @@ +% Copyright 2013 MUSIC TECHNOLOGY GROUP, UNIVERSITAT POMPEU FABRA +% +% Written by Gerard Roma <gerard.roma@upf.edu> +% +% 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 <http://www.gnu.org/licenses/>. + +function [ features ] = analyze_files( input_files, tmp_path ) + use_stored_features = 1; + + features = []; + for i=1:length(input_files) + disp(strcat('analyzing',input_files(i))); + features = [features; analyze(input_files(i))]; + end + + function ftrs = analyze(in_file) + w = 40; + h = 20; + r = 0.03; + dl = 2; + vl = 2; + [pathstr, name, ext] = fileparts(char(in_file)); + tmp_fname = strcat(tmp_path,'/',name,'.mat'); + if use_stored_features && (exist(tmp_fname)==2) + load(tmp_fname); + else + [x,fs] = wavread(char(in_file)); + x = x(:,1); + mfcc = melfcc(x,fs,'minfreq',0,'maxfreq',900,'dither',1)'; + N = max(size(mfcc)); + steps = round(N/h); + wrqa = []; + for i = 0:steps-1 + init = (i*h)+1; + final = min((i*h)+w,N); + ftr = mfcc(init:final,:); + d = squareform(pdist(ftr,'cosine')); + D = d<r; + rqa = RQA(D,dl,vl); + wrqa = [wrqa;rqa]; + end + + ftrs = [mean(mfcc,1),std(mfcc,1),mean(wrqa,1)]; + ftrs(isnan(ftrs))=0; + ftrs(isinf(ftrs))=0; + save(tmp_fname,'ftrs'); + end + end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/classify_scenes.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,107 @@ +% Copyright 2013 MUSIC TECHNOLOGY GROUP, UNIVERSITAT POMPEU FABRA +% +% Written by Gerard Roma <gerard.roma@upf.edu> +% +% 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 <http://www.gnu.org/licenses/>. + +function [ predict_label ] = classify_scenes(tmp_path, train_file,test_file, output_file, grid_search) + addpath('rastamat/'); + addpath('libsvm/') + + classes = {'bus' 'busystreet' 'office' 'openairmarket' 'park' 'quietstreet' 'restaurant' 'supermarket' 'tube' 'tubestation'}; + + [tr_names,tr_labels] = loadClassificationOutput(train_file); + [te_names,te_labels] = loadClassificationOutput(test_file); + + + %% extract features + disp('analyzing training files'); + train_features = analyze_files(tr_names,tmp_path); + disp('analyzing test files'); + test_features = analyze_files(te_names, tmp_path); + train_labels = []; + test_labels = ones(length(te_names),1); + + train_labels = get_class_indices(tr_labels); + + [train_z, mu, sigma] = zscore(train_features); + MU = repmat(mu,size(test_features,1),1); + SIGMA = repmat(sigma,size(test_features,1),1); + test_z = (test_features -MU)./SIGMA; + + + %% SVM grid Search inspired in http://labrosa.ee.columbia.edu/projects/consumervideo/ + + if grid_search + disp('performing grid search'); + tuning_data_fraction = 0.9; + train_size = round(tuning_data_fraction*size(train_z,1)); + p = randperm(size(train_z,1)); + trainX = train_z(p(1:train_size),:); + validateX = train_z(p(train_size+1:end),:); + trainY = get_class_indices(tr_labels(p(1:train_size))); + validateY = get_class_indices(tr_labels(p(train_size+1:end))); + + + gamma = 2.^[-10:1:10]; + C = 2.^[0:13]; + + params = zeros(length(gamma),length(C)); + best_a = 0; + best_g = 0; + best_C = 0; + + for gi= 1:length(gamma) + for Ci = 1:length(C) + m = svmtrain(trainY', trainX, sprintf('-c %d -g %2.5f -q',C(Ci),gamma(gi))); + [p,a] = svmpredict(validateY', validateX, m, '-q'); + if a(1) > best_a + best_C = C(Ci); + best_g = gamma(gi); + best_a = a(1); + end + end + end + disp('grid search done'); + else + best_C=70; + best_g = 0.003; + end + + + model = svmtrain(train_labels', train_z, sprintf('-c %d -g %2.5f -q',best_C,best_g)); + [predict_indices, accuracy_obj, prob_values] = svmpredict(test_labels, test_z, model,'-q'); + predict_label = classes(predict_indices); + + + outfd = fopen(output_file,'w+'); + + for i = 1:length(te_names) + fprintf(outfd,'%s\t',[char(te_names(i))]); + fprintf(outfd,'%s\n',[char(predict_label(i))]); + end + fclose(outfd); + + function idx = get_class_indices(labels) + idx = []; + for i=1:length(labels) + class = find(strcmp(classes,labels(i))); + idx(i) = class; + end + +end + + +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/foldGTlist.txt Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,5 @@ +fold1_gt.txt +fold2_gt.txt +fold3_gt.txt +fold4_gt.txt +fold5_gt.txt
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/foldRESlist.txt Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,5 @@ +fold1_result.txt +fold2_result.txt +fold3_result.txt +fold4_result.txt +fold5_result.txt
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/get_filenames.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,37 @@ +% Copyright 2013 MUSIC TECHNOLOGY GROUP, UNIVERSITAT POMPEU FABRA +% +% Written by Gerard Roma <gerard.roma@upf.edu> +% +% 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 <http://www.gnu.org/licenses/>. + +function [x,y] = get_filenames(dir) + classes = {'bus' 'busystreet' 'office' 'openairmarket' 'park' 'quietstreet' 'restaurant' 'supermarket' 'tube' 'tubestation'}; + + + names = []; + labels = []; + + for i=1:length(classes) + for j = 1:10 + suffix=sprintf('%02d.wav',j); + name = strcat(dir,classes(i),suffix); + names = [names;name]; + labels = [labels; i]; + end + end + x = names; + y = labels; +end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libsvm_linux64/Makefile Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,53 @@ +# This Makefile is used under Linux + +MATLABDIR ?= /usr/local/matlab +# for Mac +# MATLABDIR ?= /opt/local/matlab + +CXX ?= g++ +#CXX = g++-4.1 +CFLAGS = -Wall -Wconversion -O3 -fPIC -I$(MATLABDIR)/extern/include -I.. + +MEX = $(MATLABDIR)/bin/mex +MEX_OPTION = CC\#$(CXX) CXX\#$(CXX) CFLAGS\#"$(CFLAGS)" CXXFLAGS\#"$(CFLAGS)" +# comment the following line if you use MATLAB on 32-bit computer +MEX_OPTION += -largeArrayDims +MEX_EXT = $(shell $(MATLABDIR)/bin/mexext) + +OCTAVEDIR ?= /usr/include/octave +OCTAVE_MEX = env CC=$(CXX) mkoctfile +OCTAVE_MEX_OPTION = --mex +OCTAVE_MEX_EXT = mex +OCTAVE_CFLAGS = -Wall -O3 -fPIC -I$(OCTAVEDIR) -I.. + +all: matlab + +matlab: binary + +octave: + @make MEX="$(OCTAVE_MEX)" MEX_OPTION="$(OCTAVE_MEX_OPTION)" \ + MEX_EXT="$(OCTAVE_MEX_EXT)" CFLAGS="$(OCTAVE_CFLAGS)" \ + binary + +binary: svmpredict.$(MEX_EXT) svmtrain.$(MEX_EXT) libsvmread.$(MEX_EXT) libsvmwrite.$(MEX_EXT) + +svmpredict.$(MEX_EXT): svmpredict.c ../svm.h ../svm.o svm_model_matlab.o + $(MEX) $(MEX_OPTION) svmpredict.c ../svm.o svm_model_matlab.o + +svmtrain.$(MEX_EXT): svmtrain.c ../svm.h ../svm.o svm_model_matlab.o + $(MEX) $(MEX_OPTION) svmtrain.c ../svm.o svm_model_matlab.o + +libsvmread.$(MEX_EXT): libsvmread.c + $(MEX) $(MEX_OPTION) libsvmread.c + +libsvmwrite.$(MEX_EXT): libsvmwrite.c + $(MEX) $(MEX_OPTION) libsvmwrite.c + +svm_model_matlab.o: svm_model_matlab.c ../svm.h + $(CXX) $(CFLAGS) -c svm_model_matlab.c + +../svm.o: ../svm.cpp ../svm.h + make -C .. svm.o + +clean: + rm -f *~ *.o *.mex* *.obj ../svm.o
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libsvm_linux64/README Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,244 @@ +----------------------------------------- +--- MATLAB/OCTAVE interface of LIBSVM --- +----------------------------------------- + +Table of Contents +================= + +- Introduction +- Installation +- Usage +- Returned Model Structure +- Other Utilities +- Examples +- Additional Information + + +Introduction +============ + +This tool provides a simple interface to LIBSVM, a library for support vector +machines (http://www.csie.ntu.edu.tw/~cjlin/libsvm). It is very easy to use as +the usage and the way of specifying parameters are the same as that of LIBSVM. + +Installation +============ + +On Windows systems, pre-built binary files are already in the +directory '..\windows', so no need to conduct installation. Now we +provide binary files only for 64bit MATLAB on Windows. If you would +like to re-build the package, please rely on the following steps. + +We recommend using make.m on both MATLAB and OCTAVE. Just type 'make' +to build 'libsvmread.mex', 'libsvmwrite.mex', 'svmtrain.mex', and +'svmpredict.mex'. + +On MATLAB or Octave: + + >> make + +If make.m does not work on MATLAB (especially for Windows), try 'mex +-setup' to choose a suitable compiler for mex. Make sure your compiler +is accessible and workable. Then type 'make' to start the +installation. + +Example: + + matlab>> mex -setup + (ps: MATLAB will show the following messages to setup default compiler.) + Please choose your compiler for building external interface (MEX) files: + Would you like mex to locate installed compilers [y]/n? y + Select a compiler: + [1] Microsoft Visual C/C++ version 7.1 in C:\Program Files\Microsoft Visual Studio + [0] None + Compiler: 1 + Please verify your choices: + Compiler: Microsoft Visual C/C++ 7.1 + Location: C:\Program Files\Microsoft Visual Studio + Are these correct?([y]/n): y + + matlab>> make + +On Unix systems, if neither make.m nor 'mex -setup' works, please use +Makefile and type 'make' in a command window. Note that we assume +your MATLAB is installed in '/usr/local/matlab'. If not, please change +MATLABDIR in Makefile. + +Example: + linux> make + +To use octave, type 'make octave': + +Example: + linux> make octave + +For a list of supported/compatible compilers for MATLAB, please check +the following page: + +http://www.mathworks.com/support/compilers/current_release/ + +Usage +===== + +matlab> model = svmtrain(training_label_vector, training_instance_matrix [, 'libsvm_options']); + + -training_label_vector: + An m by 1 vector of training labels (type must be double). + -training_instance_matrix: + An m by n matrix of m training instances with n features. + It can be dense or sparse (type must be double). + -libsvm_options: + A string of training options in the same format as that of LIBSVM. + +matlab> [predicted_label, accuracy, decision_values/prob_estimates] = svmpredict(testing_label_vector, testing_instance_matrix, model [, 'libsvm_options']); +matlab> [predicted_label] = svmpredict(testing_label_vector, testing_instance_matrix, model [, 'libsvm_options']); + + -testing_label_vector: + An m by 1 vector of prediction labels. If labels of test + data are unknown, simply use any random values. (type must be double) + -testing_instance_matrix: + An m by n matrix of m testing instances with n features. + It can be dense or sparse. (type must be double) + -model: + The output of svmtrain. + -libsvm_options: + A string of testing options in the same format as that of LIBSVM. + +Returned Model Structure +======================== + +The 'svmtrain' function returns a model which can be used for future +prediction. It is a structure and is organized as [Parameters, nr_class, +totalSV, rho, Label, ProbA, ProbB, nSV, sv_coef, SVs]: + + -Parameters: parameters + -nr_class: number of classes; = 2 for regression/one-class svm + -totalSV: total #SV + -rho: -b of the decision function(s) wx+b + -Label: label of each class; empty for regression/one-class SVM + -ProbA: pairwise probability information; empty if -b 0 or in one-class SVM + -ProbB: pairwise probability information; empty if -b 0 or in one-class SVM + -nSV: number of SVs for each class; empty for regression/one-class SVM + -sv_coef: coefficients for SVs in decision functions + -SVs: support vectors + +If you do not use the option '-b 1', ProbA and ProbB are empty +matrices. If the '-v' option is specified, cross validation is +conducted and the returned model is just a scalar: cross-validation +accuracy for classification and mean-squared error for regression. + +More details about this model can be found in LIBSVM FAQ +(http://www.csie.ntu.edu.tw/~cjlin/libsvm/faq.html) and LIBSVM +implementation document +(http://www.csie.ntu.edu.tw/~cjlin/papers/libsvm.pdf). + +Result of Prediction +==================== + +The function 'svmpredict' has three outputs. The first one, +predictd_label, is a vector of predicted labels. The second output, +accuracy, is a vector including accuracy (for classification), mean +squared error, and squared correlation coefficient (for regression). +The third is a matrix containing decision values or probability +estimates (if '-b 1' is specified). If k is the number of classes +in training data, for decision values, each row includes results of +predicting k(k-1)/2 binary-class SVMs. For classification, k = 1 is a +special case. Decision value +1 is returned for each testing instance, +instead of an empty vector. For probabilities, each row contains k values +indicating the probability that the testing instance is in each class. +Note that the order of classes here is the same as 'Label' field +in the model structure. + +Other Utilities +=============== + +A matlab function libsvmread reads files in LIBSVM format: + +[label_vector, instance_matrix] = libsvmread('data.txt'); + +Two outputs are labels and instances, which can then be used as inputs +of svmtrain or svmpredict. + +A matlab function libsvmwrite writes Matlab matrix to a file in LIBSVM format: + +libsvmwrite('data.txt', label_vector, instance_matrix] + +The instance_matrix must be a sparse matrix. (type must be double) +For 32bit and 64bit MATLAB on Windows, pre-built binary files are ready +in the directory `..\windows', but in future releases, we will only +include 64bit MATLAB binary files. + +These codes are prepared by Rong-En Fan and Kai-Wei Chang from National +Taiwan University. + +Examples +======== + +Train and test on the provided data heart_scale: + +matlab> [heart_scale_label, heart_scale_inst] = libsvmread('../heart_scale'); +matlab> model = svmtrain(heart_scale_label, heart_scale_inst, '-c 1 -g 0.07'); +matlab> [predict_label, accuracy, dec_values] = svmpredict(heart_scale_label, heart_scale_inst, model); % test the training data + +For probability estimates, you need '-b 1' for training and testing: + +matlab> [heart_scale_label, heart_scale_inst] = libsvmread('../heart_scale'); +matlab> model = svmtrain(heart_scale_label, heart_scale_inst, '-c 1 -g 0.07 -b 1'); +matlab> [heart_scale_label, heart_scale_inst] = libsvmread('../heart_scale'); +matlab> [predict_label, accuracy, prob_estimates] = svmpredict(heart_scale_label, heart_scale_inst, model, '-b 1'); + +To use precomputed kernel, you must include sample serial number as +the first column of the training and testing data (assume your kernel +matrix is K, # of instances is n): + +matlab> K1 = [(1:n)', K]; % include sample serial number as first column +matlab> model = svmtrain(label_vector, K1, '-t 4'); +matlab> [predict_label, accuracy, dec_values] = svmpredict(label_vector, K1, model); % test the training data + +We give the following detailed example by splitting heart_scale into +150 training and 120 testing data. Constructing a linear kernel +matrix and then using the precomputed kernel gives exactly the same +testing error as using the LIBSVM built-in linear kernel. + +matlab> [heart_scale_label, heart_scale_inst] = libsvmread('../heart_scale'); +matlab> +matlab> % Split Data +matlab> train_data = heart_scale_inst(1:150,:); +matlab> train_label = heart_scale_label(1:150,:); +matlab> test_data = heart_scale_inst(151:270,:); +matlab> test_label = heart_scale_label(151:270,:); +matlab> +matlab> % Linear Kernel +matlab> model_linear = svmtrain(train_label, train_data, '-t 0'); +matlab> [predict_label_L, accuracy_L, dec_values_L] = svmpredict(test_label, test_data, model_linear); +matlab> +matlab> % Precomputed Kernel +matlab> model_precomputed = svmtrain(train_label, [(1:150)', train_data*train_data'], '-t 4'); +matlab> [predict_label_P, accuracy_P, dec_values_P] = svmpredict(test_label, [(1:120)', test_data*train_data'], model_precomputed); +matlab> +matlab> accuracy_L % Display the accuracy using linear kernel +matlab> accuracy_P % Display the accuracy using precomputed kernel + +Note that for testing, you can put anything in the +testing_label_vector. For more details of precomputed kernels, please +read the section ``Precomputed Kernels'' in the README of the LIBSVM +package. + +Additional Information +====================== + +This interface was initially written by Jun-Cheng Chen, Kuan-Jen Peng, +Chih-Yuan Yang and Chih-Huai Cheng from Department of Computer +Science, National Taiwan University. The current version was prepared +by Rong-En Fan and Ting-Fan Wu. If you find this tool useful, please +cite LIBSVM as follows + +Chih-Chung Chang and Chih-Jen Lin, LIBSVM : a library for support +vector machines. ACM Transactions on Intelligent Systems and +Technology, 2:27:1--27:27, 2011. Software available at +http://www.csie.ntu.edu.tw/~cjlin/libsvm + +For any question, please contact Chih-Jen Lin <cjlin@csie.ntu.edu.tw>, +or check the FAQ page: + +http://www.csie.ntu.edu.tw/~cjlin/libsvm/faq.html#/Q9:_MATLAB_interface
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libsvm_linux64/libsvmread.c Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,213 @@ +#include <stdio.h> +#include <string.h> +#include <stdlib.h> +#include <ctype.h> +#include <errno.h> + +#include "mex.h" + +#ifdef MX_API_VER +#if MX_API_VER < 0x07030000 +typedef int mwIndex; +#endif +#endif +#ifndef max +#define max(x,y) (((x)>(y))?(x):(y)) +#endif +#ifndef min +#define min(x,y) (((x)<(y))?(x):(y)) +#endif + +void exit_with_help() +{ + mexPrintf( + "Usage: [label_vector, instance_matrix] = libsvmread('filename');\n" + ); +} + +static void fake_answer(int nlhs, mxArray *plhs[]) +{ + int i; + for(i=0;i<nlhs;i++) + plhs[i] = mxCreateDoubleMatrix(0, 0, mxREAL); +} + +static char *line; +static int max_line_len; + +static char* readline(FILE *input) +{ + int len; + + if(fgets(line,max_line_len,input) == NULL) + return NULL; + + while(strrchr(line,'\n') == NULL) + { + max_line_len *= 2; + line = (char *) realloc(line, max_line_len); + len = (int) strlen(line); + if(fgets(line+len,max_line_len-len,input) == NULL) + break; + } + return line; +} + +// read in a problem (in libsvm format) +void read_problem(const char *filename, int nlhs, mxArray *plhs[]) +{ + int max_index, min_index, inst_max_index, i; + long elements, k; + FILE *fp = fopen(filename,"r"); + int l = 0; + char *endptr; + mwIndex *ir, *jc; + double *labels, *samples; + + if(fp == NULL) + { + mexPrintf("can't open input file %s\n",filename); + fake_answer(nlhs, plhs); + return; + } + + max_line_len = 1024; + line = (char *) malloc(max_line_len*sizeof(char)); + + max_index = 0; + min_index = 1; // our index starts from 1 + elements = 0; + while(readline(fp) != NULL) + { + char *idx, *val; + // features + int index = 0; + + inst_max_index = -1; // strtol gives 0 if wrong format, and precomputed kernel has <index> start from 0 + strtok(line," \t"); // label + while (1) + { + idx = strtok(NULL,":"); // index:value + val = strtok(NULL," \t"); + if(val == NULL) + break; + + errno = 0; + index = (int) strtol(idx,&endptr,10); + if(endptr == idx || errno != 0 || *endptr != '\0' || index <= inst_max_index) + { + mexPrintf("Wrong input format at line %d\n",l+1); + fake_answer(nlhs, plhs); + return; + } + else + inst_max_index = index; + + min_index = min(min_index, index); + elements++; + } + max_index = max(max_index, inst_max_index); + l++; + } + rewind(fp); + + // y + plhs[0] = mxCreateDoubleMatrix(l, 1, mxREAL); + // x^T + if (min_index <= 0) + plhs[1] = mxCreateSparse(max_index-min_index+1, l, elements, mxREAL); + else + plhs[1] = mxCreateSparse(max_index, l, elements, mxREAL); + + labels = mxGetPr(plhs[0]); + samples = mxGetPr(plhs[1]); + ir = mxGetIr(plhs[1]); + jc = mxGetJc(plhs[1]); + + k=0; + for(i=0;i<l;i++) + { + char *idx, *val, *label; + jc[i] = k; + + readline(fp); + + label = strtok(line," \t\n"); + if(label == NULL) + { + mexPrintf("Empty line at line %d\n",i+1); + fake_answer(nlhs, plhs); + return; + } + labels[i] = strtod(label,&endptr); + if(endptr == label || *endptr != '\0') + { + mexPrintf("Wrong input format at line %d\n",i+1); + fake_answer(nlhs, plhs); + return; + } + + // features + while(1) + { + idx = strtok(NULL,":"); + val = strtok(NULL," \t"); + if(val == NULL) + break; + + ir[k] = (mwIndex) (strtol(idx,&endptr,10) - min_index); // precomputed kernel has <index> start from 0 + + errno = 0; + samples[k] = strtod(val,&endptr); + if (endptr == val || errno != 0 || (*endptr != '\0' && !isspace(*endptr))) + { + mexPrintf("Wrong input format at line %d\n",i+1); + fake_answer(nlhs, plhs); + return; + } + ++k; + } + } + jc[l] = k; + + fclose(fp); + free(line); + + { + mxArray *rhs[1], *lhs[1]; + rhs[0] = plhs[1]; + if(mexCallMATLAB(1, lhs, 1, rhs, "transpose")) + { + mexPrintf("Error: cannot transpose problem\n"); + fake_answer(nlhs, plhs); + return; + } + plhs[1] = lhs[0]; + } +} + +void mexFunction( int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] ) +{ + char filename[256]; + + if(nrhs != 1 || nlhs != 2) + { + exit_with_help(); + fake_answer(nlhs, plhs); + return; + } + + mxGetString(prhs[0], filename, mxGetN(prhs[0]) + 1); + + if(filename == NULL) + { + mexPrintf("Error: filename is NULL\n"); + return; + } + + read_problem(filename, nlhs, plhs); + + return; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libsvm_linux64/libsvmwrite.c Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,120 @@ +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include "mex.h" + +#ifdef MX_API_VER +#if MX_API_VER < 0x07030000 +typedef int mwIndex; +#endif +#endif + +void exit_with_help() +{ + mexPrintf( + "Usage: libsvmwrite('filename', label_vector, instance_matrix);\n" + ); +} + +static void fake_answer(int nlhs, mxArray *plhs[]) +{ + int i; + for(i=0;i<nlhs;i++) + plhs[i] = mxCreateDoubleMatrix(0, 0, mxREAL); +} + +void libsvmwrite(const char *filename, const mxArray *label_vec, const mxArray *instance_mat) +{ + FILE *fp = fopen(filename,"w"); + int i, k, low, high, l; + mwIndex *ir, *jc; + int label_vector_row_num; + double *samples, *labels; + mxArray *instance_mat_col; // instance sparse matrix in column format + + if(fp ==NULL) + { + mexPrintf("can't open output file %s\n",filename); + return; + } + + // transpose instance matrix + { + mxArray *prhs[1], *plhs[1]; + prhs[0] = mxDuplicateArray(instance_mat); + if(mexCallMATLAB(1, plhs, 1, prhs, "transpose")) + { + mexPrintf("Error: cannot transpose instance matrix\n"); + return; + } + instance_mat_col = plhs[0]; + mxDestroyArray(prhs[0]); + } + + // the number of instance + l = (int) mxGetN(instance_mat_col); + label_vector_row_num = (int) mxGetM(label_vec); + + if(label_vector_row_num!=l) + { + mexPrintf("Length of label vector does not match # of instances.\n"); + return; + } + + // each column is one instance + labels = mxGetPr(label_vec); + samples = mxGetPr(instance_mat_col); + ir = mxGetIr(instance_mat_col); + jc = mxGetJc(instance_mat_col); + + for(i=0;i<l;i++) + { + fprintf(fp,"%g", labels[i]); + + low = (int) jc[i], high = (int) jc[i+1]; + for(k=low;k<high;k++) + fprintf(fp," %ld:%g", ir[k]+1, samples[k]); + + fprintf(fp,"\n"); + } + + fclose(fp); + return; +} + +void mexFunction( int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] ) +{ + if(nlhs > 0) + { + exit_with_help(); + fake_answer(nlhs, plhs); + return; + } + + // Transform the input Matrix to libsvm format + if(nrhs == 3) + { + char filename[256]; + if(!mxIsDouble(prhs[1]) || !mxIsDouble(prhs[2])) + { + mexPrintf("Error: label vector and instance matrix must be double\n"); + return; + } + + mxGetString(prhs[0], filename, mxGetN(prhs[0])+1); + + if(mxIsSparse(prhs[2])) + libsvmwrite(filename, prhs[1], prhs[2]); + else + { + mexPrintf("Instance_matrix must be sparse\n"); + return; + } + } + else + { + exit_with_help(); + return; + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libsvm_linux64/make.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,21 @@ +% This make.m is for MATLAB and OCTAVE under Windows, Mac, and Unix + +try + Type = ver; + % This part is for OCTAVE + if(strcmp(Type(1).Name, 'Octave') == 1) + mex libsvmread.c + mex libsvmwrite.c + mex svmtrain.c ../svm.cpp svm_model_matlab.c + mex svmpredict.c ../svm.cpp svm_model_matlab.c + % This part is for MATLAB + % Add -largeArrayDims on 64-bit machines of MATLAB + else + mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims libsvmread.c + mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims libsvmwrite.c + mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims svmtrain.c ../svm.cpp svm_model_matlab.c + mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims svmpredict.c ../svm.cpp svm_model_matlab.c + end +catch + fprintf('If make.m fails, please check README about detailed instructions.\n'); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libsvm_linux64/svm_model_matlab.c Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,375 @@ +#include <stdlib.h> +#include <string.h> +#include "../svm.h" + +#include "mex.h" + +#ifdef MX_API_VER +#if MX_API_VER < 0x07030000 +typedef int mwIndex; +#endif +#endif + +#define NUM_OF_RETURN_FIELD 11 + +#define Malloc(type,n) (type *)malloc((n)*sizeof(type)) + +static const char *field_names[] = { + "Parameters", + "nr_class", + "totalSV", + "rho", + "Label", + "sv_indices", + "ProbA", + "ProbB", + "nSV", + "sv_coef", + "SVs" +}; + +const char *model_to_matlab_structure(mxArray *plhs[], int num_of_feature, struct svm_model *model) +{ + int i, j, n; + double *ptr; + mxArray *return_model, **rhs; + int out_id = 0; + + rhs = (mxArray **)mxMalloc(sizeof(mxArray *)*NUM_OF_RETURN_FIELD); + + // Parameters + rhs[out_id] = mxCreateDoubleMatrix(5, 1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + ptr[0] = model->param.svm_type; + ptr[1] = model->param.kernel_type; + ptr[2] = model->param.degree; + ptr[3] = model->param.gamma; + ptr[4] = model->param.coef0; + out_id++; + + // nr_class + rhs[out_id] = mxCreateDoubleMatrix(1, 1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + ptr[0] = model->nr_class; + out_id++; + + // total SV + rhs[out_id] = mxCreateDoubleMatrix(1, 1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + ptr[0] = model->l; + out_id++; + + // rho + n = model->nr_class*(model->nr_class-1)/2; + rhs[out_id] = mxCreateDoubleMatrix(n, 1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + for(i = 0; i < n; i++) + ptr[i] = model->rho[i]; + out_id++; + + // Label + if(model->label) + { + rhs[out_id] = mxCreateDoubleMatrix(model->nr_class, 1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + for(i = 0; i < model->nr_class; i++) + ptr[i] = model->label[i]; + } + else + rhs[out_id] = mxCreateDoubleMatrix(0, 0, mxREAL); + out_id++; + + // sv_indices + if(model->sv_indices) + { + rhs[out_id] = mxCreateDoubleMatrix(model->l, 1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + for(i = 0; i < model->l; i++) + ptr[i] = model->sv_indices[i]; + } + else + rhs[out_id] = mxCreateDoubleMatrix(0, 0, mxREAL); + out_id++; + + // probA + if(model->probA != NULL) + { + rhs[out_id] = mxCreateDoubleMatrix(n, 1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + for(i = 0; i < n; i++) + ptr[i] = model->probA[i]; + } + else + rhs[out_id] = mxCreateDoubleMatrix(0, 0, mxREAL); + out_id ++; + + // probB + if(model->probB != NULL) + { + rhs[out_id] = mxCreateDoubleMatrix(n, 1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + for(i = 0; i < n; i++) + ptr[i] = model->probB[i]; + } + else + rhs[out_id] = mxCreateDoubleMatrix(0, 0, mxREAL); + out_id++; + + // nSV + if(model->nSV) + { + rhs[out_id] = mxCreateDoubleMatrix(model->nr_class, 1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + for(i = 0; i < model->nr_class; i++) + ptr[i] = model->nSV[i]; + } + else + rhs[out_id] = mxCreateDoubleMatrix(0, 0, mxREAL); + out_id++; + + // sv_coef + rhs[out_id] = mxCreateDoubleMatrix(model->l, model->nr_class-1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + for(i = 0; i < model->nr_class-1; i++) + for(j = 0; j < model->l; j++) + ptr[(i*(model->l))+j] = model->sv_coef[i][j]; + out_id++; + + // SVs + { + int ir_index, nonzero_element; + mwIndex *ir, *jc; + mxArray *pprhs[1], *pplhs[1]; + + if(model->param.kernel_type == PRECOMPUTED) + { + nonzero_element = model->l; + num_of_feature = 1; + } + else + { + nonzero_element = 0; + for(i = 0; i < model->l; i++) { + j = 0; + while(model->SV[i][j].index != -1) + { + nonzero_element++; + j++; + } + } + } + + // SV in column, easier accessing + rhs[out_id] = mxCreateSparse(num_of_feature, model->l, nonzero_element, mxREAL); + ir = mxGetIr(rhs[out_id]); + jc = mxGetJc(rhs[out_id]); + ptr = mxGetPr(rhs[out_id]); + jc[0] = ir_index = 0; + for(i = 0;i < model->l; i++) + { + if(model->param.kernel_type == PRECOMPUTED) + { + // make a (1 x model->l) matrix + ir[ir_index] = 0; + ptr[ir_index] = model->SV[i][0].value; + ir_index++; + jc[i+1] = jc[i] + 1; + } + else + { + int x_index = 0; + while (model->SV[i][x_index].index != -1) + { + ir[ir_index] = model->SV[i][x_index].index - 1; + ptr[ir_index] = model->SV[i][x_index].value; + ir_index++, x_index++; + } + jc[i+1] = jc[i] + x_index; + } + } + // transpose back to SV in row + pprhs[0] = rhs[out_id]; + if(mexCallMATLAB(1, pplhs, 1, pprhs, "transpose")) + return "cannot transpose SV matrix"; + rhs[out_id] = pplhs[0]; + out_id++; + } + + /* Create a struct matrix contains NUM_OF_RETURN_FIELD fields */ + return_model = mxCreateStructMatrix(1, 1, NUM_OF_RETURN_FIELD, field_names); + + /* Fill struct matrix with input arguments */ + for(i = 0; i < NUM_OF_RETURN_FIELD; i++) + mxSetField(return_model,0,field_names[i],mxDuplicateArray(rhs[i])); + /* return */ + plhs[0] = return_model; + mxFree(rhs); + + return NULL; +} + +struct svm_model *matlab_matrix_to_model(const mxArray *matlab_struct, const char **msg) +{ + int i, j, n, num_of_fields; + double *ptr; + int id = 0; + struct svm_node *x_space; + struct svm_model *model; + mxArray **rhs; + + num_of_fields = mxGetNumberOfFields(matlab_struct); + if(num_of_fields != NUM_OF_RETURN_FIELD) + { + *msg = "number of return field is not correct"; + return NULL; + } + rhs = (mxArray **) mxMalloc(sizeof(mxArray *)*num_of_fields); + + for(i=0;i<num_of_fields;i++) + rhs[i] = mxGetFieldByNumber(matlab_struct, 0, i); + + model = Malloc(struct svm_model, 1); + model->rho = NULL; + model->probA = NULL; + model->probB = NULL; + model->label = NULL; + model->sv_indices = NULL; + model->nSV = NULL; + model->free_sv = 1; // XXX + + ptr = mxGetPr(rhs[id]); + model->param.svm_type = (int)ptr[0]; + model->param.kernel_type = (int)ptr[1]; + model->param.degree = (int)ptr[2]; + model->param.gamma = ptr[3]; + model->param.coef0 = ptr[4]; + id++; + + ptr = mxGetPr(rhs[id]); + model->nr_class = (int)ptr[0]; + id++; + + ptr = mxGetPr(rhs[id]); + model->l = (int)ptr[0]; + id++; + + // rho + n = model->nr_class * (model->nr_class-1)/2; + model->rho = (double*) malloc(n*sizeof(double)); + ptr = mxGetPr(rhs[id]); + for(i=0;i<n;i++) + model->rho[i] = ptr[i]; + id++; + + // label + if(mxIsEmpty(rhs[id]) == 0) + { + model->label = (int*) malloc(model->nr_class*sizeof(int)); + ptr = mxGetPr(rhs[id]); + for(i=0;i<model->nr_class;i++) + model->label[i] = (int)ptr[i]; + } + id++; + + // sv_indices + if(mxIsEmpty(rhs[id]) == 0) + { + model->sv_indices = (int*) malloc(model->l*sizeof(int)); + ptr = mxGetPr(rhs[id]); + for(i=0;i<model->l;i++) + model->sv_indices[i] = (int)ptr[i]; + } + id++; + + // probA + if(mxIsEmpty(rhs[id]) == 0) + { + model->probA = (double*) malloc(n*sizeof(double)); + ptr = mxGetPr(rhs[id]); + for(i=0;i<n;i++) + model->probA[i] = ptr[i]; + } + id++; + + // probB + if(mxIsEmpty(rhs[id]) == 0) + { + model->probB = (double*) malloc(n*sizeof(double)); + ptr = mxGetPr(rhs[id]); + for(i=0;i<n;i++) + model->probB[i] = ptr[i]; + } + id++; + + // nSV + if(mxIsEmpty(rhs[id]) == 0) + { + model->nSV = (int*) malloc(model->nr_class*sizeof(int)); + ptr = mxGetPr(rhs[id]); + for(i=0;i<model->nr_class;i++) + model->nSV[i] = (int)ptr[i]; + } + id++; + + // sv_coef + ptr = mxGetPr(rhs[id]); + model->sv_coef = (double**) malloc((model->nr_class-1)*sizeof(double)); + for( i=0 ; i< model->nr_class -1 ; i++ ) + model->sv_coef[i] = (double*) malloc((model->l)*sizeof(double)); + for(i = 0; i < model->nr_class - 1; i++) + for(j = 0; j < model->l; j++) + model->sv_coef[i][j] = ptr[i*(model->l)+j]; + id++; + + // SV + { + int sr, sc, elements; + int num_samples; + mwIndex *ir, *jc; + mxArray *pprhs[1], *pplhs[1]; + + // transpose SV + pprhs[0] = rhs[id]; + if(mexCallMATLAB(1, pplhs, 1, pprhs, "transpose")) + { + svm_free_and_destroy_model(&model); + *msg = "cannot transpose SV matrix"; + return NULL; + } + rhs[id] = pplhs[0]; + + sr = (int)mxGetN(rhs[id]); + sc = (int)mxGetM(rhs[id]); + + ptr = mxGetPr(rhs[id]); + ir = mxGetIr(rhs[id]); + jc = mxGetJc(rhs[id]); + + num_samples = (int)mxGetNzmax(rhs[id]); + + elements = num_samples + sr; + + model->SV = (struct svm_node **) malloc(sr * sizeof(struct svm_node *)); + x_space = (struct svm_node *)malloc(elements * sizeof(struct svm_node)); + + // SV is in column + for(i=0;i<sr;i++) + { + int low = (int)jc[i], high = (int)jc[i+1]; + int x_index = 0; + model->SV[i] = &x_space[low+i]; + for(j=low;j<high;j++) + { + model->SV[i][x_index].index = (int)ir[j] + 1; + model->SV[i][x_index].value = ptr[j]; + x_index++; + } + model->SV[i][x_index].index = -1; + } + + id++; + } + mxFree(rhs); + + return model; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libsvm_linux64/svm_model_matlab.h Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,2 @@ +const char *model_to_matlab_structure(mxArray *plhs[], int num_of_feature, struct svm_model *model); +struct svm_model *matlab_matrix_to_model(const mxArray *matlab_struct, const char **error_message);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libsvm_linux64/svmpredict.c Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,370 @@ +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include "../svm.h" + +#include "mex.h" +#include "svm_model_matlab.h" + +#ifdef MX_API_VER +#if MX_API_VER < 0x07030000 +typedef int mwIndex; +#endif +#endif + +#define CMD_LEN 2048 + +int print_null(const char *s,...) {} +int (*info)(const char *fmt,...) = &mexPrintf; + +void read_sparse_instance(const mxArray *prhs, int index, struct svm_node *x) +{ + int i, j, low, high; + mwIndex *ir, *jc; + double *samples; + + ir = mxGetIr(prhs); + jc = mxGetJc(prhs); + samples = mxGetPr(prhs); + + // each column is one instance + j = 0; + low = (int)jc[index], high = (int)jc[index+1]; + for(i=low;i<high;i++) + { + x[j].index = (int)ir[i] + 1; + x[j].value = samples[i]; + j++; + } + x[j].index = -1; +} + +static void fake_answer(int nlhs, mxArray *plhs[]) +{ + int i; + for(i=0;i<nlhs;i++) + plhs[i] = mxCreateDoubleMatrix(0, 0, mxREAL); +} + +void predict(int nlhs, mxArray *plhs[], const mxArray *prhs[], struct svm_model *model, const int predict_probability) +{ + int label_vector_row_num, label_vector_col_num; + int feature_number, testing_instance_number; + int instance_index; + double *ptr_instance, *ptr_label, *ptr_predict_label; + double *ptr_prob_estimates, *ptr_dec_values, *ptr; + struct svm_node *x; + mxArray *pplhs[1]; // transposed instance sparse matrix + mxArray *tplhs[3]; // temporary storage for plhs[] + + int correct = 0; + int total = 0; + double error = 0; + double sump = 0, sumt = 0, sumpp = 0, sumtt = 0, sumpt = 0; + + int svm_type=svm_get_svm_type(model); + int nr_class=svm_get_nr_class(model); + double *prob_estimates=NULL; + + // prhs[1] = testing instance matrix + feature_number = (int)mxGetN(prhs[1]); + testing_instance_number = (int)mxGetM(prhs[1]); + label_vector_row_num = (int)mxGetM(prhs[0]); + label_vector_col_num = (int)mxGetN(prhs[0]); + + if(label_vector_row_num!=testing_instance_number) + { + mexPrintf("Length of label vector does not match # of instances.\n"); + fake_answer(nlhs, plhs); + return; + } + if(label_vector_col_num!=1) + { + mexPrintf("label (1st argument) should be a vector (# of column is 1).\n"); + fake_answer(nlhs, plhs); + return; + } + + ptr_instance = mxGetPr(prhs[1]); + ptr_label = mxGetPr(prhs[0]); + + // transpose instance matrix + if(mxIsSparse(prhs[1])) + { + if(model->param.kernel_type == PRECOMPUTED) + { + // precomputed kernel requires dense matrix, so we make one + mxArray *rhs[1], *lhs[1]; + rhs[0] = mxDuplicateArray(prhs[1]); + if(mexCallMATLAB(1, lhs, 1, rhs, "full")) + { + mexPrintf("Error: cannot full testing instance matrix\n"); + fake_answer(nlhs, plhs); + return; + } + ptr_instance = mxGetPr(lhs[0]); + mxDestroyArray(rhs[0]); + } + else + { + mxArray *pprhs[1]; + pprhs[0] = mxDuplicateArray(prhs[1]); + if(mexCallMATLAB(1, pplhs, 1, pprhs, "transpose")) + { + mexPrintf("Error: cannot transpose testing instance matrix\n"); + fake_answer(nlhs, plhs); + return; + } + } + } + + if(predict_probability) + { + if(svm_type==NU_SVR || svm_type==EPSILON_SVR) + info("Prob. model for test data: target value = predicted value + z,\nz: Laplace distribution e^(-|z|/sigma)/(2sigma),sigma=%g\n",svm_get_svr_probability(model)); + else + prob_estimates = (double *) malloc(nr_class*sizeof(double)); + } + + tplhs[0] = mxCreateDoubleMatrix(testing_instance_number, 1, mxREAL); + if(predict_probability) + { + // prob estimates are in plhs[2] + if(svm_type==C_SVC || svm_type==NU_SVC) + tplhs[2] = mxCreateDoubleMatrix(testing_instance_number, nr_class, mxREAL); + else + tplhs[2] = mxCreateDoubleMatrix(0, 0, mxREAL); + } + else + { + // decision values are in plhs[2] + if(svm_type == ONE_CLASS || + svm_type == EPSILON_SVR || + svm_type == NU_SVR || + nr_class == 1) // if only one class in training data, decision values are still returned. + tplhs[2] = mxCreateDoubleMatrix(testing_instance_number, 1, mxREAL); + else + tplhs[2] = mxCreateDoubleMatrix(testing_instance_number, nr_class*(nr_class-1)/2, mxREAL); + } + + ptr_predict_label = mxGetPr(tplhs[0]); + ptr_prob_estimates = mxGetPr(tplhs[2]); + ptr_dec_values = mxGetPr(tplhs[2]); + x = (struct svm_node*)malloc((feature_number+1)*sizeof(struct svm_node) ); + for(instance_index=0;instance_index<testing_instance_number;instance_index++) + { + int i; + double target_label, predict_label; + + target_label = ptr_label[instance_index]; + + if(mxIsSparse(prhs[1]) && model->param.kernel_type != PRECOMPUTED) // prhs[1]^T is still sparse + read_sparse_instance(pplhs[0], instance_index, x); + else + { + for(i=0;i<feature_number;i++) + { + x[i].index = i+1; + x[i].value = ptr_instance[testing_instance_number*i+instance_index]; + } + x[feature_number].index = -1; + } + + if(predict_probability) + { + if(svm_type==C_SVC || svm_type==NU_SVC) + { + predict_label = svm_predict_probability(model, x, prob_estimates); + ptr_predict_label[instance_index] = predict_label; + for(i=0;i<nr_class;i++) + ptr_prob_estimates[instance_index + i * testing_instance_number] = prob_estimates[i]; + } else { + predict_label = svm_predict(model,x); + ptr_predict_label[instance_index] = predict_label; + } + } + else + { + if(svm_type == ONE_CLASS || + svm_type == EPSILON_SVR || + svm_type == NU_SVR) + { + double res; + predict_label = svm_predict_values(model, x, &res); + ptr_dec_values[instance_index] = res; + } + else + { + double *dec_values = (double *) malloc(sizeof(double) * nr_class*(nr_class-1)/2); + predict_label = svm_predict_values(model, x, dec_values); + if(nr_class == 1) + ptr_dec_values[instance_index] = 1; + else + for(i=0;i<(nr_class*(nr_class-1))/2;i++) + ptr_dec_values[instance_index + i * testing_instance_number] = dec_values[i]; + free(dec_values); + } + ptr_predict_label[instance_index] = predict_label; + } + + if(predict_label == target_label) + ++correct; + error += (predict_label-target_label)*(predict_label-target_label); + sump += predict_label; + sumt += target_label; + sumpp += predict_label*predict_label; + sumtt += target_label*target_label; + sumpt += predict_label*target_label; + ++total; + } + if(svm_type==NU_SVR || svm_type==EPSILON_SVR) + { + info("Mean squared error = %g (regression)\n",error/total); + info("Squared correlation coefficient = %g (regression)\n", + ((total*sumpt-sump*sumt)*(total*sumpt-sump*sumt))/ + ((total*sumpp-sump*sump)*(total*sumtt-sumt*sumt)) + ); + } + else + info("Accuracy = %g%% (%d/%d) (classification)\n", + (double)correct/total*100,correct,total); + + // return accuracy, mean squared error, squared correlation coefficient + tplhs[1] = mxCreateDoubleMatrix(3, 1, mxREAL); + ptr = mxGetPr(tplhs[1]); + ptr[0] = (double)correct/total*100; + ptr[1] = error/total; + ptr[2] = ((total*sumpt-sump*sumt)*(total*sumpt-sump*sumt))/ + ((total*sumpp-sump*sump)*(total*sumtt-sumt*sumt)); + + free(x); + if(prob_estimates != NULL) + free(prob_estimates); + + switch(nlhs) + { + case 3: + plhs[2] = tplhs[2]; + plhs[1] = tplhs[1]; + case 1: + case 0: + plhs[0] = tplhs[0]; + } +} + +void exit_with_help() +{ + mexPrintf( + "Usage: [predicted_label, accuracy, decision_values/prob_estimates] = svmpredict(testing_label_vector, testing_instance_matrix, model, 'libsvm_options')\n" + " [predicted_label] = svmpredict(testing_label_vector, testing_instance_matrix, model, 'libsvm_options')\n" + "Parameters:\n" + " model: SVM model structure from svmtrain.\n" + " libsvm_options:\n" + " -b probability_estimates: whether to predict probability estimates, 0 or 1 (default 0); one-class SVM not supported yet\n" + " -q : quiet mode (no outputs)\n" + "Returns:\n" + " predicted_label: SVM prediction output vector.\n" + " accuracy: a vector with accuracy, mean squared error, squared correlation coefficient.\n" + " prob_estimates: If selected, probability estimate vector.\n" + ); +} + +void mexFunction( int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] ) +{ + int prob_estimate_flag = 0; + struct svm_model *model; + info = &mexPrintf; + + if(nlhs == 2 || nlhs > 3 || nrhs > 4 || nrhs < 3) + { + exit_with_help(); + fake_answer(nlhs, plhs); + return; + } + + if(!mxIsDouble(prhs[0]) || !mxIsDouble(prhs[1])) { + mexPrintf("Error: label vector and instance matrix must be double\n"); + fake_answer(nlhs, plhs); + return; + } + + if(mxIsStruct(prhs[2])) + { + const char *error_msg; + + // parse options + if(nrhs==4) + { + int i, argc = 1; + char cmd[CMD_LEN], *argv[CMD_LEN/2]; + + // put options in argv[] + mxGetString(prhs[3], cmd, mxGetN(prhs[3]) + 1); + if((argv[argc] = strtok(cmd, " ")) != NULL) + while((argv[++argc] = strtok(NULL, " ")) != NULL) + ; + + for(i=1;i<argc;i++) + { + if(argv[i][0] != '-') break; + if((++i>=argc) && argv[i-1][1] != 'q') + { + exit_with_help(); + fake_answer(nlhs, plhs); + return; + } + switch(argv[i-1][1]) + { + case 'b': + prob_estimate_flag = atoi(argv[i]); + break; + case 'q': + i--; + info = &print_null; + break; + default: + mexPrintf("Unknown option: -%c\n", argv[i-1][1]); + exit_with_help(); + fake_answer(nlhs, plhs); + return; + } + } + } + + model = matlab_matrix_to_model(prhs[2], &error_msg); + if (model == NULL) + { + mexPrintf("Error: can't read model: %s\n", error_msg); + fake_answer(nlhs, plhs); + return; + } + + if(prob_estimate_flag) + { + if(svm_check_probability_model(model)==0) + { + mexPrintf("Model does not support probabiliy estimates\n"); + fake_answer(nlhs, plhs); + svm_free_and_destroy_model(&model); + return; + } + } + else + { + if(svm_check_probability_model(model)!=0) + info("Model supports probability estimates, but disabled in predicton.\n"); + } + + predict(nlhs, plhs, prhs, model, prob_estimate_flag); + // destroy model + svm_free_and_destroy_model(&model); + } + else + { + mexPrintf("model file should be a struct array\n"); + fake_answer(nlhs, plhs); + } + + return; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libsvm_linux64/svmtrain.c Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,483 @@ +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <ctype.h> +#include "../svm.h" + +#include "mex.h" +#include "svm_model_matlab.h" + +#ifdef MX_API_VER +#if MX_API_VER < 0x07030000 +typedef int mwIndex; +#endif +#endif + +#define CMD_LEN 2048 +#define Malloc(type,n) (type *)malloc((n)*sizeof(type)) + +void print_null(const char *s) {} +void print_string_matlab(const char *s) {mexPrintf(s);} + +void exit_with_help() +{ + mexPrintf( + "Usage: model = svmtrain(training_label_vector, training_instance_matrix, 'libsvm_options');\n" + "libsvm_options:\n" + "-s svm_type : set type of SVM (default 0)\n" + " 0 -- C-SVC (multi-class classification)\n" + " 1 -- nu-SVC (multi-class classification)\n" + " 2 -- one-class SVM\n" + " 3 -- epsilon-SVR (regression)\n" + " 4 -- nu-SVR (regression)\n" + "-t kernel_type : set type of kernel function (default 2)\n" + " 0 -- linear: u'*v\n" + " 1 -- polynomial: (gamma*u'*v + coef0)^degree\n" + " 2 -- radial basis function: exp(-gamma*|u-v|^2)\n" + " 3 -- sigmoid: tanh(gamma*u'*v + coef0)\n" + " 4 -- precomputed kernel (kernel values in training_instance_matrix)\n" + "-d degree : set degree in kernel function (default 3)\n" + "-g gamma : set gamma in kernel function (default 1/num_features)\n" + "-r coef0 : set coef0 in kernel function (default 0)\n" + "-c cost : set the parameter C of C-SVC, epsilon-SVR, and nu-SVR (default 1)\n" + "-n nu : set the parameter nu of nu-SVC, one-class SVM, and nu-SVR (default 0.5)\n" + "-p epsilon : set the epsilon in loss function of epsilon-SVR (default 0.1)\n" + "-m cachesize : set cache memory size in MB (default 100)\n" + "-e epsilon : set tolerance of termination criterion (default 0.001)\n" + "-h shrinking : whether to use the shrinking heuristics, 0 or 1 (default 1)\n" + "-b probability_estimates : whether to train a SVC or SVR model for probability estimates, 0 or 1 (default 0)\n" + "-wi weight : set the parameter C of class i to weight*C, for C-SVC (default 1)\n" + "-v n : n-fold cross validation mode\n" + "-q : quiet mode (no outputs)\n" + ); +} + +// svm arguments +struct svm_parameter param; // set by parse_command_line +struct svm_problem prob; // set by read_problem +struct svm_model *model; +struct svm_node *x_space; +int cross_validation; +int nr_fold; + + +double do_cross_validation() +{ + int i; + int total_correct = 0; + double total_error = 0; + double sumv = 0, sumy = 0, sumvv = 0, sumyy = 0, sumvy = 0; + double *target = Malloc(double,prob.l); + double retval = 0.0; + + svm_cross_validation(&prob,¶m,nr_fold,target); + if(param.svm_type == EPSILON_SVR || + param.svm_type == NU_SVR) + { + for(i=0;i<prob.l;i++) + { + double y = prob.y[i]; + double v = target[i]; + total_error += (v-y)*(v-y); + sumv += v; + sumy += y; + sumvv += v*v; + sumyy += y*y; + sumvy += v*y; + } + mexPrintf("Cross Validation Mean squared error = %g\n",total_error/prob.l); + mexPrintf("Cross Validation Squared correlation coefficient = %g\n", + ((prob.l*sumvy-sumv*sumy)*(prob.l*sumvy-sumv*sumy))/ + ((prob.l*sumvv-sumv*sumv)*(prob.l*sumyy-sumy*sumy)) + ); + retval = total_error/prob.l; + } + else + { + for(i=0;i<prob.l;i++) + if(target[i] == prob.y[i]) + ++total_correct; + mexPrintf("Cross Validation Accuracy = %g%%\n",100.0*total_correct/prob.l); + retval = 100.0*total_correct/prob.l; + } + free(target); + return retval; +} + +// nrhs should be 3 +int parse_command_line(int nrhs, const mxArray *prhs[], char *model_file_name) +{ + int i, argc = 1; + char cmd[CMD_LEN]; + char *argv[CMD_LEN/2]; + void (*print_func)(const char *) = print_string_matlab; // default printing to matlab display + + // default values + param.svm_type = C_SVC; + param.kernel_type = RBF; + param.degree = 3; + param.gamma = 0; // 1/num_features + param.coef0 = 0; + param.nu = 0.5; + param.cache_size = 100; + param.C = 1; + param.eps = 1e-3; + param.p = 0.1; + param.shrinking = 1; + param.probability = 0; + param.nr_weight = 0; + param.weight_label = NULL; + param.weight = NULL; + cross_validation = 0; + + if(nrhs <= 1) + return 1; + + if(nrhs > 2) + { + // put options in argv[] + mxGetString(prhs[2], cmd, mxGetN(prhs[2]) + 1); + if((argv[argc] = strtok(cmd, " ")) != NULL) + while((argv[++argc] = strtok(NULL, " ")) != NULL) + ; + } + + // parse options + for(i=1;i<argc;i++) + { + if(argv[i][0] != '-') break; + ++i; + if(i>=argc && argv[i-1][1] != 'q') // since option -q has no parameter + return 1; + switch(argv[i-1][1]) + { + case 's': + param.svm_type = atoi(argv[i]); + break; + case 't': + param.kernel_type = atoi(argv[i]); + break; + case 'd': + param.degree = atoi(argv[i]); + break; + case 'g': + param.gamma = atof(argv[i]); + break; + case 'r': + param.coef0 = atof(argv[i]); + break; + case 'n': + param.nu = atof(argv[i]); + break; + case 'm': + param.cache_size = atof(argv[i]); + break; + case 'c': + param.C = atof(argv[i]); + break; + case 'e': + param.eps = atof(argv[i]); + break; + case 'p': + param.p = atof(argv[i]); + break; + case 'h': + param.shrinking = atoi(argv[i]); + break; + case 'b': + param.probability = atoi(argv[i]); + break; + case 'q': + print_func = &print_null; + i--; + break; + case 'v': + cross_validation = 1; + nr_fold = atoi(argv[i]); + if(nr_fold < 2) + { + mexPrintf("n-fold cross validation: n must >= 2\n"); + return 1; + } + break; + case 'w': + ++param.nr_weight; + param.weight_label = (int *)realloc(param.weight_label,sizeof(int)*param.nr_weight); + param.weight = (double *)realloc(param.weight,sizeof(double)*param.nr_weight); + param.weight_label[param.nr_weight-1] = atoi(&argv[i-1][2]); + param.weight[param.nr_weight-1] = atof(argv[i]); + break; + default: + mexPrintf("Unknown option -%c\n", argv[i-1][1]); + return 1; + } + } + + svm_set_print_string_function(print_func); + + return 0; +} + +// read in a problem (in svmlight format) +int read_problem_dense(const mxArray *label_vec, const mxArray *instance_mat) +{ + int i, j, k; + int elements, max_index, sc, label_vector_row_num; + double *samples, *labels; + + prob.x = NULL; + prob.y = NULL; + x_space = NULL; + + labels = mxGetPr(label_vec); + samples = mxGetPr(instance_mat); + sc = (int)mxGetN(instance_mat); + + elements = 0; + // the number of instance + prob.l = (int)mxGetM(instance_mat); + label_vector_row_num = (int)mxGetM(label_vec); + + if(label_vector_row_num!=prob.l) + { + mexPrintf("Length of label vector does not match # of instances.\n"); + return -1; + } + + if(param.kernel_type == PRECOMPUTED) + elements = prob.l * (sc + 1); + else + { + for(i = 0; i < prob.l; i++) + { + for(k = 0; k < sc; k++) + if(samples[k * prob.l + i] != 0) + elements++; + // count the '-1' element + elements++; + } + } + + prob.y = Malloc(double,prob.l); + prob.x = Malloc(struct svm_node *,prob.l); + x_space = Malloc(struct svm_node, elements); + + max_index = sc; + j = 0; + for(i = 0; i < prob.l; i++) + { + prob.x[i] = &x_space[j]; + prob.y[i] = labels[i]; + + for(k = 0; k < sc; k++) + { + if(param.kernel_type == PRECOMPUTED || samples[k * prob.l + i] != 0) + { + x_space[j].index = k + 1; + x_space[j].value = samples[k * prob.l + i]; + j++; + } + } + x_space[j++].index = -1; + } + + if(param.gamma == 0 && max_index > 0) + param.gamma = 1.0/max_index; + + if(param.kernel_type == PRECOMPUTED) + for(i=0;i<prob.l;i++) + { + if((int)prob.x[i][0].value <= 0 || (int)prob.x[i][0].value > max_index) + { + mexPrintf("Wrong input format: sample_serial_number out of range\n"); + return -1; + } + } + + return 0; +} + +int read_problem_sparse(const mxArray *label_vec, const mxArray *instance_mat) +{ + int i, j, k, low, high; + mwIndex *ir, *jc; + int elements, max_index, num_samples, label_vector_row_num; + double *samples, *labels; + mxArray *instance_mat_col; // transposed instance sparse matrix + + prob.x = NULL; + prob.y = NULL; + x_space = NULL; + + // transpose instance matrix + { + mxArray *prhs[1], *plhs[1]; + prhs[0] = mxDuplicateArray(instance_mat); + if(mexCallMATLAB(1, plhs, 1, prhs, "transpose")) + { + mexPrintf("Error: cannot transpose training instance matrix\n"); + return -1; + } + instance_mat_col = plhs[0]; + mxDestroyArray(prhs[0]); + } + + // each column is one instance + labels = mxGetPr(label_vec); + samples = mxGetPr(instance_mat_col); + ir = mxGetIr(instance_mat_col); + jc = mxGetJc(instance_mat_col); + + num_samples = (int)mxGetNzmax(instance_mat_col); + + // the number of instance + prob.l = (int)mxGetN(instance_mat_col); + label_vector_row_num = (int)mxGetM(label_vec); + + if(label_vector_row_num!=prob.l) + { + mexPrintf("Length of label vector does not match # of instances.\n"); + return -1; + } + + elements = num_samples + prob.l; + max_index = (int)mxGetM(instance_mat_col); + + prob.y = Malloc(double,prob.l); + prob.x = Malloc(struct svm_node *,prob.l); + x_space = Malloc(struct svm_node, elements); + + j = 0; + for(i=0;i<prob.l;i++) + { + prob.x[i] = &x_space[j]; + prob.y[i] = labels[i]; + low = (int)jc[i], high = (int)jc[i+1]; + for(k=low;k<high;k++) + { + x_space[j].index = (int)ir[k] + 1; + x_space[j].value = samples[k]; + j++; + } + x_space[j++].index = -1; + } + + if(param.gamma == 0 && max_index > 0) + param.gamma = 1.0/max_index; + + return 0; +} + +static void fake_answer(int nlhs, mxArray *plhs[]) +{ + int i; + for(i=0;i<nlhs;i++) + plhs[i] = mxCreateDoubleMatrix(0, 0, mxREAL); +} + +// Interface function of matlab +// now assume prhs[0]: label prhs[1]: features +void mexFunction( int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] ) +{ + const char *error_msg; + + // fix random seed to have same results for each run + // (for cross validation and probability estimation) + srand(1); + + if(nlhs > 1) + { + exit_with_help(); + fake_answer(nlhs, plhs); + return; + } + + // Transform the input Matrix to libsvm format + if(nrhs > 1 && nrhs < 4) + { + int err; + + if(!mxIsDouble(prhs[0]) || !mxIsDouble(prhs[1])) { + mexPrintf("Error: label vector and instance matrix must be double\n"); + fake_answer(nlhs, plhs); + return; + } + + if(parse_command_line(nrhs, prhs, NULL)) + { + exit_with_help(); + svm_destroy_param(¶m); + fake_answer(nlhs, plhs); + return; + } + + if(mxIsSparse(prhs[1])) + { + if(param.kernel_type == PRECOMPUTED) + { + // precomputed kernel requires dense matrix, so we make one + mxArray *rhs[1], *lhs[1]; + + rhs[0] = mxDuplicateArray(prhs[1]); + if(mexCallMATLAB(1, lhs, 1, rhs, "full")) + { + mexPrintf("Error: cannot generate a full training instance matrix\n"); + svm_destroy_param(¶m); + fake_answer(nlhs, plhs); + return; + } + err = read_problem_dense(prhs[0], lhs[0]); + mxDestroyArray(lhs[0]); + mxDestroyArray(rhs[0]); + } + else + err = read_problem_sparse(prhs[0], prhs[1]); + } + else + err = read_problem_dense(prhs[0], prhs[1]); + + // svmtrain's original code + error_msg = svm_check_parameter(&prob, ¶m); + + if(err || error_msg) + { + if (error_msg != NULL) + mexPrintf("Error: %s\n", error_msg); + svm_destroy_param(¶m); + free(prob.y); + free(prob.x); + free(x_space); + fake_answer(nlhs, plhs); + return; + } + + if(cross_validation) + { + double *ptr; + plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL); + ptr = mxGetPr(plhs[0]); + ptr[0] = do_cross_validation(); + } + else + { + int nr_feat = (int)mxGetN(prhs[1]); + const char *error_msg; + model = svm_train(&prob, ¶m); + error_msg = model_to_matlab_structure(plhs, nr_feat, model); + if(error_msg) + mexPrintf("Error: can't convert libsvm model to matrix structure: %s\n", error_msg); + svm_free_and_destroy_model(&model); + } + svm_destroy_param(¶m); + free(prob.y); + free(prob.x); + free(x_space); + } + else + { + exit_with_help(); + fake_answer(nlhs, plhs); + return; + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libsvm_osx64/Makefile Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,53 @@ +# This Makefile is used under Linux + +MATLABDIR ?= /usr/local/matlab +# for Mac +# MATLABDIR ?= /opt/local/matlab + +CXX ?= g++ +#CXX = g++-4.1 +CFLAGS = -Wall -Wconversion -O3 -fPIC -I$(MATLABDIR)/extern/include -I.. + +MEX = $(MATLABDIR)/bin/mex +MEX_OPTION = CC\#$(CXX) CXX\#$(CXX) CFLAGS\#"$(CFLAGS)" CXXFLAGS\#"$(CFLAGS)" +# comment the following line if you use MATLAB on 32-bit computer +MEX_OPTION += -largeArrayDims +MEX_EXT = $(shell $(MATLABDIR)/bin/mexext) + +OCTAVEDIR ?= /usr/include/octave +OCTAVE_MEX = env CC=$(CXX) mkoctfile +OCTAVE_MEX_OPTION = --mex +OCTAVE_MEX_EXT = mex +OCTAVE_CFLAGS = -Wall -O3 -fPIC -I$(OCTAVEDIR) -I.. + +all: matlab + +matlab: binary + +octave: + @make MEX="$(OCTAVE_MEX)" MEX_OPTION="$(OCTAVE_MEX_OPTION)" \ + MEX_EXT="$(OCTAVE_MEX_EXT)" CFLAGS="$(OCTAVE_CFLAGS)" \ + binary + +binary: svmpredict.$(MEX_EXT) svmtrain.$(MEX_EXT) libsvmread.$(MEX_EXT) libsvmwrite.$(MEX_EXT) + +svmpredict.$(MEX_EXT): svmpredict.c ../svm.h ../svm.o svm_model_matlab.o + $(MEX) $(MEX_OPTION) svmpredict.c ../svm.o svm_model_matlab.o + +svmtrain.$(MEX_EXT): svmtrain.c ../svm.h ../svm.o svm_model_matlab.o + $(MEX) $(MEX_OPTION) svmtrain.c ../svm.o svm_model_matlab.o + +libsvmread.$(MEX_EXT): libsvmread.c + $(MEX) $(MEX_OPTION) libsvmread.c + +libsvmwrite.$(MEX_EXT): libsvmwrite.c + $(MEX) $(MEX_OPTION) libsvmwrite.c + +svm_model_matlab.o: svm_model_matlab.c ../svm.h + $(CXX) $(CFLAGS) -c svm_model_matlab.c + +../svm.o: ../svm.cpp ../svm.h + make -C .. svm.o + +clean: + rm -f *~ *.o *.mex* *.obj ../svm.o
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libsvm_osx64/README Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,243 @@ +----------------------------------------- +--- MATLAB/OCTAVE interface of LIBSVM --- +----------------------------------------- + +Table of Contents +================= + +- Introduction +- Installation +- Usage +- Returned Model Structure +- Other Utilities +- Examples +- Additional Information + + +Introduction +============ + +This tool provides a simple interface to LIBSVM, a library for support vector +machines (http://www.csie.ntu.edu.tw/~cjlin/libsvm). It is very easy to use as +the usage and the way of specifying parameters are the same as that of LIBSVM. + +Installation +============ + +On Windows systems, pre-built binary files are already in the +directory '..\windows', so no need to conduct installation. Now we +provide binary files only for 64bit MATLAB on Windows. If you would +like to re-build the package, please rely on the following steps. + +We recommend using make.m on both MATLAB and OCTAVE. Just type 'make' +to build 'libsvmread.mex', 'libsvmwrite.mex', 'svmtrain.mex', and +'svmpredict.mex'. + +On MATLAB or Octave: + + >> make + +If make.m does not work on MATLAB (especially for Windows), try 'mex +-setup' to choose a suitable compiler for mex. Make sure your compiler +is accessible and workable. Then type 'make' to start the +installation. + +Example: + + matlab>> mex -setup + (ps: MATLAB will show the following messages to setup default compiler.) + Please choose your compiler for building external interface (MEX) files: + Would you like mex to locate installed compilers [y]/n? y + Select a compiler: + [1] Microsoft Visual C/C++ version 7.1 in C:\Program Files\Microsoft Visual Studio + [0] None + Compiler: 1 + Please verify your choices: + Compiler: Microsoft Visual C/C++ 7.1 + Location: C:\Program Files\Microsoft Visual Studio + Are these correct?([y]/n): y + + matlab>> make + +On Unix systems, if neither make.m nor 'mex -setup' works, please use +Makefile and type 'make' in a command window. Note that we assume +your MATLAB is installed in '/usr/local/matlab'. If not, please change +MATLABDIR in Makefile. + +Example: + linux> make + +To use octave, type 'make octave': + +Example: + linux> make octave + +For a list of supported/compatible compilers for MATLAB, please check +the following page: + +http://www.mathworks.com/support/compilers/current_release/ + +Usage +===== + +matlab> model = svmtrain(training_label_vector, training_instance_matrix [, 'libsvm_options']); + + -training_label_vector: + An m by 1 vector of training labels (type must be double). + -training_instance_matrix: + An m by n matrix of m training instances with n features. + It can be dense or sparse (type must be double). + -libsvm_options: + A string of training options in the same format as that of LIBSVM. + +matlab> [predicted_label, accuracy, decision_values/prob_estimates] = svmpredict(testing_label_vector, testing_instance_matrix, model [, 'libsvm_options']); + + -testing_label_vector: + An m by 1 vector of prediction labels. If labels of test + data are unknown, simply use any random values. (type must be double) + -testing_instance_matrix: + An m by n matrix of m testing instances with n features. + It can be dense or sparse. (type must be double) + -model: + The output of svmtrain. + -libsvm_options: + A string of testing options in the same format as that of LIBSVM. + +Returned Model Structure +======================== + +The 'svmtrain' function returns a model which can be used for future +prediction. It is a structure and is organized as [Parameters, nr_class, +totalSV, rho, Label, ProbA, ProbB, nSV, sv_coef, SVs]: + + -Parameters: parameters + -nr_class: number of classes; = 2 for regression/one-class svm + -totalSV: total #SV + -rho: -b of the decision function(s) wx+b + -Label: label of each class; empty for regression/one-class SVM + -ProbA: pairwise probability information; empty if -b 0 or in one-class SVM + -ProbB: pairwise probability information; empty if -b 0 or in one-class SVM + -nSV: number of SVs for each class; empty for regression/one-class SVM + -sv_coef: coefficients for SVs in decision functions + -SVs: support vectors + +If you do not use the option '-b 1', ProbA and ProbB are empty +matrices. If the '-v' option is specified, cross validation is +conducted and the returned model is just a scalar: cross-validation +accuracy for classification and mean-squared error for regression. + +More details about this model can be found in LIBSVM FAQ +(http://www.csie.ntu.edu.tw/~cjlin/libsvm/faq.html) and LIBSVM +implementation document +(http://www.csie.ntu.edu.tw/~cjlin/papers/libsvm.pdf). + +Result of Prediction +==================== + +The function 'svmpredict' has three outputs. The first one, +predictd_label, is a vector of predicted labels. The second output, +accuracy, is a vector including accuracy (for classification), mean +squared error, and squared correlation coefficient (for regression). +The third is a matrix containing decision values or probability +estimates (if '-b 1' is specified). If k is the number of classes +in training data, for decision values, each row includes results of +predicting k(k-1)/2 binary-class SVMs. For classification, k = 1 is a +special case. Decision value +1 is returned for each testing instance, +instead of an empty vector. For probabilities, each row contains k values +indicating the probability that the testing instance is in each class. +Note that the order of classes here is the same as 'Label' field +in the model structure. + +Other Utilities +=============== + +A matlab function libsvmread reads files in LIBSVM format: + +[label_vector, instance_matrix] = libsvmread('data.txt'); + +Two outputs are labels and instances, which can then be used as inputs +of svmtrain or svmpredict. + +A matlab function libsvmwrite writes Matlab matrix to a file in LIBSVM format: + +libsvmwrite('data.txt', label_vector, instance_matrix] + +The instance_matrix must be a sparse matrix. (type must be double) +For 32bit and 64bit MATLAB on Windows, pre-built binary files are ready +in the directory `..\windows', but in future releases, we will only +include 64bit MATLAB binary files. + +These codes are prepared by Rong-En Fan and Kai-Wei Chang from National +Taiwan University. + +Examples +======== + +Train and test on the provided data heart_scale: + +matlab> [heart_scale_label, heart_scale_inst] = libsvmread('../heart_scale'); +matlab> model = svmtrain(heart_scale_label, heart_scale_inst, '-c 1 -g 0.07'); +matlab> [predict_label, accuracy, dec_values] = svmpredict(heart_scale_label, heart_scale_inst, model); % test the training data + +For probability estimates, you need '-b 1' for training and testing: + +matlab> [heart_scale_label, heart_scale_inst] = libsvmread('../heart_scale'); +matlab> model = svmtrain(heart_scale_label, heart_scale_inst, '-c 1 -g 0.07 -b 1'); +matlab> [heart_scale_label, heart_scale_inst] = libsvmread('../heart_scale'); +matlab> [predict_label, accuracy, prob_estimates] = svmpredict(heart_scale_label, heart_scale_inst, model, '-b 1'); + +To use precomputed kernel, you must include sample serial number as +the first column of the training and testing data (assume your kernel +matrix is K, # of instances is n): + +matlab> K1 = [(1:n)', K]; % include sample serial number as first column +matlab> model = svmtrain(label_vector, K1, '-t 4'); +matlab> [predict_label, accuracy, dec_values] = svmpredict(label_vector, K1, model); % test the training data + +We give the following detailed example by splitting heart_scale into +150 training and 120 testing data. Constructing a linear kernel +matrix and then using the precomputed kernel gives exactly the same +testing error as using the LIBSVM built-in linear kernel. + +matlab> [heart_scale_label, heart_scale_inst] = libsvmread('../heart_scale'); +matlab> +matlab> % Split Data +matlab> train_data = heart_scale_inst(1:150,:); +matlab> train_label = heart_scale_label(1:150,:); +matlab> test_data = heart_scale_inst(151:270,:); +matlab> test_label = heart_scale_label(151:270,:); +matlab> +matlab> % Linear Kernel +matlab> model_linear = svmtrain(train_label, train_data, '-t 0'); +matlab> [predict_label_L, accuracy_L, dec_values_L] = svmpredict(test_label, test_data, model_linear); +matlab> +matlab> % Precomputed Kernel +matlab> model_precomputed = svmtrain(train_label, [(1:150)', train_data*train_data'], '-t 4'); +matlab> [predict_label_P, accuracy_P, dec_values_P] = svmpredict(test_label, [(1:120)', test_data*train_data'], model_precomputed); +matlab> +matlab> accuracy_L % Display the accuracy using linear kernel +matlab> accuracy_P % Display the accuracy using precomputed kernel + +Note that for testing, you can put anything in the +testing_label_vector. For more details of precomputed kernels, please +read the section ``Precomputed Kernels'' in the README of the LIBSVM +package. + +Additional Information +====================== + +This interface was initially written by Jun-Cheng Chen, Kuan-Jen Peng, +Chih-Yuan Yang and Chih-Huai Cheng from Department of Computer +Science, National Taiwan University. The current version was prepared +by Rong-En Fan and Ting-Fan Wu. If you find this tool useful, please +cite LIBSVM as follows + +Chih-Chung Chang and Chih-Jen Lin, LIBSVM : a library for support +vector machines. ACM Transactions on Intelligent Systems and +Technology, 2:27:1--27:27, 2011. Software available at +http://www.csie.ntu.edu.tw/~cjlin/libsvm + +For any question, please contact Chih-Jen Lin <cjlin@csie.ntu.edu.tw>, +or check the FAQ page: + +http://www.csie.ntu.edu.tw/~cjlin/libsvm/faq.html#/Q9:_MATLAB_interface
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libsvm_osx64/libsvmread.c Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,212 @@ +#include <stdio.h> +#include <string.h> +#include <stdlib.h> +#include <ctype.h> +#include <errno.h> + +#include "mex.h" + +#ifdef MX_API_VER +#if MX_API_VER < 0x07030000 +typedef int mwIndex; +#endif +#endif +#ifndef max +#define max(x,y) (((x)>(y))?(x):(y)) +#endif +#ifndef min +#define min(x,y) (((x)<(y))?(x):(y)) +#endif + +void exit_with_help() +{ + mexPrintf( + "Usage: [label_vector, instance_matrix] = libsvmread('filename');\n" + ); +} + +static void fake_answer(mxArray *plhs[]) +{ + plhs[0] = mxCreateDoubleMatrix(0, 0, mxREAL); + plhs[1] = mxCreateDoubleMatrix(0, 0, mxREAL); +} + +static char *line; +static int max_line_len; + +static char* readline(FILE *input) +{ + int len; + + if(fgets(line,max_line_len,input) == NULL) + return NULL; + + while(strrchr(line,'\n') == NULL) + { + max_line_len *= 2; + line = (char *) realloc(line, max_line_len); + len = (int) strlen(line); + if(fgets(line+len,max_line_len-len,input) == NULL) + break; + } + return line; +} + +// read in a problem (in libsvm format) +void read_problem(const char *filename, mxArray *plhs[]) +{ + int max_index, min_index, inst_max_index, i; + long elements, k; + FILE *fp = fopen(filename,"r"); + int l = 0; + char *endptr; + mwIndex *ir, *jc; + double *labels, *samples; + + if(fp == NULL) + { + mexPrintf("can't open input file %s\n",filename); + fake_answer(plhs); + return; + } + + max_line_len = 1024; + line = (char *) malloc(max_line_len*sizeof(char)); + + max_index = 0; + min_index = 1; // our index starts from 1 + elements = 0; + while(readline(fp) != NULL) + { + char *idx, *val; + // features + int index = 0; + + inst_max_index = -1; // strtol gives 0 if wrong format, and precomputed kernel has <index> start from 0 + strtok(line," \t"); // label + while (1) + { + idx = strtok(NULL,":"); // index:value + val = strtok(NULL," \t"); + if(val == NULL) + break; + + errno = 0; + index = (int) strtol(idx,&endptr,10); + if(endptr == idx || errno != 0 || *endptr != '\0' || index <= inst_max_index) + { + mexPrintf("Wrong input format at line %d\n",l+1); + fake_answer(plhs); + return; + } + else + inst_max_index = index; + + min_index = min(min_index, index); + elements++; + } + max_index = max(max_index, inst_max_index); + l++; + } + rewind(fp); + + // y + plhs[0] = mxCreateDoubleMatrix(l, 1, mxREAL); + // x^T + if (min_index <= 0) + plhs[1] = mxCreateSparse(max_index-min_index+1, l, elements, mxREAL); + else + plhs[1] = mxCreateSparse(max_index, l, elements, mxREAL); + + labels = mxGetPr(plhs[0]); + samples = mxGetPr(plhs[1]); + ir = mxGetIr(plhs[1]); + jc = mxGetJc(plhs[1]); + + k=0; + for(i=0;i<l;i++) + { + char *idx, *val, *label; + jc[i] = k; + + readline(fp); + + label = strtok(line," \t\n"); + if(label == NULL) + { + mexPrintf("Empty line at line %d\n",i+1); + fake_answer(plhs); + return; + } + labels[i] = strtod(label,&endptr); + if(endptr == label || *endptr != '\0') + { + mexPrintf("Wrong input format at line %d\n",i+1); + fake_answer(plhs); + return; + } + + // features + while(1) + { + idx = strtok(NULL,":"); + val = strtok(NULL," \t"); + if(val == NULL) + break; + + ir[k] = (mwIndex) (strtol(idx,&endptr,10) - min_index); // precomputed kernel has <index> start from 0 + + errno = 0; + samples[k] = strtod(val,&endptr); + if (endptr == val || errno != 0 || (*endptr != '\0' && !isspace(*endptr))) + { + mexPrintf("Wrong input format at line %d\n",i+1); + fake_answer(plhs); + return; + } + ++k; + } + } + jc[l] = k; + + fclose(fp); + free(line); + + { + mxArray *rhs[1], *lhs[1]; + rhs[0] = plhs[1]; + if(mexCallMATLAB(1, lhs, 1, rhs, "transpose")) + { + mexPrintf("Error: cannot transpose problem\n"); + fake_answer(plhs); + return; + } + plhs[1] = lhs[0]; + } +} + +void mexFunction( int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] ) +{ + if(nrhs == 1) + { + char filename[256]; + + mxGetString(prhs[0], filename, mxGetN(prhs[0]) + 1); + + if(filename == NULL) + { + mexPrintf("Error: filename is NULL\n"); + return; + } + + read_problem(filename, plhs); + } + else + { + exit_with_help(); + fake_answer(plhs); + return; + } +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libsvm_osx64/libsvmwrite.c Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,106 @@ +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include "mex.h" + +#ifdef MX_API_VER +#if MX_API_VER < 0x07030000 +typedef int mwIndex; +#endif +#endif + +void exit_with_help() +{ + mexPrintf( + "Usage: libsvmwrite('filename', label_vector, instance_matrix);\n" + ); +} + +void libsvmwrite(const char *filename, const mxArray *label_vec, const mxArray *instance_mat) +{ + FILE *fp = fopen(filename,"w"); + int i, k, low, high, l; + mwIndex *ir, *jc; + int label_vector_row_num; + double *samples, *labels; + mxArray *instance_mat_col; // instance sparse matrix in column format + + if(fp ==NULL) + { + mexPrintf("can't open output file %s\n",filename); + return; + } + + // transpose instance matrix + { + mxArray *prhs[1], *plhs[1]; + prhs[0] = mxDuplicateArray(instance_mat); + if(mexCallMATLAB(1, plhs, 1, prhs, "transpose")) + { + mexPrintf("Error: cannot transpose instance matrix\n"); + return; + } + instance_mat_col = plhs[0]; + mxDestroyArray(prhs[0]); + } + + // the number of instance + l = (int) mxGetN(instance_mat_col); + label_vector_row_num = (int) mxGetM(label_vec); + + if(label_vector_row_num!=l) + { + mexPrintf("Length of label vector does not match # of instances.\n"); + return; + } + + // each column is one instance + labels = mxGetPr(label_vec); + samples = mxGetPr(instance_mat_col); + ir = mxGetIr(instance_mat_col); + jc = mxGetJc(instance_mat_col); + + for(i=0;i<l;i++) + { + fprintf(fp,"%g", labels[i]); + + low = (int) jc[i], high = (int) jc[i+1]; + for(k=low;k<high;k++) + fprintf(fp," %ld:%g", ir[k]+1, samples[k]); + + fprintf(fp,"\n"); + } + + fclose(fp); + return; +} + +void mexFunction( int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] ) +{ + // Transform the input Matrix to libsvm format + if(nrhs == 3) + { + char filename[256]; + if(!mxIsDouble(prhs[1]) || !mxIsDouble(prhs[2])) + { + mexPrintf("Error: label vector and instance matrix must be double\n"); + return; + } + + mxGetString(prhs[0], filename, mxGetN(prhs[0])+1); + + if(mxIsSparse(prhs[2])) + libsvmwrite(filename, prhs[1], prhs[2]); + else + { + mexPrintf("Instance_matrix must be sparse\n"); + return; + } + } + else + { + exit_with_help(); + return; + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libsvm_osx64/make.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,21 @@ +% This make.m is for MATLAB and OCTAVE under Windows, Mac, and Unix + +try + Type = ver; + % This part is for OCTAVE + if(strcmp(Type(1).Name, 'Octave') == 1) + mex libsvmread.c + mex libsvmwrite.c + mex svmtrain.c ../svm.cpp svm_model_matlab.c + mex svmpredict.c ../svm.cpp svm_model_matlab.c + % This part is for MATLAB + % Add -largeArrayDims on 64-bit machines of MATLAB + else + mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims libsvmread.c + mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims libsvmwrite.c + mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims svmtrain.c ../svm.cpp svm_model_matlab.c + mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims svmpredict.c ../svm.cpp svm_model_matlab.c + end +catch + fprintf('If make.m fails, please check README about detailed instructions.\n'); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libsvm_osx64/svm_model_matlab.c Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,375 @@ +#include <stdlib.h> +#include <string.h> +#include "../svm.h" + +#include "mex.h" + +#ifdef MX_API_VER +#if MX_API_VER < 0x07030000 +typedef int mwIndex; +#endif +#endif + +#define NUM_OF_RETURN_FIELD 11 + +#define Malloc(type,n) (type *)malloc((n)*sizeof(type)) + +static const char *field_names[] = { + "Parameters", + "nr_class", + "totalSV", + "rho", + "Label", + "sv_indices", + "ProbA", + "ProbB", + "nSV", + "sv_coef", + "SVs" +}; + +const char *model_to_matlab_structure(mxArray *plhs[], int num_of_feature, struct svm_model *model) +{ + int i, j, n; + double *ptr; + mxArray *return_model, **rhs; + int out_id = 0; + + rhs = (mxArray **)mxMalloc(sizeof(mxArray *)*NUM_OF_RETURN_FIELD); + + // Parameters + rhs[out_id] = mxCreateDoubleMatrix(5, 1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + ptr[0] = model->param.svm_type; + ptr[1] = model->param.kernel_type; + ptr[2] = model->param.degree; + ptr[3] = model->param.gamma; + ptr[4] = model->param.coef0; + out_id++; + + // nr_class + rhs[out_id] = mxCreateDoubleMatrix(1, 1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + ptr[0] = model->nr_class; + out_id++; + + // total SV + rhs[out_id] = mxCreateDoubleMatrix(1, 1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + ptr[0] = model->l; + out_id++; + + // rho + n = model->nr_class*(model->nr_class-1)/2; + rhs[out_id] = mxCreateDoubleMatrix(n, 1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + for(i = 0; i < n; i++) + ptr[i] = model->rho[i]; + out_id++; + + // Label + if(model->label) + { + rhs[out_id] = mxCreateDoubleMatrix(model->nr_class, 1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + for(i = 0; i < model->nr_class; i++) + ptr[i] = model->label[i]; + } + else + rhs[out_id] = mxCreateDoubleMatrix(0, 0, mxREAL); + out_id++; + + // sv_indices + if(model->sv_indices) + { + rhs[out_id] = mxCreateDoubleMatrix(model->l, 1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + for(i = 0; i < model->l; i++) + ptr[i] = model->sv_indices[i]; + } + else + rhs[out_id] = mxCreateDoubleMatrix(0, 0, mxREAL); + out_id++; + + // probA + if(model->probA != NULL) + { + rhs[out_id] = mxCreateDoubleMatrix(n, 1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + for(i = 0; i < n; i++) + ptr[i] = model->probA[i]; + } + else + rhs[out_id] = mxCreateDoubleMatrix(0, 0, mxREAL); + out_id ++; + + // probB + if(model->probB != NULL) + { + rhs[out_id] = mxCreateDoubleMatrix(n, 1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + for(i = 0; i < n; i++) + ptr[i] = model->probB[i]; + } + else + rhs[out_id] = mxCreateDoubleMatrix(0, 0, mxREAL); + out_id++; + + // nSV + if(model->nSV) + { + rhs[out_id] = mxCreateDoubleMatrix(model->nr_class, 1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + for(i = 0; i < model->nr_class; i++) + ptr[i] = model->nSV[i]; + } + else + rhs[out_id] = mxCreateDoubleMatrix(0, 0, mxREAL); + out_id++; + + // sv_coef + rhs[out_id] = mxCreateDoubleMatrix(model->l, model->nr_class-1, mxREAL); + ptr = mxGetPr(rhs[out_id]); + for(i = 0; i < model->nr_class-1; i++) + for(j = 0; j < model->l; j++) + ptr[(i*(model->l))+j] = model->sv_coef[i][j]; + out_id++; + + // SVs + { + int ir_index, nonzero_element; + mwIndex *ir, *jc; + mxArray *pprhs[1], *pplhs[1]; + + if(model->param.kernel_type == PRECOMPUTED) + { + nonzero_element = model->l; + num_of_feature = 1; + } + else + { + nonzero_element = 0; + for(i = 0; i < model->l; i++) { + j = 0; + while(model->SV[i][j].index != -1) + { + nonzero_element++; + j++; + } + } + } + + // SV in column, easier accessing + rhs[out_id] = mxCreateSparse(num_of_feature, model->l, nonzero_element, mxREAL); + ir = mxGetIr(rhs[out_id]); + jc = mxGetJc(rhs[out_id]); + ptr = mxGetPr(rhs[out_id]); + jc[0] = ir_index = 0; + for(i = 0;i < model->l; i++) + { + if(model->param.kernel_type == PRECOMPUTED) + { + // make a (1 x model->l) matrix + ir[ir_index] = 0; + ptr[ir_index] = model->SV[i][0].value; + ir_index++; + jc[i+1] = jc[i] + 1; + } + else + { + int x_index = 0; + while (model->SV[i][x_index].index != -1) + { + ir[ir_index] = model->SV[i][x_index].index - 1; + ptr[ir_index] = model->SV[i][x_index].value; + ir_index++, x_index++; + } + jc[i+1] = jc[i] + x_index; + } + } + // transpose back to SV in row + pprhs[0] = rhs[out_id]; + if(mexCallMATLAB(1, pplhs, 1, pprhs, "transpose")) + return "cannot transpose SV matrix"; + rhs[out_id] = pplhs[0]; + out_id++; + } + + /* Create a struct matrix contains NUM_OF_RETURN_FIELD fields */ + return_model = mxCreateStructMatrix(1, 1, NUM_OF_RETURN_FIELD, field_names); + + /* Fill struct matrix with input arguments */ + for(i = 0; i < NUM_OF_RETURN_FIELD; i++) + mxSetField(return_model,0,field_names[i],mxDuplicateArray(rhs[i])); + /* return */ + plhs[0] = return_model; + mxFree(rhs); + + return NULL; +} + +struct svm_model *matlab_matrix_to_model(const mxArray *matlab_struct, const char **msg) +{ + int i, j, n, num_of_fields; + double *ptr; + int id = 0; + struct svm_node *x_space; + struct svm_model *model; + mxArray **rhs; + + num_of_fields = mxGetNumberOfFields(matlab_struct); + if(num_of_fields != NUM_OF_RETURN_FIELD) + { + *msg = "number of return field is not correct"; + return NULL; + } + rhs = (mxArray **) mxMalloc(sizeof(mxArray *)*num_of_fields); + + for(i=0;i<num_of_fields;i++) + rhs[i] = mxGetFieldByNumber(matlab_struct, 0, i); + + model = Malloc(struct svm_model, 1); + model->rho = NULL; + model->probA = NULL; + model->probB = NULL; + model->label = NULL; + model->sv_indices = NULL; + model->nSV = NULL; + model->free_sv = 1; // XXX + + ptr = mxGetPr(rhs[id]); + model->param.svm_type = (int)ptr[0]; + model->param.kernel_type = (int)ptr[1]; + model->param.degree = (int)ptr[2]; + model->param.gamma = ptr[3]; + model->param.coef0 = ptr[4]; + id++; + + ptr = mxGetPr(rhs[id]); + model->nr_class = (int)ptr[0]; + id++; + + ptr = mxGetPr(rhs[id]); + model->l = (int)ptr[0]; + id++; + + // rho + n = model->nr_class * (model->nr_class-1)/2; + model->rho = (double*) malloc(n*sizeof(double)); + ptr = mxGetPr(rhs[id]); + for(i=0;i<n;i++) + model->rho[i] = ptr[i]; + id++; + + // label + if(mxIsEmpty(rhs[id]) == 0) + { + model->label = (int*) malloc(model->nr_class*sizeof(int)); + ptr = mxGetPr(rhs[id]); + for(i=0;i<model->nr_class;i++) + model->label[i] = (int)ptr[i]; + } + id++; + + // sv_indices + if(mxIsEmpty(rhs[id]) == 0) + { + model->sv_indices = (int*) malloc(model->l*sizeof(int)); + ptr = mxGetPr(rhs[id]); + for(i=0;i<model->l;i++) + model->sv_indices[i] = (int)ptr[i]; + } + id++; + + // probA + if(mxIsEmpty(rhs[id]) == 0) + { + model->probA = (double*) malloc(n*sizeof(double)); + ptr = mxGetPr(rhs[id]); + for(i=0;i<n;i++) + model->probA[i] = ptr[i]; + } + id++; + + // probB + if(mxIsEmpty(rhs[id]) == 0) + { + model->probB = (double*) malloc(n*sizeof(double)); + ptr = mxGetPr(rhs[id]); + for(i=0;i<n;i++) + model->probB[i] = ptr[i]; + } + id++; + + // nSV + if(mxIsEmpty(rhs[id]) == 0) + { + model->nSV = (int*) malloc(model->nr_class*sizeof(int)); + ptr = mxGetPr(rhs[id]); + for(i=0;i<model->nr_class;i++) + model->nSV[i] = (int)ptr[i]; + } + id++; + + // sv_coef + ptr = mxGetPr(rhs[id]); + model->sv_coef = (double**) malloc((model->nr_class-1)*sizeof(double)); + for( i=0 ; i< model->nr_class -1 ; i++ ) + model->sv_coef[i] = (double*) malloc((model->l)*sizeof(double)); + for(i = 0; i < model->nr_class - 1; i++) + for(j = 0; j < model->l; j++) + model->sv_coef[i][j] = ptr[i*(model->l)+j]; + id++; + + // SV + { + int sr, sc, elements; + int num_samples; + mwIndex *ir, *jc; + mxArray *pprhs[1], *pplhs[1]; + + // transpose SV + pprhs[0] = rhs[id]; + if(mexCallMATLAB(1, pplhs, 1, pprhs, "transpose")) + { + svm_free_and_destroy_model(&model); + *msg = "cannot transpose SV matrix"; + return NULL; + } + rhs[id] = pplhs[0]; + + sr = (int)mxGetN(rhs[id]); + sc = (int)mxGetM(rhs[id]); + + ptr = mxGetPr(rhs[id]); + ir = mxGetIr(rhs[id]); + jc = mxGetJc(rhs[id]); + + num_samples = (int)mxGetNzmax(rhs[id]); + + elements = num_samples + sr; + + model->SV = (struct svm_node **) malloc(sr * sizeof(struct svm_node *)); + x_space = (struct svm_node *)malloc(elements * sizeof(struct svm_node)); + + // SV is in column + for(i=0;i<sr;i++) + { + int low = (int)jc[i], high = (int)jc[i+1]; + int x_index = 0; + model->SV[i] = &x_space[low+i]; + for(j=low;j<high;j++) + { + model->SV[i][x_index].index = (int)ir[j] + 1; + model->SV[i][x_index].value = ptr[j]; + x_index++; + } + model->SV[i][x_index].index = -1; + } + + id++; + } + mxFree(rhs); + + return model; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libsvm_osx64/svm_model_matlab.h Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,2 @@ +const char *model_to_matlab_structure(mxArray *plhs[], int num_of_feature, struct svm_model *model); +struct svm_model *matlab_matrix_to_model(const mxArray *matlab_struct, const char **error_message);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libsvm_osx64/svmpredict.c Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,358 @@ +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include "../svm.h" + +#include "mex.h" +#include "svm_model_matlab.h" + +#ifdef MX_API_VER +#if MX_API_VER < 0x07030000 +typedef int mwIndex; +#endif +#endif + +#define CMD_LEN 2048 + +int print_null(const char *s,...) {} +int (*info)(const char *fmt,...) = &mexPrintf; + +void read_sparse_instance(const mxArray *prhs, int index, struct svm_node *x) +{ + int i, j, low, high; + mwIndex *ir, *jc; + double *samples; + + ir = mxGetIr(prhs); + jc = mxGetJc(prhs); + samples = mxGetPr(prhs); + + // each column is one instance + j = 0; + low = (int)jc[index], high = (int)jc[index+1]; + for(i=low;i<high;i++) + { + x[j].index = (int)ir[i] + 1; + x[j].value = samples[i]; + j++; + } + x[j].index = -1; +} + +static void fake_answer(mxArray *plhs[]) +{ + plhs[0] = mxCreateDoubleMatrix(0, 0, mxREAL); + plhs[1] = mxCreateDoubleMatrix(0, 0, mxREAL); + plhs[2] = mxCreateDoubleMatrix(0, 0, mxREAL); +} + +void predict(mxArray *plhs[], const mxArray *prhs[], struct svm_model *model, const int predict_probability) +{ + int label_vector_row_num, label_vector_col_num; + int feature_number, testing_instance_number; + int instance_index; + double *ptr_instance, *ptr_label, *ptr_predict_label; + double *ptr_prob_estimates, *ptr_dec_values, *ptr; + struct svm_node *x; + mxArray *pplhs[1]; // transposed instance sparse matrix + + int correct = 0; + int total = 0; + double error = 0; + double sump = 0, sumt = 0, sumpp = 0, sumtt = 0, sumpt = 0; + + int svm_type=svm_get_svm_type(model); + int nr_class=svm_get_nr_class(model); + double *prob_estimates=NULL; + + // prhs[1] = testing instance matrix + feature_number = (int)mxGetN(prhs[1]); + testing_instance_number = (int)mxGetM(prhs[1]); + label_vector_row_num = (int)mxGetM(prhs[0]); + label_vector_col_num = (int)mxGetN(prhs[0]); + + if(label_vector_row_num!=testing_instance_number) + { + mexPrintf("Length of label vector does not match # of instances.\n"); + fake_answer(plhs); + return; + } + if(label_vector_col_num!=1) + { + mexPrintf("label (1st argument) should be a vector (# of column is 1).\n"); + fake_answer(plhs); + return; + } + + ptr_instance = mxGetPr(prhs[1]); + ptr_label = mxGetPr(prhs[0]); + + // transpose instance matrix + if(mxIsSparse(prhs[1])) + { + if(model->param.kernel_type == PRECOMPUTED) + { + // precomputed kernel requires dense matrix, so we make one + mxArray *rhs[1], *lhs[1]; + rhs[0] = mxDuplicateArray(prhs[1]); + if(mexCallMATLAB(1, lhs, 1, rhs, "full")) + { + mexPrintf("Error: cannot full testing instance matrix\n"); + fake_answer(plhs); + return; + } + ptr_instance = mxGetPr(lhs[0]); + mxDestroyArray(rhs[0]); + } + else + { + mxArray *pprhs[1]; + pprhs[0] = mxDuplicateArray(prhs[1]); + if(mexCallMATLAB(1, pplhs, 1, pprhs, "transpose")) + { + mexPrintf("Error: cannot transpose testing instance matrix\n"); + fake_answer(plhs); + return; + } + } + } + + if(predict_probability) + { + if(svm_type==NU_SVR || svm_type==EPSILON_SVR) + info("Prob. model for test data: target value = predicted value + z,\nz: Laplace distribution e^(-|z|/sigma)/(2sigma),sigma=%g\n",svm_get_svr_probability(model)); + else + prob_estimates = (double *) malloc(nr_class*sizeof(double)); + } + + plhs[0] = mxCreateDoubleMatrix(testing_instance_number, 1, mxREAL); + if(predict_probability) + { + // prob estimates are in plhs[2] + if(svm_type==C_SVC || svm_type==NU_SVC) + plhs[2] = mxCreateDoubleMatrix(testing_instance_number, nr_class, mxREAL); + else + plhs[2] = mxCreateDoubleMatrix(0, 0, mxREAL); + } + else + { + // decision values are in plhs[2] + if(svm_type == ONE_CLASS || + svm_type == EPSILON_SVR || + svm_type == NU_SVR || + nr_class == 1) // if only one class in training data, decision values are still returned. + plhs[2] = mxCreateDoubleMatrix(testing_instance_number, 1, mxREAL); + else + plhs[2] = mxCreateDoubleMatrix(testing_instance_number, nr_class*(nr_class-1)/2, mxREAL); + } + + ptr_predict_label = mxGetPr(plhs[0]); + ptr_prob_estimates = mxGetPr(plhs[2]); + ptr_dec_values = mxGetPr(plhs[2]); + x = (struct svm_node*)malloc((feature_number+1)*sizeof(struct svm_node) ); + for(instance_index=0;instance_index<testing_instance_number;instance_index++) + { + int i; + double target_label, predict_label; + + target_label = ptr_label[instance_index]; + + if(mxIsSparse(prhs[1]) && model->param.kernel_type != PRECOMPUTED) // prhs[1]^T is still sparse + read_sparse_instance(pplhs[0], instance_index, x); + else + { + for(i=0;i<feature_number;i++) + { + x[i].index = i+1; + x[i].value = ptr_instance[testing_instance_number*i+instance_index]; + } + x[feature_number].index = -1; + } + + if(predict_probability) + { + if(svm_type==C_SVC || svm_type==NU_SVC) + { + predict_label = svm_predict_probability(model, x, prob_estimates); + ptr_predict_label[instance_index] = predict_label; + for(i=0;i<nr_class;i++) + ptr_prob_estimates[instance_index + i * testing_instance_number] = prob_estimates[i]; + } else { + predict_label = svm_predict(model,x); + ptr_predict_label[instance_index] = predict_label; + } + } + else + { + if(svm_type == ONE_CLASS || + svm_type == EPSILON_SVR || + svm_type == NU_SVR) + { + double res; + predict_label = svm_predict_values(model, x, &res); + ptr_dec_values[instance_index] = res; + } + else + { + double *dec_values = (double *) malloc(sizeof(double) * nr_class*(nr_class-1)/2); + predict_label = svm_predict_values(model, x, dec_values); + if(nr_class == 1) + ptr_dec_values[instance_index] = 1; + else + for(i=0;i<(nr_class*(nr_class-1))/2;i++) + ptr_dec_values[instance_index + i * testing_instance_number] = dec_values[i]; + free(dec_values); + } + ptr_predict_label[instance_index] = predict_label; + } + + if(predict_label == target_label) + ++correct; + error += (predict_label-target_label)*(predict_label-target_label); + sump += predict_label; + sumt += target_label; + sumpp += predict_label*predict_label; + sumtt += target_label*target_label; + sumpt += predict_label*target_label; + ++total; + } + if(svm_type==NU_SVR || svm_type==EPSILON_SVR) + { + info("Mean squared error = %g (regression)\n",error/total); + info("Squared correlation coefficient = %g (regression)\n", + ((total*sumpt-sump*sumt)*(total*sumpt-sump*sumt))/ + ((total*sumpp-sump*sump)*(total*sumtt-sumt*sumt)) + ); + } + else + info("Accuracy = %g%% (%d/%d) (classification)\n", + (double)correct/total*100,correct,total); + + // return accuracy, mean squared error, squared correlation coefficient + plhs[1] = mxCreateDoubleMatrix(3, 1, mxREAL); + ptr = mxGetPr(plhs[1]); + ptr[0] = (double)correct/total*100; + ptr[1] = error/total; + ptr[2] = ((total*sumpt-sump*sumt)*(total*sumpt-sump*sumt))/ + ((total*sumpp-sump*sump)*(total*sumtt-sumt*sumt)); + + free(x); + if(prob_estimates != NULL) + free(prob_estimates); +} + +void exit_with_help() +{ + mexPrintf( + "Usage: [predicted_label, accuracy, decision_values/prob_estimates] = svmpredict(testing_label_vector, testing_instance_matrix, model, 'libsvm_options')\n" + "Parameters:\n" + " model: SVM model structure from svmtrain.\n" + " libsvm_options:\n" + " -b probability_estimates: whether to predict probability estimates, 0 or 1 (default 0); one-class SVM not supported yet\n" + " -q : quiet mode (no outputs)\n" + "Returns:\n" + " predicted_label: SVM prediction output vector.\n" + " accuracy: a vector with accuracy, mean squared error, squared correlation coefficient.\n" + " prob_estimates: If selected, probability estimate vector.\n" + ); +} + +void mexFunction( int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] ) +{ + int prob_estimate_flag = 0; + struct svm_model *model; + info = &mexPrintf; + + if(nrhs > 4 || nrhs < 3) + { + exit_with_help(); + fake_answer(plhs); + return; + } + + if(!mxIsDouble(prhs[0]) || !mxIsDouble(prhs[1])) { + mexPrintf("Error: label vector and instance matrix must be double\n"); + fake_answer(plhs); + return; + } + + if(mxIsStruct(prhs[2])) + { + const char *error_msg; + + // parse options + if(nrhs==4) + { + int i, argc = 1; + char cmd[CMD_LEN], *argv[CMD_LEN/2]; + + // put options in argv[] + mxGetString(prhs[3], cmd, mxGetN(prhs[3]) + 1); + if((argv[argc] = strtok(cmd, " ")) != NULL) + while((argv[++argc] = strtok(NULL, " ")) != NULL) + ; + + for(i=1;i<argc;i++) + { + if(argv[i][0] != '-') break; + if((++i>=argc) && argv[i-1][1] != 'q') + { + exit_with_help(); + fake_answer(plhs); + return; + } + switch(argv[i-1][1]) + { + case 'b': + prob_estimate_flag = atoi(argv[i]); + break; + case 'q': + i--; + info = &print_null; + break; + default: + mexPrintf("Unknown option: -%c\n", argv[i-1][1]); + exit_with_help(); + fake_answer(plhs); + return; + } + } + } + + model = matlab_matrix_to_model(prhs[2], &error_msg); + if (model == NULL) + { + mexPrintf("Error: can't read model: %s\n", error_msg); + fake_answer(plhs); + return; + } + + if(prob_estimate_flag) + { + if(svm_check_probability_model(model)==0) + { + mexPrintf("Model does not support probabiliy estimates\n"); + fake_answer(plhs); + svm_free_and_destroy_model(&model); + return; + } + } + else + { + if(svm_check_probability_model(model)!=0) + info("Model supports probability estimates, but disabled in predicton.\n"); + } + + predict(plhs, prhs, model, prob_estimate_flag); + // destroy model + svm_free_and_destroy_model(&model); + } + else + { + mexPrintf("model file should be a struct array\n"); + fake_answer(plhs); + } + + return; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libsvm_osx64/svmtrain.c Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,474 @@ +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <ctype.h> +#include "../svm.h" + +#include "mex.h" +#include "svm_model_matlab.h" + +#ifdef MX_API_VER +#if MX_API_VER < 0x07030000 +typedef int mwIndex; +#endif +#endif + +#define CMD_LEN 2048 +#define Malloc(type,n) (type *)malloc((n)*sizeof(type)) + +void print_null(const char *s) {} +void print_string_matlab(const char *s) {mexPrintf(s);} + +void exit_with_help() +{ + mexPrintf( + "Usage: model = svmtrain(training_label_vector, training_instance_matrix, 'libsvm_options');\n" + "libsvm_options:\n" + "-s svm_type : set type of SVM (default 0)\n" + " 0 -- C-SVC (multi-class classification)\n" + " 1 -- nu-SVC (multi-class classification)\n" + " 2 -- one-class SVM\n" + " 3 -- epsilon-SVR (regression)\n" + " 4 -- nu-SVR (regression)\n" + "-t kernel_type : set type of kernel function (default 2)\n" + " 0 -- linear: u'*v\n" + " 1 -- polynomial: (gamma*u'*v + coef0)^degree\n" + " 2 -- radial basis function: exp(-gamma*|u-v|^2)\n" + " 3 -- sigmoid: tanh(gamma*u'*v + coef0)\n" + " 4 -- precomputed kernel (kernel values in training_instance_matrix)\n" + "-d degree : set degree in kernel function (default 3)\n" + "-g gamma : set gamma in kernel function (default 1/num_features)\n" + "-r coef0 : set coef0 in kernel function (default 0)\n" + "-c cost : set the parameter C of C-SVC, epsilon-SVR, and nu-SVR (default 1)\n" + "-n nu : set the parameter nu of nu-SVC, one-class SVM, and nu-SVR (default 0.5)\n" + "-p epsilon : set the epsilon in loss function of epsilon-SVR (default 0.1)\n" + "-m cachesize : set cache memory size in MB (default 100)\n" + "-e epsilon : set tolerance of termination criterion (default 0.001)\n" + "-h shrinking : whether to use the shrinking heuristics, 0 or 1 (default 1)\n" + "-b probability_estimates : whether to train a SVC or SVR model for probability estimates, 0 or 1 (default 0)\n" + "-wi weight : set the parameter C of class i to weight*C, for C-SVC (default 1)\n" + "-v n : n-fold cross validation mode\n" + "-q : quiet mode (no outputs)\n" + ); +} + +// svm arguments +struct svm_parameter param; // set by parse_command_line +struct svm_problem prob; // set by read_problem +struct svm_model *model; +struct svm_node *x_space; +int cross_validation; +int nr_fold; + + +double do_cross_validation() +{ + int i; + int total_correct = 0; + double total_error = 0; + double sumv = 0, sumy = 0, sumvv = 0, sumyy = 0, sumvy = 0; + double *target = Malloc(double,prob.l); + double retval = 0.0; + + svm_cross_validation(&prob,¶m,nr_fold,target); + if(param.svm_type == EPSILON_SVR || + param.svm_type == NU_SVR) + { + for(i=0;i<prob.l;i++) + { + double y = prob.y[i]; + double v = target[i]; + total_error += (v-y)*(v-y); + sumv += v; + sumy += y; + sumvv += v*v; + sumyy += y*y; + sumvy += v*y; + } + mexPrintf("Cross Validation Mean squared error = %g\n",total_error/prob.l); + mexPrintf("Cross Validation Squared correlation coefficient = %g\n", + ((prob.l*sumvy-sumv*sumy)*(prob.l*sumvy-sumv*sumy))/ + ((prob.l*sumvv-sumv*sumv)*(prob.l*sumyy-sumy*sumy)) + ); + retval = total_error/prob.l; + } + else + { + for(i=0;i<prob.l;i++) + if(target[i] == prob.y[i]) + ++total_correct; + mexPrintf("Cross Validation Accuracy = %g%%\n",100.0*total_correct/prob.l); + retval = 100.0*total_correct/prob.l; + } + free(target); + return retval; +} + +// nrhs should be 3 +int parse_command_line(int nrhs, const mxArray *prhs[], char *model_file_name) +{ + int i, argc = 1; + char cmd[CMD_LEN]; + char *argv[CMD_LEN/2]; + void (*print_func)(const char *) = print_string_matlab; // default printing to matlab display + + // default values + param.svm_type = C_SVC; + param.kernel_type = RBF; + param.degree = 3; + param.gamma = 0; // 1/num_features + param.coef0 = 0; + param.nu = 0.5; + param.cache_size = 100; + param.C = 1; + param.eps = 1e-3; + param.p = 0.1; + param.shrinking = 1; + param.probability = 0; + param.nr_weight = 0; + param.weight_label = NULL; + param.weight = NULL; + cross_validation = 0; + + if(nrhs <= 1) + return 1; + + if(nrhs > 2) + { + // put options in argv[] + mxGetString(prhs[2], cmd, mxGetN(prhs[2]) + 1); + if((argv[argc] = strtok(cmd, " ")) != NULL) + while((argv[++argc] = strtok(NULL, " ")) != NULL) + ; + } + + // parse options + for(i=1;i<argc;i++) + { + if(argv[i][0] != '-') break; + ++i; + if(i>=argc && argv[i-1][1] != 'q') // since option -q has no parameter + return 1; + switch(argv[i-1][1]) + { + case 's': + param.svm_type = atoi(argv[i]); + break; + case 't': + param.kernel_type = atoi(argv[i]); + break; + case 'd': + param.degree = atoi(argv[i]); + break; + case 'g': + param.gamma = atof(argv[i]); + break; + case 'r': + param.coef0 = atof(argv[i]); + break; + case 'n': + param.nu = atof(argv[i]); + break; + case 'm': + param.cache_size = atof(argv[i]); + break; + case 'c': + param.C = atof(argv[i]); + break; + case 'e': + param.eps = atof(argv[i]); + break; + case 'p': + param.p = atof(argv[i]); + break; + case 'h': + param.shrinking = atoi(argv[i]); + break; + case 'b': + param.probability = atoi(argv[i]); + break; + case 'q': + print_func = &print_null; + i--; + break; + case 'v': + cross_validation = 1; + nr_fold = atoi(argv[i]); + if(nr_fold < 2) + { + mexPrintf("n-fold cross validation: n must >= 2\n"); + return 1; + } + break; + case 'w': + ++param.nr_weight; + param.weight_label = (int *)realloc(param.weight_label,sizeof(int)*param.nr_weight); + param.weight = (double *)realloc(param.weight,sizeof(double)*param.nr_weight); + param.weight_label[param.nr_weight-1] = atoi(&argv[i-1][2]); + param.weight[param.nr_weight-1] = atof(argv[i]); + break; + default: + mexPrintf("Unknown option -%c\n", argv[i-1][1]); + return 1; + } + } + + svm_set_print_string_function(print_func); + + return 0; +} + +// read in a problem (in svmlight format) +int read_problem_dense(const mxArray *label_vec, const mxArray *instance_mat) +{ + int i, j, k; + int elements, max_index, sc, label_vector_row_num; + double *samples, *labels; + + prob.x = NULL; + prob.y = NULL; + x_space = NULL; + + labels = mxGetPr(label_vec); + samples = mxGetPr(instance_mat); + sc = (int)mxGetN(instance_mat); + + elements = 0; + // the number of instance + prob.l = (int)mxGetM(instance_mat); + label_vector_row_num = (int)mxGetM(label_vec); + + if(label_vector_row_num!=prob.l) + { + mexPrintf("Length of label vector does not match # of instances.\n"); + return -1; + } + + if(param.kernel_type == PRECOMPUTED) + elements = prob.l * (sc + 1); + else + { + for(i = 0; i < prob.l; i++) + { + for(k = 0; k < sc; k++) + if(samples[k * prob.l + i] != 0) + elements++; + // count the '-1' element + elements++; + } + } + + prob.y = Malloc(double,prob.l); + prob.x = Malloc(struct svm_node *,prob.l); + x_space = Malloc(struct svm_node, elements); + + max_index = sc; + j = 0; + for(i = 0; i < prob.l; i++) + { + prob.x[i] = &x_space[j]; + prob.y[i] = labels[i]; + + for(k = 0; k < sc; k++) + { + if(param.kernel_type == PRECOMPUTED || samples[k * prob.l + i] != 0) + { + x_space[j].index = k + 1; + x_space[j].value = samples[k * prob.l + i]; + j++; + } + } + x_space[j++].index = -1; + } + + if(param.gamma == 0 && max_index > 0) + param.gamma = 1.0/max_index; + + if(param.kernel_type == PRECOMPUTED) + for(i=0;i<prob.l;i++) + { + if((int)prob.x[i][0].value <= 0 || (int)prob.x[i][0].value > max_index) + { + mexPrintf("Wrong input format: sample_serial_number out of range\n"); + return -1; + } + } + + return 0; +} + +int read_problem_sparse(const mxArray *label_vec, const mxArray *instance_mat) +{ + int i, j, k, low, high; + mwIndex *ir, *jc; + int elements, max_index, num_samples, label_vector_row_num; + double *samples, *labels; + mxArray *instance_mat_col; // transposed instance sparse matrix + + prob.x = NULL; + prob.y = NULL; + x_space = NULL; + + // transpose instance matrix + { + mxArray *prhs[1], *plhs[1]; + prhs[0] = mxDuplicateArray(instance_mat); + if(mexCallMATLAB(1, plhs, 1, prhs, "transpose")) + { + mexPrintf("Error: cannot transpose training instance matrix\n"); + return -1; + } + instance_mat_col = plhs[0]; + mxDestroyArray(prhs[0]); + } + + // each column is one instance + labels = mxGetPr(label_vec); + samples = mxGetPr(instance_mat_col); + ir = mxGetIr(instance_mat_col); + jc = mxGetJc(instance_mat_col); + + num_samples = (int)mxGetNzmax(instance_mat_col); + + // the number of instance + prob.l = (int)mxGetN(instance_mat_col); + label_vector_row_num = (int)mxGetM(label_vec); + + if(label_vector_row_num!=prob.l) + { + mexPrintf("Length of label vector does not match # of instances.\n"); + return -1; + } + + elements = num_samples + prob.l; + max_index = (int)mxGetM(instance_mat_col); + + prob.y = Malloc(double,prob.l); + prob.x = Malloc(struct svm_node *,prob.l); + x_space = Malloc(struct svm_node, elements); + + j = 0; + for(i=0;i<prob.l;i++) + { + prob.x[i] = &x_space[j]; + prob.y[i] = labels[i]; + low = (int)jc[i], high = (int)jc[i+1]; + for(k=low;k<high;k++) + { + x_space[j].index = (int)ir[k] + 1; + x_space[j].value = samples[k]; + j++; + } + x_space[j++].index = -1; + } + + if(param.gamma == 0 && max_index > 0) + param.gamma = 1.0/max_index; + + return 0; +} + +static void fake_answer(mxArray *plhs[]) +{ + plhs[0] = mxCreateDoubleMatrix(0, 0, mxREAL); +} + +// Interface function of matlab +// now assume prhs[0]: label prhs[1]: features +void mexFunction( int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] ) +{ + const char *error_msg; + + // fix random seed to have same results for each run + // (for cross validation and probability estimation) + srand(1); + + // Transform the input Matrix to libsvm format + if(nrhs > 1 && nrhs < 4) + { + int err; + + if(!mxIsDouble(prhs[0]) || !mxIsDouble(prhs[1])) { + mexPrintf("Error: label vector and instance matrix must be double\n"); + fake_answer(plhs); + return; + } + + if(parse_command_line(nrhs, prhs, NULL)) + { + exit_with_help(); + svm_destroy_param(¶m); + fake_answer(plhs); + return; + } + + if(mxIsSparse(prhs[1])) + { + if(param.kernel_type == PRECOMPUTED) + { + // precomputed kernel requires dense matrix, so we make one + mxArray *rhs[1], *lhs[1]; + + rhs[0] = mxDuplicateArray(prhs[1]); + if(mexCallMATLAB(1, lhs, 1, rhs, "full")) + { + mexPrintf("Error: cannot generate a full training instance matrix\n"); + svm_destroy_param(¶m); + fake_answer(plhs); + return; + } + err = read_problem_dense(prhs[0], lhs[0]); + mxDestroyArray(lhs[0]); + mxDestroyArray(rhs[0]); + } + else + err = read_problem_sparse(prhs[0], prhs[1]); + } + else + err = read_problem_dense(prhs[0], prhs[1]); + + // svmtrain's original code + error_msg = svm_check_parameter(&prob, ¶m); + + if(err || error_msg) + { + if (error_msg != NULL) + mexPrintf("Error: %s\n", error_msg); + svm_destroy_param(¶m); + free(prob.y); + free(prob.x); + free(x_space); + fake_answer(plhs); + return; + } + + if(cross_validation) + { + double *ptr; + plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL); + ptr = mxGetPr(plhs[0]); + ptr[0] = do_cross_validation(); + } + else + { + int nr_feat = (int)mxGetN(prhs[1]); + const char *error_msg; + model = svm_train(&prob, ¶m); + error_msg = model_to_matlab_structure(plhs, nr_feat, model); + if(error_msg) + mexPrintf("Error: can't convert libsvm model to matrix structure: %s\n", error_msg); + svm_free_and_destroy_model(&model); + } + svm_destroy_param(¶m); + free(prob.y); + free(prob.x); + free(x_space); + } + else + { + exit_with_help(); + fake_answer(plhs); + return; + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/loadClassificationOutput.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,27 @@ +% from https://soundsoftware.ac.uk/projects/aasp-d-case-metrics +function [fileID,classID] = loadClassificationOutput(filename) + +% Open raw file +fid = fopen(filename,'r+'); + +% Read 1st line +tline = fgetl(fid); +fileID{1} = char(sscanf(tline, '%s\t%*s')); +classID{1} = char(sscanf(tline, '%*s\t%s')); + +% Read rest of the lines +i=1; +while ischar(tline) + i = i+1; + tline = fgetl(fid); + if (ischar(tline)) + + fileID{i} = char(sscanf(tline, '%s\t%*s')); + classID{i} = char(sscanf(tline, '%*s\t%s')); + + end; +end + +% Close file +fclose(fid); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/make_list_files.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,48 @@ +% Copyright 2013 MUSIC TECHNOLOGY GROUP, UNIVERSITAT POMPEU FABRA +% +% Written by Gerard Roma <gerard.roma@upf.edu> +% +% 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 <http://www.gnu.org/licenses/>. + +function make_list_files + NUM_FOLDS = 5; + classes = {'bus' 'busystreet' 'office' 'openairmarket' 'park' 'quietstreet' 'restaurant' 'supermarket' 'tube' 'tubestation'}; + [names,labels] = get_filenames('path_to_files'); + cp = cvpartition(labels,'k',NUM_FOLDS); + for i = 1:NUM_FOLDS + tr_fnames = names(cp.training(i)); + tr_classes = labels(cp.training(i)); + te_fnames = names(cp.test(i)); + te_classes = labels(cp.test(i)); + train_filename = strcat('fold',num2str(i),'_train.txt'); + train_fid = fopen(train_filename,'w+'); + for j = 1:length(tr_fnames) + fprintf(train_fid,'%s\t',char(tr_fnames(j))); + fprintf(train_fid,'%s\n',char(classes(tr_classes(j)))); + end + fclose(train_fid); + + test_filename = strcat('fold',num2str(i),'_test.txt'); + test_fid = fopen(test_filename,'w+'); + gt_filename = strcat('fold',num2str(i),'_gt.txt'); + gt_fid = fopen(gt_filename,'w+'); + for j = 1:length(te_fnames) + fprintf(test_fid,'%s\n',[char(te_fnames(j))]); + fprintf(gt_fid,'%s\t',[char(te_fnames(j))]); + fprintf(gt_fid,'%s\n',[char(classes(te_classes(j)))]); + end + fclose(gt_fid); + fclose(test_fid); + + end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/audspec.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,40 @@ +function [aspectrum,wts] = audspec(pspectrum, sr, nfilts, fbtype, minfreq, maxfreq, sumpower, bwidth) +%[aspectrum,wts] = audspec(pspectrum, sr, nfilts, fbtype, minfreq, maxfreq, sumpower, bwidth) +% +% perform critical band analysis (see PLP) +% takes power spectrogram as input + +if nargin < 2; sr = 16000; end +if nargin < 3; nfilts = ceil(hz2bark(sr/2))+1; end +if nargin < 4; fbtype = 'bark'; end +if nargin < 5; minfreq = 0; end +if nargin < 6; maxfreq = sr/2; end +if nargin < 7; sumpower = 1; end +if nargin < 8; bwidth = 1.0; end + +[nfreqs,nframes] = size(pspectrum); + +nfft = (nfreqs-1)*2; + +if strcmp(fbtype, 'bark') + wts = fft2barkmx(nfft, sr, nfilts, bwidth, minfreq, maxfreq); +elseif strcmp(fbtype, 'mel') + wts = fft2melmx(nfft, sr, nfilts, bwidth, minfreq, maxfreq); +elseif strcmp(fbtype, 'htkmel') + wts = fft2melmx(nfft, sr, nfilts, bwidth, minfreq, maxfreq, 1, 1); +elseif strcmp(fbtype, 'fcmel') + wts = fft2melmx(nfft, sr, nfilts, bwidth, minfreq, maxfreq, 1, 0); +else + disp(['fbtype ', fbtype, ' not recognized']); + error; +end + +wts = wts(:, 1:nfreqs); + +% Integrate FFT bins into Mel bins, in abs or abs^2 domains: +if (sumpower) + aspectrum = wts * pspectrum; +else + aspectrum = (wts * sqrt(pspectrum)).^2; +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/bark2hz.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,18 @@ +function hz=bark2hz(z) +% BARK2HZ Converts frequencies Bark to Hertz (Hz) +% function hz=bark2hz(z) +% Traunmueller-formula for z > 2 Bark +% linear mapping for z <= 2 Bark +% +% Author: Kyrill, Oct. 1996 +% Kyrill, March 1997 (linear mapping added) +% + +%hz_gt_2 = 1960 .* (z + 0.53) ./ (26.28 - z); +%hz_le_2 = z .* 102.9; +% +%hz = (z>2) .* hz_gt_2 + (z<=2) .* hz_le_2; + +% Hynek's formula (2003-04-11 dpwe@ee.columbia.edu) +% (taken from rasta/audspec.c) +hz = 600 * sinh(z/6);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/cep2spec.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,45 @@ +function [spec,idctm] = cep2spec(cep, nfreq, type) +% spec = cep2spec(cep, nfreq, type) +% Reverse the cepstrum to recover a spectrum. +% i.e. converse of spec2cep +% nfreq is how many points to reconstruct in spec +% 2005-05-15 dpwe@ee.columbia.edu + +if nargin < 2; nfreq = 21; end +if nargin < 3; type = 2; end % type of DCT + +[ncep,ncol] = size(cep); + +% Make the DCT matrix +dctm = zeros(ncep, nfreq); +idctm = zeros(nfreq, ncep); +if type == 2 || type == 3 + % this is the orthogonal one, so inv matrix is same as fwd matrix + for i = 1:ncep + dctm(i,:) = cos((i-1)*[1:2:(2*nfreq-1)]/(2*nfreq)*pi) * sqrt(2/nfreq); + end + if type == 2 + % make it unitary! (but not for HTK type 3) + dctm(1,:) = dctm(1,:)/sqrt(2); + else + dctm(1,:) = dctm(1,:)/2; + end + idctm = dctm'; +elseif type == 4 % type 1 with implicit repetition of first, last bins + % so all we do is reconstruct the middle nfreq rows of an nfreq+2 row idctm + for i = 1:ncep + % 2x to compensate for fact that only getting +ve freq half + idctm(:,i) = 2*cos((i-1)*[1:nfreq]'/(nfreq+1)*pi); + end + % fixup 'non-repeated' basis fns + idctm(:, [1 ncep]) = idctm(:, [1 ncep])/2; +else % dpwe type 1 - idft of cosine terms + for i = 1:ncep + % 2x to compensate for fact that only getting +ve freq half + idctm(:,i) = 2*cos((i-1)*[0:(nfreq-1)]'/(nfreq-1)*pi); + end + % fixup 'non-repeated' basis fns + idctm(:, [1 ncep]) = 0.5* idctm(:, [1 ncep]); +end + +spec = exp(idctm*cep);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/deltas.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,35 @@ +function d = deltas(x, w) +% D = deltas(X,W) Calculate the deltas (derivatives) of a sequence +% Use a W-point window (W odd, default 9) to calculate deltas using a +% simple linear slope. This mirrors the delta calculation performed +% in feacalc etc. Each row of X is filtered separately. +% 2003-06-30 dpwe@ee.columbia.edu + +if nargin < 2 + w = 9; +end + +[nr,nc] = size(x); + +if nc == 0 + % empty vector passed in; return empty vector + d = x; + +else + % actually calculate deltas + + % Define window shape + hlen = floor(w/2); + w = 2*hlen + 1; + win = hlen:-1:-hlen; + + % pad data by repeating first and last columns + xx = [repmat(x(:,1),1,hlen),x,repmat(x(:,end),1,hlen)]; + + % Apply the delta filter + d = filter(win, 1, xx, [], 2); % filter along dim 2 (rows) + + % Trim edges + d = d(:,2*hlen + [1:nc]); + +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/dolpc.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,30 @@ +function y = dolpc(x,modelorder) +%y = dolpc(x,modelorder) +% +% compute autoregressive model from spectral magnitude samples +% +% rows(x) = critical band +% col(x) = frame +% +% row(y) = lpc a_i coeffs, scaled by gain +% col(y) = frame +% +% modelorder is order of model, defaults to 8 +% 2003-04-12 dpwe@ee.columbia.edu after shire@icsi.berkeley.edu + +[nbands,nframes] = size(x); + +if nargin < 2 + modelorder = 8; +end + +% Calculate autocorrelation +r = real(ifft([x;x([(nbands-1):-1:2],:)])); +% First half only +r = r(1:nbands,:); + +% Find LPC coeffs by durbin +[y,e] = levinson(r, modelorder); + +% Normalize each poly by gain +y = y'./repmat(e',(modelorder+1),1);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/fft2barkmx.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,42 @@ +function wts = fft2barkmx(nfft, sr, nfilts, width, minfreq, maxfreq) +% wts = fft2barkmx(nfft, sr, nfilts, width, minfreq, maxfreq) +% Generate a matrix of weights to combine FFT bins into Bark +% bins. nfft defines the source FFT size at sampling rate sr. +% Optional nfilts specifies the number of output bands required +% (else one per bark), and width is the constant width of each +% band in Bark (default 1). +% While wts has nfft columns, the second half are all zero. +% Hence, Bark spectrum is fft2barkmx(nfft,sr)*abs(fft(xincols,nfft)); +% 2004-09-05 dpwe@ee.columbia.edu based on rastamat/audspec.m + +if nargin < 3; nfilts = 0; end +if nargin < 4; width = 1.0; end +if nargin < 5; minfreq = 0; end +if nargin < 6; maxfreq = sr/2; end + +min_bark = hz2bark(minfreq); +nyqbark = hz2bark(maxfreq) - min_bark; +if nfilts == 0 + nfilts = ceil(nyqbark)+1; +end + +wts = zeros(nfilts, nfft); + +% bark per filt +step_barks = nyqbark/(nfilts-1); + +% Frequency of each FFT bin in Bark +binbarks = hz2bark([0:nfft/2]*sr/nfft); + +for i = 1:nfilts + f_bark_mid = min_bark + (i-1) * step_barks; + % Linear slopes in log-space (i.e. dB) intersect to trapezoidal window + lof = (binbarks - f_bark_mid - 0.5); + hif = (binbarks - f_bark_mid + 0.5); + wts(i,1:(nfft/2+1)) = 10.^(min(0, min([hif; -2.5*lof])/width)); +end + +function z = hz2bark(f) +% HZ2BARK Converts frequencies Hertz (Hz) to Bark +% taken from rastamat +z = 6 * asinh(f/600);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/fft2melmx.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,140 @@ +function [wts,binfrqs] = fft2melmx(nfft, sr, nfilts, width, minfrq, maxfrq, htkmel, constamp) +% [wts,frqs] = fft2melmx(nfft, sr, nfilts, width, minfrq, maxfrq, htkmel, constamp) +% Generate a matrix of weights to combine FFT bins into Mel +% bins. nfft defines the source FFT size at sampling rate sr. +% Optional nfilts specifies the number of output bands required +% (else one per "mel/width"), and width is the constant width of each +% band relative to standard Mel (default 1). +% While wts has nfft columns, the second half are all zero. +% Hence, Mel spectrum is fft2melmx(nfft,sr)*abs(fft(xincols,nfft)); +% minfrq is the frequency (in Hz) of the lowest band edge; +% default is 0, but 133.33 is a common standard (to skip LF). +% maxfrq is frequency in Hz of upper edge; default sr/2. +% You can exactly duplicate the mel matrix in Slaney's mfcc.m +% as fft2melmx(512, 8000, 40, 1, 133.33, 6855.5, 0); +% htkmel=1 means use HTK's version of the mel curve, not Slaney's. +% constamp=1 means make integration windows peak at 1, not sum to 1. +% frqs returns bin center frqs. +% 2004-09-05 dpwe@ee.columbia.edu based on fft2barkmx + +if nargin < 2; sr = 8000; end +if nargin < 3; nfilts = 0; end +if nargin < 4; width = 1.0; end +if nargin < 5; minfrq = 0; end % default bottom edge at 0 +if nargin < 6; maxfrq = sr/2; end % default top edge at nyquist +if nargin < 7; htkmel = 0; end +if nargin < 8; constamp = 0; end + +if nfilts == 0 + nfilts = ceil(hz2mel(maxfrq, htkmel)/2); +end + +wts = zeros(nfilts, nfft); + +% Center freqs of each FFT bin +fftfrqs = [0:(nfft/2)]/nfft*sr; + +% 'Center freqs' of mel bands - uniformly spaced between limits +minmel = hz2mel(minfrq, htkmel); +maxmel = hz2mel(maxfrq, htkmel); +binfrqs = mel2hz(minmel+[0:(nfilts+1)]/(nfilts+1)*(maxmel-minmel), htkmel); + +binbin = round(binfrqs/sr*(nfft-1)); + +for i = 1:nfilts +% fs = mel2hz(i + [-1 0 1], htkmel); + fs = binfrqs(i+[0 1 2]); + % scale by width + fs = fs(2)+width*(fs - fs(2)); + % lower and upper slopes for all bins + loslope = (fftfrqs - fs(1))/(fs(2) - fs(1)); + hislope = (fs(3) - fftfrqs)/(fs(3) - fs(2)); + % .. then intersect them with each other and zero +% wts(i,:) = 2/(fs(3)-fs(1))*max(0,min(loslope, hislope)); + wts(i,1+[0:(nfft/2)]) = max(0,min(loslope, hislope)); + + % actual algo and weighting in feacalc (more or less) +% wts(i,:) = 0; +% ww = binbin(i+2)-binbin(i); +% usl = binbin(i+1)-binbin(i); +% wts(i,1+binbin(i)+[1:usl]) = 2/ww * [1:usl]/usl; +% dsl = binbin(i+2)-binbin(i+1); +% wts(i,1+binbin(i+1)+[1:(dsl-1)]) = 2/ww * [(dsl-1):-1:1]/dsl; +% need to disable weighting below if you use this one + +end + +if (constamp == 0) + % Slaney-style mel is scaled to be approx constant E per channel + wts = diag(2./(binfrqs(2+[1:nfilts])-binfrqs(1:nfilts)))*wts; +end + +% Make sure 2nd half of FFT is zero +wts(:,(nfft/2+2):nfft) = 0; +% seems like a good idea to avoid aliasing + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function f = mel2hz(z, htk) +% f = mel2hz(z, htk) +% Convert 'mel scale' frequencies into Hz +% Optional htk = 1 means use the HTK formula +% else use the formula from Slaney's mfcc.m +% 2005-04-19 dpwe@ee.columbia.edu + +if nargin < 2 + htk = 0; +end + +if htk == 1 + f = 700*(10.^(z/2595)-1); +else + + f_0 = 0; % 133.33333; + f_sp = 200/3; % 66.66667; + brkfrq = 1000; + brkpt = (brkfrq - f_0)/f_sp; % starting mel value for log region + logstep = exp(log(6.4)/27); % the magic 1.0711703 which is the ratio needed to get from 1000 Hz to 6400 Hz in 27 steps, and is *almost* the ratio between 1000 Hz and the preceding linear filter center at 933.33333 Hz (actually 1000/933.33333 = 1.07142857142857 and exp(log(6.4)/27) = 1.07117028749447) + + linpts = (z < brkpt); + + f = 0*z; + + % fill in parts separately + f(linpts) = f_0 + f_sp*z(linpts); + f(~linpts) = brkfrq*exp(log(logstep)*(z(~linpts)-brkpt)); + +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function z = hz2mel(f,htk) +% z = hz2mel(f,htk) +% Convert frequencies f (in Hz) to mel 'scale'. +% Optional htk = 1 uses the mel axis defined in the HTKBook +% otherwise use Slaney's formula +% 2005-04-19 dpwe@ee.columbia.edu + +if nargin < 2 + htk = 0; +end + +if htk == 1 + z = 2595 * log10(1+f/700); +else + % Mel fn to match Slaney's Auditory Toolbox mfcc.m + + f_0 = 0; % 133.33333; + f_sp = 200/3; % 66.66667; + brkfrq = 1000; + brkpt = (brkfrq - f_0)/f_sp; % starting mel value for log region + logstep = exp(log(6.4)/27); % the magic 1.0711703 which is the ratio needed to get from 1000 Hz to 6400 Hz in 27 steps, and is *almost* the ratio between 1000 Hz and the preceding linear filter center at 933.33333 Hz (actually 1000/933.33333 = 1.07142857142857 and exp(log(6.4)/27) = 1.07117028749447) + + linpts = (f < brkfrq); + + z = 0*f; + + % fill in parts separately + z(linpts) = (f(linpts) - f_0)/f_sp; + z(~linpts) = brkpt+(log(f(~linpts)/brkfrq))./log(logstep); + +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/hz2bark.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,22 @@ +function z= hz2bark(f) +% HZ2BARK Converts frequencies Hertz (Hz) to Bark +% function z= hz2bark(f) +% Uses +% Traunmueller-formula for f > 200 Hz +% linear mapping for f <= 200 Hz +% +% Author: Kyrill, Oct. 1996 +% Kyrill, March 1997 (linear mapping added) +% + +%z_gt_200 = 26.81 .* f ./ (1960 + f) - 0.53; +%z_le_200 = f ./ 102.9; +% +%z = (f>200) .* z_gt_200 + (f<=200) .* z_le_200; + +% Inverse of Hynek's formula (see bark2hz) +z = 6 * asinh(f/600); + +% Formula used in rasta/rasta.h +%z = 6 * log(f/600 + sqrt(1+ ((f/600).^2))); +% They are the same! \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/hz2mel.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,31 @@ +function z = hz2mel(f,htk) +% z = hz2mel(f,htk) +% Convert frequencies f (in Hz) to mel 'scale'. +% Optional htk = 1 uses the mel axis defined in the HTKBook +% otherwise use Slaney's formula +% 2005-04-19 dpwe@ee.columbia.edu + +if nargin < 2 + htk = 0; +end + +if htk == 1 + z = 2595 * log10(1+f/700); +else + % Mel fn to match Slaney's Auditory Toolbox mfcc.m + + f_0 = 0; % 133.33333; + f_sp = 200/3; % 66.66667; + brkfrq = 1000; + brkpt = (brkfrq - f_0)/f_sp; % starting mel value for log region + logstep = exp(log(6.4)/27); % the magic 1.0711703 which is the ratio needed to get from 1000 Hz to 6400 Hz in 27 steps, and is *almost* the ratio between 1000 Hz and the preceding linear filter center at 933.33333 Hz (actually 1000/933.33333 = 1.07142857142857 and exp(log(6.4)/27) = 1.07117028749447) + + linpts = (f < brkfrq); + + z = 0*f; + + % fill in parts separately + z(linpts) = (f(linpts) - f_0)/f_sp; + z(~linpts) = brkpt+(log(f(~linpts)/brkfrq))./log(logstep); + +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/invaudspec.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,43 @@ +function [spec,wts,iwts] = invaudspec(aspectrum, sr, nfft, fbtype, minfreq, maxfreq, sumpower, bwidth) +%pspectrum = invaudspec(aspectrum, sr, nfft, fbtype, minfreq, maxfreq, sumpower, bwidth) +% +% Invert (as best we can) the effects of audspec() +% +% 2004-02-04 dpwe@ee.columbia.edu + +if nargin < 2; sr = 16000; end +if nargin < 3; nfft = 512; end +if nargin < 4; fbtype = 'bark'; end +if nargin < 5; minfreq = 0; end +if nargin < 6; maxfreq = sr/2; end +if nargin < 7; sumpower = 1; end +if nargin < 8; bwidth = 1.0; end + +[nfilts,nframes] = size(aspectrum); + +if strcmp(fbtype, 'bark') + wts = fft2barkmx(nfft, sr, nfilts, bwidth, minfreq, maxfreq); +elseif strcmp(fbtype, 'mel') + wts = fft2melmx(nfft, sr, nfilts, bwidth, minfreq, maxfreq); +elseif strcmp(fbtype, 'htkmel') + wts = fft2melmx(nfft, sr, nfilts, bwidth, minfreq, maxfreq, 1, 1); +elseif strcmp(fbtype, 'fcmel') + wts = fft2melmx(nfft, sr, nfilts, bwidth, minfreq, maxfreq, 1); +else + disp(['fbtype ', fbtype, ' not recognized']); + error; +end + +% Cut off 2nd half +wts = wts(:,1:((nfft/2)+1)); + +% Just transpose, fix up +ww = wts'*wts; +iwts = wts'./(repmat(max(mean(diag(ww))/100, sum(ww))',1,nfilts)); +% Apply weights +if (sumpower) + spec = iwts * aspectrum; +else + spec = (iwts * sqrt(aspectrum)).^2; +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/invmelfcc.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,43 @@ +function [x,aspc,spec] = invmelfcc(cep, sr, varargin) +% [x,aspc,spec] = invmelfcc(cep, sr[, opts ...]) +% Attempt to invert plp cepstra back to a full spectrum +% and even a waveform. Takes all the same options as melfcc. +% x is (noise-excited) time domain waveform; aspc is the +% auditory spectrogram, spec is the |STFT| spectrogram. +% 2005-05-15 dpwe@ee.columbia.edu + +% Parse out the optional arguments +[wintime, hoptime, numcep, lifterexp, sumpower, preemph, dither, ... + minfreq, maxfreq, nbands, bwidth, dcttype, fbtype, usecmp, modelorder, broaden, excitation] = ... + process_options(varargin, 'wintime', 0.025, 'hoptime', 0.010, ... + 'numcep', 13, 'lifterexp', 0.6, 'sumpower', 1, 'preemph', 0.97, ... + 'dither', 0, 'minfreq', 0, 'maxfreq', 4000, ... + 'nbands', 40, 'bwidth', 1.0, 'dcttype', 2, ... + 'fbtype', 'mel', 'usecmp', 0, 'modelorder', 0, 'broaden', ... + 0, 'excitation', []); + +winpts = round(wintime*sr); +nfft = 2^(ceil(log(winpts)/log(2))); + +cep = lifter(cep, lifterexp, 1); % 3rd arg nonzero means undo liftering + +% Need to reconstruct the two extra flanking bands for invpostaud to delete +% (if we're doing usecmp) +pspc = cep2spec(cep, nbands+2*broaden, dcttype); + +if (usecmp) + aspc = invpostaud(pspc, maxfreq, fbtype, broaden); +else + aspc = pspc; +end + +% Undo the auditory spectrum +spec = invaudspec(aspc, sr, nfft, fbtype, minfreq, maxfreq, sumpower, bwidth); + +% Back to waveform (modulate white noise, or specified excitation) +x = invpowspec(spec, sr, wintime, hoptime, excitation); + +if preemph ~= 0 + % Undo the original preemphasis + x = filter(1, [1 -preemph], x); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/invpostaud.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,56 @@ +function [x,eql] = invpostaud(y,fmax,fbtype,broaden) +%x = postaud(y,fmax,fbtype) +% +% invert the effects of postaud (loudness equalization and cube +% root compression) +% y = postaud output +% x = reconstructed critical band filters +% rows = critical bands +% cols = frames + +if nargin < 3 + fbtype = 'bark'; +end +if nargin < 4 + % By default, don't add extra flanking bands + broaden = 0; +end + +[nbands,nframes] = size(y); + +% equal loundness weights stolen from rasta code +%eql = [0.000479 0.005949 0.021117 0.044806 0.073345 0.104417 0.137717 ... +% 0.174255 0.215590 0.263260 0.318302 0.380844 0.449798 0.522813 0.596597]; + +if strcmp(fbtype, 'bark') + bandcfhz = bark2hz(linspace(0, hz2bark(fmax), nbands)); +elseif strcmp(fbtype, 'mel') + bandcfhz = mel2hz(linspace(0, hz2mel(fmax), nbands)); +elseif strcmp(fbtype, 'htkmel') || strcmp(fbtype, 'fcmel') + bandcfhz = mel2hz(linspace(0, hz2mel(fmax,1), nbands),1); +else + disp(['unknown fbtype', fbtype]); + error; +end + +% Remove extremal bands (the ones that got duplicated) +bandcfhz = bandcfhz((1+broaden):(nbands-broaden)); + +% Hynek's magic equal-loudness-curve formula +fsq = bandcfhz.^2; +ftmp = fsq + 1.6e5; +eql = ((fsq./ftmp).^2) .* ((fsq + 1.44e6)./(fsq + 9.61e6)); + + +% cube expand +x = y .^ (1/.33); + +%% squash the zero in the eql curve +if eql(1) == 0 % or maybe always + eql(1) = eql(2); + eql(end) = eql(end-1); +end + +% weight the critical bands +x = x((1+broaden):(nbands-broaden),:)./repmat(eql',1,nframes); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/invpowspec.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,64 @@ +function x = invpowspec(y, sr, wintime, steptime, excit) +%x = invpowspec(y, sr, wintime, steptime, excit) +% +% Attempt to go back from specgram-like powerspec to audio waveform +% by scaling specgram of white noise +% +% default values: +% sr = 8000Hz +% wintime = 25ms (200 samps) +% steptime = 10ms (80 samps) +% which means use 256 point fft +% hamming window +% +% excit is input excitation; white noise is used if not specified + +% for sr = 8000 +%NFFT = 256; +%NOVERLAP = 120; +%SAMPRATE = 8000; +%WINDOW = hamming(200); + +[nrow, ncol] = size(y); + +if nargin < 2 + sr = 8000; +end +if nargin < 3 + wintime = 0.025; +end +if nargin < 4 + steptime = 0.010; +end +if nargin < 5 + r = []; +else + r = excit; +end + +winpts = round(wintime*sr); +steppts = round(steptime*sr); + +NFFT = 2^(ceil(log(winpts)/log(2))); + +if NFFT ~= 2*(nrow-1) + disp('Inferred FFT size doesn''t match specgram'); +end + +NOVERLAP = winpts - steppts; +SAMPRATE = sr; + +% Values coming out of rasta treat samples as integers, +% not range -1..1, hence scale up here to match (approx) +%y = abs(specgram(x*32768,NFFT,SAMPRATE,WINDOW,NOVERLAP)).^2; + +xlen = winpts + steppts*(ncol - 1); + +if length(r) == 0 + r = randn(xlen,1); +end +r = r(1:xlen); + +R = specgram(r/32768/12, NFFT, SAMPRATE, winpts, NOVERLAP); +R = R .* sqrt(y);; +x = ispecgram(R, NFFT, SAMPRATE, winpts, NOVERLAP);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/ispecgram.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,50 @@ +function x = ispecgram(d, ftsize, sr, win, nov) +% X = ispecgram(D, F, SR, WIN, NOV) Inverse specgram +% Overlap-add the inverse of the output of specgram +% ftsize is implied by sizeof d, sr is ignored, nov defaults to ftsize/2 +% dpwe 2005may16. after istft + +[nspec,ncol] = size(d); + +if nargin < 2 + ftsize = 2*(nspec-1); +end +if nargin < 3 + % who cares? +end +if nargin < 4 + win = ftsize; % doesn't matter either - assume it added up OK +end +if nargin < 5 + nov = ftsize/2; +end + +hop = win - nov; + +if nspec ~= (ftsize/2)+1 + error('number of rows should be fftsize/2+1') +end + +xlen = ftsize + (ncol-1) * hop; +x = zeros(xlen,1); + +halff = ftsize/2; % midpoint of win + +% No reconstruction win (for now...) + +for c = 1:ncol + ft = d(:,c); + ft = [ft(1:(ftsize/2+1)); conj(ft([(ftsize/2):-1:2]))]; + + if max(imag(ifft(ft))) > 1e-5 + disp('imag oflow'); + end + + px = real(ifft(ft)); % no shift in specgram + + b = (c-1)*hop; + x(b+[1:ftsize]) = x(b+[1:ftsize]) + px; +end; + +x = x * win/ftsize; % scale amplitude +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/lifter.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,38 @@ +function y = lifter(x, lift, invs) +% y = lifter(x, lift, invs) +% Apply lifter to matrix of cepstra (one per column) +% lift = exponent of x i^n liftering +% or, as a negative integer, the length of HTK-style sin-curve liftering. +% If inverse == 1 (default 0), undo the liftering. +% 2005-05-19 dpwe@ee.columbia.edu + +if nargin < 2; lift = 0.6; end % liftering exponent +if nargin < 3; invs = 0; end % flag to undo liftering + +[ncep, nfrm] = size(x); + +if lift == 0 + y = x; +else + + if lift > 0 + if lift > 10 + disp(['Unlikely lift exponent of ', num2str(lift),' (did you mean -ve?)']); + end + liftwts = [1, ([1:(ncep-1)].^lift)]; + elseif lift < 0 + % Hack to support HTK liftering + L = -lift; + if (L ~= round(L)) + disp(['HTK liftering value ', num2str(L),' must be integer']); + end + liftwts = [1, (1+L/2*sin([1:(ncep-1)]*pi/L))]; + end + + if (invs) + liftwts = 1./liftwts; + end + + y = diag(liftwts)*x; + +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/lpc2cep.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,35 @@ +function features = lpc2cep(a,nout) +% features = lpc2cep(lpcas,nout) +% Convert the LPC 'a' coefficients in each column of lpcas +% into frames of cepstra. +% nout is number of cepstra to produce, defaults to size(lpcas,1) +% 2003-04-11 dpwe@ee.columbia.edu + +[nin, ncol] = size(a); + +order = nin - 1; + +if nargin < 2 + nout = order + 1; +end + +c = zeros(nout, ncol); + +% Code copied from HSigP.c: LPC2Cepstrum + +% First cep is log(Error) from Durbin +c(1,:) = -log(a(1,:)); + +% Renormalize lpc A coeffs +a = a ./ repmat(a(1,:), nin, 1); + +for n = 2:nout + sum = 0; + for m = 2:n + sum = sum + (n - m) * a(m,:) .* c(n - m + 1, :); + end + c(n,:) = -(a(n,:) + sum / (n-1)); +end + +features = c; +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/lpc2spec.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,40 @@ +function [features,F,M] = lpc2spec(lpcas, nout) +% [features,F,M] = lpc2spec(lpcas,nout) +% Convert LPC coeffs back into spectra +% nout is number of freq channels, default 17 (i.e. for 8 kHz) +% 2003-04-11 dpwe@ee.columbia.edu part of rastamat + +if nargin < 2 + nout = 17; +end + +[rows, cols] = size(lpcas); +order = rows - 1; + +gg = lpcas(1,:); +aa = lpcas./repmat(gg,rows,1); + +% Calculate the actual z-plane polyvals: nout points around unit circle +zz = exp((-j*[0:(nout-1)]'*pi/(nout-1))*[0:order]); + +% Actual polyvals, in power (mag^2) +features = ((1./abs(zz*aa)).^2)./repmat(gg,nout,1); + +F = zeros(cols, floor(rows/2)); +M = F; + +for c = 1:cols; + aaa = aa(:,c); + rr = roots(aaa'); + ff = angle(rr'); +% size(ff) +% size(aaa) + zz = exp(j*ff'*[0:(length(aaa)-1)]); + mags = sqrt(((1./abs(zz*aaa)).^2)/gg(c))'; + + [dummy, ix] = sort(ff); + keep = ff(ix) > 0; + ix = ix(keep); + F(c,1:length(ix)) = ff(ix); + M(c,1:length(ix)) = mags(ix); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/mel2hz.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,30 @@ +function f = mel2hz(z, htk) +% f = mel2hz(z, htk) +% Convert 'mel scale' frequencies into Hz +% Optional htk = 1 means use the HTK formula +% else use the formula from Slaney's mfcc.m +% 2005-04-19 dpwe@ee.columbia.edu + +if nargin < 2 + htk = 0; +end + +if htk == 1 + f = 700*(10.^(z/2595)-1); +else + + f_0 = 0; % 133.33333; + f_sp = 200/3; % 66.66667; + brkfrq = 1000; + brkpt = (brkfrq - f_0)/f_sp; % starting mel value for log region + logstep = exp(log(6.4)/27); % the magic 1.0711703 which is the ratio needed to get from 1000 Hz to 6400 Hz in 27 steps, and is *almost* the ratio between 1000 Hz and the preceding linear filter center at 933.33333 Hz (actually 1000/933.33333 = 1.07142857142857 and exp(log(6.4)/27) = 1.07117028749447) + + linpts = (z < brkpt); + + f = 0*z; + + % fill in parts separately + f(linpts) = f_0 + f_sp*z(linpts); + f(~linpts) = brkfrq*exp(log(logstep)*(z(~linpts)-brkpt)); + +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/melfcc.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,107 @@ +function [cepstra,aspectrum,pspectrum] = melfcc(samples, sr, varargin) +%[cepstra,aspectrum,pspectrum] = melfcc(samples, sr[, opts ...]) +% Calculate Mel-frequency cepstral coefficients by: +% - take the absolute value of the STFT +% - warp to a Mel frequency scale +% - take the DCT of the log-Mel-spectrum +% - return the first <ncep> components +% This version allows a lot of options to be controlled, as optional +% 'name', value pairs from the 3rd argument on: (defaults in parens) +% 'wintime' (0.025): window length in sec +% 'hoptime' (0.010): step between successive windows in sec +% 'numcep' (13): number of cepstra to return +% 'lifterexp' (0.6): exponent for liftering; 0 = none; < 0 = HTK sin lifter +% 'sumpower' (1): 1 = sum abs(fft)^2; 0 = sum abs(fft) +% 'preemph' (0.97): apply pre-emphasis filter [1 -preemph] (0 = none) +% 'dither' (0): 1 = add offset to spectrum as if dither noise +% 'minfreq' (0): lowest band edge of mel filters (Hz) +% 'maxfreq' (4000): highest band edge of mel filters (Hz) +% 'nbands' (40): number of warped spectral bands to use +% 'bwidth' (1.0): width of aud spec filters relative to default +% 'dcttype' (2): type of DCT used - 1 or 2 (or 3 for HTK or 4 for feac) +% 'fbtype' ('mel'): frequency warp: 'mel','bark','htkmel','fcmel' +% 'usecmp' (0): apply equal-loudness weighting and cube-root compr. +% 'modelorder' (0): if > 0, fit a PLP model of this order +% 'broaden' (0): flag to retain the (useless?) first and last bands +% 'useenergy' (0): overwrite C0 with true log energy +% The following non-default values nearly duplicate Malcolm Slaney's mfcc +% (i.e. melfcc(d,16000,opts...) =~= log(10)*2*mfcc(d*(2^17),16000) ) +% 'wintime': 0.016 +% 'lifterexp': 0 +% 'minfreq': 133.33 +% 'maxfreq': 6855.6 +% 'sumpower': 0 +% The following non-default values nearly duplicate HTK's MFCC +% (i.e. melfcc(d,16000,opts...) =~= 2*htkmelfcc(:,[13,[1:12]])' +% where HTK config has PREEMCOEF = 0.97, NUMCHANS = 20, CEPLIFTER = 22, +% NUMCEPS = 12, WINDOWSIZE = 250000.0, USEHAMMING = T, TARGETKIND = MFCC_0) +% 'lifterexp': -22 +% 'nbands': 20 +% 'maxfreq': 8000 +% 'sumpower': 0 +% 'fbtype': 'htkmel' +% 'dcttype': 3 +% For more detail on reproducing other programs' outputs, see +% http://www.ee.columbia.edu/~dpwe/resources/matlab/rastamat/mfccs.html +% +% 2005-04-19 dpwe@ee.columbia.edu after rastaplp.m. +% Uses Mark Paskin's process_options.m from KPMtools +% $Header: /Users/dpwe/matlab/rastamat/RCS/melfcc.m,v 1.3 2012/09/03 14:01:26 dpwe Exp dpwe $ + +if nargin < 2; sr = 16000; end + +% Parse out the optional arguments +[wintime, hoptime, numcep, lifterexp, sumpower, preemph, dither, ... + minfreq, maxfreq, nbands, bwidth, dcttype, fbtype, usecmp, modelorder, ... + broaden, useenergy] = ... + process_options(varargin, 'wintime', 0.025, 'hoptime', 0.010, ... + 'numcep', 13, 'lifterexp', 0.6, 'sumpower', 1, 'preemph', 0.97, ... + 'dither', 0, 'minfreq', 0, 'maxfreq', 4000, ... + 'nbands', 40, 'bwidth', 1.0, 'dcttype', 2, ... + 'fbtype', 'mel', 'usecmp', 0, 'modelorder', 0, ... + 'broaden', 0, 'useenergy', 0); + +if preemph ~= 0 + samples = filter([1 -preemph], 1, samples); +end + +% Compute FFT power spectrum +[pspectrum,logE] = powspec(samples, sr, wintime, hoptime, dither); + +aspectrum = audspec(pspectrum, sr, nbands, fbtype, minfreq, maxfreq, sumpower, bwidth); + +if (usecmp) + % PLP-like weighting/compression + aspectrum = postaud(aspectrum, maxfreq, fbtype, broaden); +end + +if modelorder > 0 + + if (dcttype ~= 1) + disp(['warning: plp cepstra are implicitly dcttype 1 (not ', num2str(dcttype), ')']); + end + + % LPC analysis + lpcas = dolpc(aspectrum, modelorder); + + % convert lpc to cepstra + cepstra = lpc2cep(lpcas, numcep); + + % Return the auditory spectrum corresponding to the cepstra? +% aspectrum = lpc2spec(lpcas, nbands); + % else return the aspectrum that the cepstra are based on, prior to PLP + +else + + % Convert to cepstra via DCT + cepstra = spec2cep(aspectrum, numcep, dcttype); + +end + +cepstra = lifter(cepstra, lifterexp); + +if useenergy + cepstra(1,:) = logE; +end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/postaud.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,60 @@ +function [y,eql] = postaud(x,fmax,fbtype,broaden) +%y = postaud(x, fmax, fbtype) +% +% do loudness equalization and cube root compression +% x = critical band filters +% rows = critical bands +% cols = frames + +if nargin < 3 + fbtype = 'bark'; +end +if nargin < 4 + % By default, don't add extra flanking bands + broaden = 0; +end + + +[nbands,nframes] = size(x); + +% equal loundness weights stolen from rasta code +%eql = [0.000479 0.005949 0.021117 0.044806 0.073345 0.104417 0.137717 ... +% 0.174255 0.215590 0.263260 0.318302 0.380844 0.449798 0.522813 +% 0.596597]; + +% Include frequency points at extremes, discard later +nfpts = nbands+2*broaden; + +if strcmp(fbtype, 'bark') + bandcfhz = bark2hz(linspace(0, hz2bark(fmax), nfpts)); +elseif strcmp(fbtype, 'mel') + bandcfhz = mel2hz(linspace(0, hz2mel(fmax), nfpts)); +elseif strcmp(fbtype, 'htkmel') || strcmp(fbtype, 'fcmel') + bandcfhz = mel2hz(linspace(0, hz2mel(fmax,1), nfpts),1); +else + disp(['unknown fbtype', fbtype]); + error; +end + +% Remove extremal bands (the ones that will be duplicated) +bandcfhz = bandcfhz((1+broaden):(nfpts-broaden)); + +% Hynek's magic equal-loudness-curve formula +fsq = bandcfhz.^2; +ftmp = fsq + 1.6e5; +eql = ((fsq./ftmp).^2) .* ((fsq + 1.44e6)./(fsq + 9.61e6)); + +% weight the critical bands +z = repmat(eql',1,nframes).*x; + +% cube root compress +z = z .^ (.33); + +% replicate first and last band (because they are unreliable as calculated) +if (broaden) + y = z([1,1:nbands,nbands],:); +else + y = z([2,2:(nbands-1),nbands-1],:); +end +%y = z([1,1:nbands-2,nbands-2],:); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/powspec.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,65 @@ +function [y,e] = powspec(x, sr, wintime, steptime, dither) +%[y,e] = powspec(x, sr, wintime, steptime, sumlin, dither) +% +% compute the powerspectrum and frame energy of the input signal. +% basically outputs a power spectrogram +% +% each column represents a power spectrum for a given frame +% each row represents a frequency +% +% default values: +% sr = 8000Hz +% wintime = 25ms (200 samps) +% steptime = 10ms (80 samps) +% which means use 256 point fft +% hamming window +% +% $Header: /Users/dpwe/matlab/rastamat/RCS/powspec.m,v 1.3 2012/09/03 14:02:01 dpwe Exp dpwe $ + +% for sr = 8000 +%NFFT = 256; +%NOVERLAP = 120; +%SAMPRATE = 8000; +%WINDOW = hamming(200); + +if nargin < 2 + sr = 8000; +end +if nargin < 3 + wintime = 0.025; +end +if nargin < 4 + steptime = 0.010; +end +if nargin < 5 + dither = 1; +end + +winpts = round(wintime*sr); +steppts = round(steptime*sr); + +NFFT = 2^(ceil(log(winpts)/log(2))); +%WINDOW = hamming(winpts); +%WINDOW = [0,hanning(winpts)']; +WINDOW = [hanning(winpts)]; +% hanning gives much less noisy sidelobes +NOVERLAP = winpts - steppts; +SAMPRATE = sr; + +% Values coming out of rasta treat samples as integers, +% not range -1..1, hence scale up here to match (approx) +y = abs(specgram(x*32768,NFFT,SAMPRATE,WINDOW,NOVERLAP)).^2; + +% imagine we had random dither that had a variance of 1 sample +% step and a white spectrum. That's like (in expectation, anyway) +% adding a constant value to every bin (to avoid digital zero) +if (dither) + y = y + winpts; +end +% ignoring the hamming window, total power would be = #pts +% I think this doesn't quite make sense, but it's what rasta/powspec.c does + +% that's all she wrote + +% 2012-09-03 Calculate log energy - after windowing, by parseval +e = log(sum(y));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/process_options.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,146 @@ +% PROCESS_OPTIONS - Processes options passed to a Matlab function. +% This function provides a simple means of +% parsing attribute-value options. Each option is +% named by a unique string and is given a default +% value. +% +% Usage: [var1, var2, ..., varn[, unused]] = ... +% process_options(args, ... +% str1, def1, str2, def2, ..., strn, defn) +% +% Arguments: +% args - a cell array of input arguments, such +% as that provided by VARARGIN. Its contents +% should alternate between strings and +% values. +% str1, ..., strn - Strings that are associated with a +% particular variable +% def1, ..., defn - Default values returned if no option +% is supplied +% NOTE: this version (modified by dpwe) will run str2num on +% string arguments if the default value is numeric +% (to support passing in numbers from the command line). +% +% Returns: +% var1, ..., varn - values to be assigned to variables +% unused - an optional cell array of those +% string-value pairs that were unused; +% if this is not supplied, then a +% warning will be issued for each +% option in args that lacked a match. +% +% Examples: +% +% Suppose we wish to define a Matlab function 'func' that has +% required parameters x and y, and optional arguments 'u' and 'v'. +% With the definition +% +% function y = func(x, y, varargin) +% +% [u, v] = process_options(varargin, 'u', 0, 'v', 1); +% +% calling func(0, 1, 'v', 2) will assign 0 to x, 1 to y, 0 to u, and 2 +% to v. The parameter names are insensitive to case; calling +% func(0, 1, 'V', 2) has the same effect. The function call +% +% func(0, 1, 'u', 5, 'z', 2); +% +% will result in u having the value 5 and v having value 1, but +% will issue a warning that the 'z' option has not been used. On +% the other hand, if func is defined as +% +% function y = func(x, y, varargin) +% +% [u, v, unused_args] = process_options(varargin, 'u', 0, 'v', 1); +% +% then the call func(0, 1, 'u', 5, 'z', 2) will yield no warning, +% and unused_args will have the value {'z', 2}. This behaviour is +% useful for functions with options that invoke other functions +% with options; all options can be passed to the outer function and +% its unprocessed arguments can be passed to the inner function. + +% Copyright (C) 2002 Mark A. Paskin +% +% This program is free software; you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation; either version 2 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 +% General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program; if not, write to the Free Software +% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +% USA. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function [varargout] = process_options(args, varargin) + +% Check the number of input arguments +n = length(varargin); +if (mod(n, 2)) + error('Each option must be a string/value pair.'); +end + +% Check the number of supplied output arguments +if (nargout < (n / 2)) + error('Insufficient number of output arguments given'); +elseif (nargout == (n / 2)) + warn = 1; + nout = n / 2; +else + warn = 0; + nout = n / 2 + 1; +end + +% Set outputs to be defaults +varargout = cell(1, nout); +for i=2:2:n + varargout{i/2} = varargin{i}; +end + +% Now process all arguments +nunused = 0; +for i=1:2:length(args) + found = 0; + for j=1:2:n + if strcmpi(args{i}, varargin{j}) + % dpwe: promote string args to numeric if default arg is + % numeric + if isnumeric(varargin{j+1}) && ischar(args{i+1}) + varargout{(j + 1)/2} = str2num(args{i + 1}); + % this occurs when specifying arguments from command line + else + % in all other cases, take what you get + varargout{(j + 1)/2} = args{i + 1}; + end + found = 1; + break; + end + end + if (~found) + if (warn) + warning(sprintf('Option ''%s'' not used.', args{i})); + args{i} + else + nunused = nunused + 1; + unused{2 * nunused - 1} = args{i}; + if length(args) > i + % don't demand a value for the last, unused tag (e.g. -help) + unused{2 * nunused} = args{i + 1}; + end + end + end +end + +% Assign the unused arguments +if (~warn) + if (nunused) + varargout{nout} = unused; + else + varargout{nout} = cell(0); + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/rastafilt.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,31 @@ +function y = rastafilt(x) +% y = rastafilt(x) +% +% rows of x = critical bands, cols of x = frame +% same for y but after filtering +% +% default filter is single pole at 0.94 + +% rasta filter +numer = [-2:2]; +numer = -numer ./ sum(numer.*numer); +denom = [1 -0.94]; + +% Initialize the state. This avoids a big spike at the beginning +% resulting from the dc offset level in each band. +% (this is effectively what rasta/rasta_filt.c does). +% Because Matlab uses a DF2Trans implementation, we have to +% specify the FIR part to get the state right (but not the IIR part) +%[y,z] = filter(numer, 1, x(:,1:4)'); +% 2010-08-12 filter() chooses the wrong dimension for very short +% signals (thanks to Benjamin K otric1@gmail.com) +[y,z] = filter(numer, 1, x(:,1:4)',[],1); +% .. but don't keep any of these values, just output zero at the beginning +y = 0*y'; + +size(z) +size(x) + +% Apply the full filter to the rest of the signal, append it +y = [y,filter(numer, denom, x(:,5:end)',z,1)']; +%y = [y,filter(numer, denom, x(:,5:end)',z)'];
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/rastaplp.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,74 @@ +function [cepstra, spectra, pspectrum, lpcas, F, M] = rastaplp(samples, sr, dorasta, modelorder) +%[cepstra, spectra, lpcas] = rastaplp(samples, sr, dorasta, modelorder) +% +% cheap version of log rasta with fixed parameters +% +% output is matrix of features, row = feature, col = frame +% +% sr is sampling rate of samples, defaults to 8000 +% dorasta defaults to 1; if 0, just calculate PLP +% modelorder is order of PLP model, defaults to 8. 0 -> no PLP +% +% rastaplp(d, sr, 0, 12) is pretty close to the unix command line +% feacalc -dith -delta 0 -ras no -plp 12 -dom cep ... +% except during very quiet areas, where our approach of adding noise +% in the time domain is different from rasta's approach +% +% 2003-04-12 dpwe@ee.columbia.edu after shire@icsi.berkeley.edu's version + +if nargin < 2 + sr = 8000; +end +if nargin < 3 + dorasta = 1; +end +if nargin < 4 + modelorder = 8; +end + +% add miniscule amount of noise +%samples = samples + randn(size(samples))*0.0001; + +% first compute power spectrum +pspectrum = powspec(samples, sr); + +% next group to critical bands +aspectrum = audspec(pspectrum, sr); +nbands = size(aspectrum,1); + +if dorasta ~= 0 + + % put in log domain + nl_aspectrum = log(aspectrum); + + % next do rasta filtering + ras_nl_aspectrum = rastafilt(nl_aspectrum); + + % do inverse log + aspectrum = exp(ras_nl_aspectrum); + +end + +% do final auditory compressions +postspectrum = postaud(aspectrum, sr/2); % 2012-09-03 bug: was sr + +if modelorder > 0 + + % LPC analysis + lpcas = dolpc(postspectrum, modelorder); + + % convert lpc to cepstra + cepstra = lpc2cep(lpcas, modelorder+1); + + % .. or to spectra + [spectra,F,M] = lpc2spec(lpcas, nbands); + +else + + % No LPC smoothing of spectrum + spectra = postspectrum; + cepstra = spec2cep(spectra); + +end + +cepstra = lifter(cepstra, 0.6);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/readhtk.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,92 @@ +function [d,fp,dt,tc,t]=readhtk(file) +%READHTK read an HTK parameter file [D,FP,DT,TC,T]=(FILE) +% +% d is data, fp is frame period in seconds +% dt is data type, tc is full type code, t is a text version of the full typecode +% tc is the sum of the following values: +% 0 WAVEFORM +% 1 LPC +% 2 LPREFC +% 3 LPCEPSTRA +% 4 LPDELCEP +% 5 IREFC +% 6 MFCC +% 7 FBANK +% 8 MELSPEC +% 9 USER +% 10 DISCRETE +% 11 PLP +% 64 -E Includes energy terms +% 128 _N Suppress absolute energy +% 256 _D Include delta coefs +% 512 _A Include acceleration coefs +% 1024 _C Compressed +% 2048 _Z Zero mean static coefs +% 4096 _K CRC checksum (not implemented yet) +% 8192 _0 Include 0'th cepstral coef + + +% Copyright (C) Mike Brookes 1997 +% +% This version modified to read HTK's compressed feature files +% 2005-05-18 dpwe@ee.columbia.edu +% +% VOICEBOX home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% This program is free software; you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation; either version 2 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 General Public License for more details. +% +% You can obtain a copy of the GNU General Public License from +% ftp://prep.ai.mit.edu/pub/gnu/COPYING-2.0 or by writing to +% Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +fid=fopen(file,'r','b'); +if fid < 0 + error(sprintf('Cannot read from file %s',file)); +end +nf=fread(fid,1,'long'); +fp=fread(fid,1,'long')*1.E-7; +by=fread(fid,1,'short'); +tc=fread(fid,1,'short'); +hb=floor(tc*pow2(-14:-6)); +hd=hb(9:-1:2)-2*hb(8:-1:1); +dt=tc-64*hb(9); + +% hd(7)=1 CRC check +% hd(5)=1 compressed data + +if ( dt == 0 ), + d=fread(fid,Inf,'short'); +else + if (hd(5) == 1) + % compressed data - first read scales + nf = nf - 4; + ncol = by / 2; + scales = fread(fid, ncol, 'float'); + biases = fread(fid, ncol, 'float'); + d = fread(fid,[ncol, nf], 'short'); + d = repmat(1./scales,1,nf).*(d+repmat(biases,1,nf)); + d = d.'; + else + d=fread(fid,[by/4,nf],'float').'; + end +end; +fclose(fid); +if nargout > 2 + hd(7)=0; + hd(5)=0; + ns=sum(hd); + kinds=['WAVEFORM ';'LPC ';'LPREFC ';'LPCEPSTRA ';'LPDELCEP ';'IREFC ';'MFCC ';'FBANK ';'MELSPEC ';'USER ';'DISCRETE ';'PLP ';'??? ']; + kind=kinds(min(dt+1,size(kinds,1)),:); + cc='ENDACZK0'; + t=[kind(1:min(find(kind==' '))-1) reshape(['_'*ones(1,ns);cc(hd>0)],1,2*ns)]; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rastamat/spec2cep.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,51 @@ +function [cep,dctm] = spec2cep(spec, ncep, type) +% [cep,dctm] = spec2cep(spec, ncep, type) +% Calculate cepstra from spectral samples (in columns of spec) +% Return ncep cepstral rows (defaults to 9) +% This one does type II dct, or type I if type is specified as 1 +% dctm returns the DCT matrix that spec was multiplied by to give cep. +% 2005-04-19 dpwe@ee.columbia.edu for mfcc_dpwe + +if nargin < 2; ncep = 13; end +if nargin < 3; type = 2; end % type of DCT + +[nrow, ncol] = size(spec); + +% Make the DCT matrix +dctm = zeros(ncep, nrow); +if type == 2 || type == 3 + % this is the orthogonal one, the one you want + for i = 1:ncep + dctm(i,:) = cos((i-1)*[1:2:(2*nrow-1)]/(2*nrow)*pi) * sqrt(2/nrow); + end + if type == 2 + % make it unitary! (but not for HTK type 3) + dctm(1,:) = dctm(1,:)/sqrt(2); + end +elseif type == 4 % type 1 with implicit repeating of first, last bins + % Deep in the heart of the rasta/feacalc code, there is the logic + % that the first and last auditory bands extend beyond the edge of + % the actual spectra, and they are thus copied from their neighbors. + % Normally, we just ignore those bands and take the 19 in the middle, + % but when feacalc calculates mfccs, it actually takes the cepstrum + % over the spectrum *including* the repeated bins at each end. + % Here, we simulate 'repeating' the bins and an nrow+2-length + % spectrum by adding in extra DCT weight to the first and last + % bins. + for i = 1:ncep + dctm(i,:) = cos((i-1)*[1:nrow]/(nrow+1)*pi) * 2; + % Add in edge points at ends (includes fixup scale) + dctm(i,1) = dctm(i,1) + 1; + dctm(i,nrow) = dctm(i,nrow) + ((-1)^(i-1)); + end + dctm = dctm / (2*(nrow+1)); +else % dpwe type 1 - same as old spec2cep that expanded & used fft + for i = 1:ncep + dctm(i,:) = cos((i-1)*[0:(nrow-1)]/(nrow-1)*pi) * 2 / (2*(nrow-1)); + end + % fixup 'non-repeated' points + dctm(:,[1 nrow]) = dctm(:, [1 nrow])/2; +end + +cep = dctm*log(spec); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sceneClassificationMetrics_eval.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,50 @@ +% from https://soundsoftware.ac.uk/projects/aasp-d-case-metrics + +function [confusionMat, AccFolds, Acc, Std] = sceneClassificationMetrics_eval(numfolds, foldsGTlist, foldstestlist) + +% Function takes as input a number of folds (numfolds), a text file +% containing a list of filenames with the ground truth per fold (foldsGTlist), +% and a text file containing a list of filenames with the classification +% output (foldstestlist) +% +% e.g. [confusionMat, AccFolds, Acc, Std] = sceneClassificationMetrics_eval(5, 'foldGTlist.txt', 'foldtestlist.txt'); + + +% Initialize +classList = {'bus','busystreet','office','openairmarket','park','quietstreet','restaurant','supermarket','tube','tubestation'}; +confusionMat = zeros(10,10); +AccFolds = zeros(1,numfolds); + + +% For each fold +fid1 = fopen(foldsGTlist,'r+'); +fid2 = fopen(foldstestlist,'r+'); +for i=1:numfolds + + % Load classification output and ground truth + tline1 = fgetl(fid1); + [fileIDGT,classIDGT] = loadClassificationOutput(tline1); + tline2 = fgetl(fid2); + [fileID,classID] = loadClassificationOutput(tline2); + + % Compute confusion matrix per fold + confusionMatTemp = zeros(10,10); + for j=1:length(classIDGT) + pos = strmatch(fileIDGT{j}, fileID, 'exact'); + posClassGT = strmatch(classIDGT{j}, classList, 'exact'); + posClass = strmatch(classID{pos}, classList, 'exact'); + confusionMatTemp(posClassGT,posClass) = confusionMatTemp(posClassGT,posClass) + 1; + end + + % Compute accuracy per fold + AccFolds(i) = sum(diag(confusionMatTemp))/sum(sum(confusionMatTemp)); + confusionMat = confusionMat + confusionMatTemp; + +end +fclose(fid1); +fclose(fid2); + + +% Compute global accuracy and std +Acc = mean(AccFolds); +Std = std(AccFolds);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test.m Fri Nov 01 08:48:19 2013 +0000 @@ -0,0 +1,7 @@ +make_list_files; +classify_scenes('/tmp', 'fold1_train.txt','fold1_test.txt', 'fold1_result.txt', 1) +classify_scenes('/tmp', 'fold2_train.txt','fold2_test.txt', 'fold2_result.txt', 1) +classify_scenes('/tmp', 'fold3_train.txt','fold3_test.txt', 'fold3_result.txt', 1) +classify_scenes('/tmp', 'fold4_train.txt','fold4_test.txt', 'fold4_result.txt', 1) +classify_scenes('/tmp', 'fold5_train.txt','fold5_test.txt', 'fold5_result.txt', 1) +[confusionMat, AccFolds, Acc, Std] = sceneClassificationMetrics_eval(5, 'foldGTlist.txt', 'foldRESlist.txt')