aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/cosine.cpp131
-rw-r--r--src/haar.cpp418
-rw-r--r--src/opcodes.cpp220
-rw-r--r--src/walsh.cpp1081
-rw-r--r--src/walsh2.cpp58
-rw-r--r--src/windowing.cpp208
6 files changed, 2116 insertions, 0 deletions
diff --git a/src/cosine.cpp b/src/cosine.cpp
new file mode 100644
index 0000000..78cafd8
--- /dev/null
+++ b/src/cosine.cpp
@@ -0,0 +1,131 @@
+/*
+ cosine.cpp
+ Copyright (C) 2015 John Burkardt
+ Copyright (C) 2019 Richard Knight
+
+
+ This program is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 3 of the License, or (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with this program; if not, write to the Free Software Foundation,
+ Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+
+ */
+
+#include <cmath>
+#include <cstdlib>
+#include <iostream>
+#include <ctime>
+#include <plugin.h>
+
+using namespace std;
+
+#include "cosine.hpp"
+
+//****************************************************************************80
+
+void cosine_transform_data (csnd::Vector<MYFLT> d, csnd::Vector<MYFLT> c)
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// COSINE_TRANSFORM_DATA does a cosine transform on a vector of data.
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 27 August 2015
+//
+// Author:
+//
+// John Burkardt
+//
+// Parameters:
+//
+// Input, integer N, the number of data points.
+//
+// Input, double D[N], the vector of data.
+//
+// Output, double COSINE_TRANSFORM_DATA[N], the transform coefficients.
+//
+{
+ int n = d.len();
+ double angle;
+ int i;
+ int j;
+ const double r8_pi = 3.141592653589793;
+
+
+ for ( i = 0; i < n; i++ )
+ {
+ c[i] = 0.0;
+ for ( j = 0; j < n; j++ )
+ {
+ angle = r8_pi * ( double ) ( i * ( 2 * j + 1 ) ) / ( double ) ( 2 * n );
+ c[i] = c[i] + cos ( angle ) * d[j];
+ }
+ c[i] = c[i] * sqrt ( 2.0 / ( double ) ( n ) );
+ }
+}
+//****************************************************************************80
+
+void cosine_transform_inverse (csnd::Vector<MYFLT> c, csnd::Vector<MYFLT> d)
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// COSINE_TRANSFORM_INVERSE does an inverse cosine transform.
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 27 August 2015
+//
+// Author:
+//
+// John Burkardt
+//
+// Parameters:
+//
+// Input, integer N, the number of data points.
+//
+// Input, double C[N], the vector of transform coefficients
+//
+// Output, double COSINE_TRANSFORM_INVERSE[N], the original data.
+//
+{
+ int n = c.len();
+ double angle;
+ int i;
+ int j;
+ double r8_pi = 3.141592653589793;
+
+
+ for ( i = 0; i < n; i++ )
+ {
+ d[i] = c[0] / 2.0;
+ for ( j = 1; j < n; j++ )
+ {
+ angle = r8_pi * ( double ) ( ( 2 * i + 1 ) * j ) / ( double ) ( 2 * n );
+ d[i] = d[i] + cos ( angle ) * c[j];
+ }
+ d[i] = d[i] * sqrt ( 2.0 / ( double ) ( n ) );
+ }
+}
+//****************************************************************************80
diff --git a/src/haar.cpp b/src/haar.cpp
new file mode 100644
index 0000000..36d39da
--- /dev/null
+++ b/src/haar.cpp
@@ -0,0 +1,418 @@
+/*
+ haar.cpp
+ Copyright (C) 2014 John Burkardt
+ Copyright (C) 2019 Richard Knight
+
+
+ This program is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 3 of the License, or (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with this program; if not, write to the Free Software Foundation,
+ Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+
+ */
+
+#include <cstdlib>
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+#include <ctime>
+#include <plugin.h>
+
+using namespace std;
+
+#include "haar.hpp"
+
+//****************************************************************************80
+
+void haar_1d (csnd::Csound* csound, csnd::Vector<MYFLT> x)
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// HAAR_1D computes the Haar transform of a vector.
+//
+// Discussion:
+//
+// For the classical Haar transform, N should be a power of 2.
+// However, this is not required here.
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 06 March 2014
+//
+// Author:
+//
+// John Burkardt
+//
+// Parameters:
+//
+// Input, int N, the dimension of the vector.
+//
+// Input/output, double X[N], on input, the vector to be transformed.
+// On output, the transformed vector.
+//
+{
+ int n = x.len();
+ int i;
+ int k;
+ double s;
+ double *y;
+
+ s = sqrt ( 2.0 );
+
+ //y = new double[n];
+ y = (double*) csound->malloc(sizeof(double) * n);
+//
+// Initialize.
+//
+ for ( i = 0; i < n; i++ )
+ {
+ y[i] = 0.0;
+ }
+//
+// Determine K, the largest power of 2 such that K <= N.
+//
+ k = 1;
+ while ( k * 2 <= n )
+ {
+ k = k * 2;
+ }
+
+ while ( 1 < k )
+ {
+ k = k / 2;
+ for ( i = 0; i < k; i++ )
+ {
+ y[i] = ( x[2*i] + x[2*i+1] ) / s;
+ y[i+k] = ( x[2*i] - x[2*i+1] ) / s;
+ }
+ for ( i = 0; i < k * 2; i++ )
+ {
+ x[i] = y[i];
+ }
+ }
+
+ //delete [] y;
+ csound->free(y);
+
+ return;
+}
+//****************************************************************************80
+
+void haar_1d_inverse (csnd::Csound* csound, csnd::Vector<MYFLT> x)
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// HAAR_1D_INVERSE computes the inverse Haar transform of a vector.
+//
+// Discussion:
+//
+// For the classical Haar transform, N should be a power of 2.
+// However, this is not required here.
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 06 March 2014
+//
+// Author:
+//
+// John Burkardt
+//
+// Parameters:
+//
+// Input, int N, the dimension of the vector.
+//
+// Input/output, double X[N], on input, the vector to be transformed.
+// On output, the transformed vector.
+//
+{
+
+ int n = x.len();
+ int i;
+ int k;
+ double s;
+ double *y;
+
+ s = sqrt ( 2.0 );
+
+ //y = new double[n];
+ y = (double*) csound->malloc(sizeof(double) * n);
+//
+// Initialize.
+//
+ for ( i = 0; i < n; i++ )
+ {
+ y[i] = 0.0;
+ }
+
+ k = 1;
+ while ( k * 2 <= n )
+ {
+ for ( i = 0; i < k; i++ )
+ {
+ y[2*i] = ( x[i] + x[i+k] ) / s;
+ y[2*i+1] = ( x[i] - x[i+k] ) / s;
+ }
+ for ( i = 0; i < k * 2; i++ )
+ {
+ x[i] = y[i];
+ }
+ k = k * 2;
+ }
+
+ //delete [] y;
+ csound->free(y);
+ return;
+}
+
+
+
+/*
+
+
+//****************************************************************************80
+
+void haar_2d ( int m, int n, double u[] )
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// HAAR_2D computes the Haar transform of an array.
+//
+// Discussion:
+//
+// For the classical Haar transform, M and N should be a power of 2.
+// However, this is not required here.
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 06 March 2014
+//
+// Author:
+//
+// John Burkardt
+//
+// Parameters:
+//
+// Input, int M, N, the dimensions of the array.
+//
+// Input/output, double U[M*N], the array to be transformed.
+//
+{
+ int i;
+ int j;
+ int k;
+ double s;
+ double *v;
+
+ s = sqrt ( 2.0 );
+
+ v = new double[m*n];
+
+ for ( j = 0; j < n; j++ )
+ {
+ for ( i = 0; i < m; i++ )
+ {
+ v[i+j*m] = u[i+j*m];
+ }
+ }
+//
+// Determine K, the largest power of 2 such that K <= M.
+//
+ k = 1;
+ while ( k * 2 <= m )
+ {
+ k = k * 2;
+ }
+//
+// Transform all columns.
+//
+ while ( 1 < k )
+ {
+ k = k / 2;
+
+ for ( j = 0; j < n; j++ )
+ {
+ for ( i = 0; i < k; i++ )
+ {
+ v[i +j*m] = ( u[2*i+j*m] + u[2*i+1+j*m] ) / s;
+ v[k+i+j*m] = ( u[2*i+j*m] - u[2*i+1+j*m] ) / s;
+ }
+ }
+ for ( j = 0; j < n; j++ )
+ {
+ for ( i = 0; i < 2 * k; i++ )
+ {
+ u[i+j*m] = v[i+j*m];
+ }
+ }
+ }
+//
+// Determine K, the largest power of 2 such that K <= N.
+//
+ k = 1;
+ while ( k * 2 <= n )
+ {
+ k = k * 2;
+ }
+//
+// Transform all rows.
+//
+ while ( 1 < k )
+ {
+ k = k / 2;
+
+ for ( j = 0; j < k; j++ )
+ {
+ for ( i = 0; i < m; i++ )
+ {
+ v[i+( j)*m] = ( u[i+2*j*m] + u[i+(2*j+1)*m] ) / s;
+ v[i+(k+j)*m] = ( u[i+2*j*m] - u[i+(2*j+1)*m] ) / s;
+ }
+ }
+
+ for ( j = 0; j < 2 * k; j++ )
+ {
+ for ( i = 0; i < m; i++ )
+ {
+ u[i+j*m] = v[i+j*m];
+ }
+ }
+ }
+ delete [] v;
+
+ return;
+}
+//****************************************************************************80
+
+void haar_2d_inverse ( int m, int n, double u[] )
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// HAAR_2D_INVERSE inverts the Haar transform of an array.
+//
+// Discussion:
+//
+// For the classical Haar transform, M and N should be a power of 2.
+// However, this is not required here.
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 06 March 2014
+//
+// Author:
+//
+// John Burkardt
+//
+// Parameters:
+//
+// Input, int M, N, the dimensions of the array.
+//
+// Input/output, double U[M*N], the array to be transformed.
+//
+{
+ int i;
+ int j;
+ int k;
+ double s;
+ double *v;
+
+ s = sqrt ( 2.0 );
+
+ v = new double[m*n];
+
+ for ( j = 0; j < n; j++ )
+ {
+ for ( i = 0; i < m; i++ )
+ {
+ v[i+j*m] = u[i+j*m];
+ }
+ }
+//
+// Inverse transform of all rows.
+//
+ k = 1;
+
+ while ( k * 2 <= n )
+ {
+ for ( j = 0; j < k; j++ )
+ {
+ for ( i = 0; i < m; i++ )
+ {
+ v[i+(2*j )*m] = ( u[i+j*m] + u[i+(k+j)*m] ) / s;
+ v[i+(2*j+1)*m] = ( u[i+j*m] - u[i+(k+j)*m] ) / s;
+ }
+ }
+
+ for ( j = 0; j < 2 * k; j++ )
+ {
+ for ( i = 0; i < m; i++ )
+ {
+ u[i+j*m] = v[i+j*m];
+ }
+ }
+ k = k * 2;
+ }
+//
+// Inverse transform of all columns.
+//
+ k = 1;
+
+ while ( k * 2 <= m )
+ {
+ for ( j = 0; j < n; j++ )
+ {
+ for ( i = 0; i < k; i++ )
+ {
+ v[2*i +j*m] = ( u[i+j*m] + u[k+i+j*m] ) / s;
+ v[2*i+1+j*m] = ( u[i+j*m] - u[k+i+j*m] ) / s;
+ }
+ }
+
+ for ( j = 0; j < n; j++ )
+ {
+ for ( i = 0; i < 2 * k; i++ )
+ {
+ u[i+j*m] = v[i+j*m];
+ }
+ }
+ k = k * 2;
+ }
+ delete [] v;
+
+ return;
+}
+
+ *
+ */ \ No newline at end of file
diff --git a/src/opcodes.cpp b/src/opcodes.cpp
new file mode 100644
index 0000000..ef4f2b1
--- /dev/null
+++ b/src/opcodes.cpp
@@ -0,0 +1,220 @@
+/*
+ opcodes.cpp
+ Copyright (C) 2019 Richard Knight
+
+
+ This program is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 3 of the License, or (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with this program; if not, write to the Free Software Foundation,
+ Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+
+ */
+
+#include <plugin.h>
+#include "cosine.hpp"
+#include "haar.hpp"
+#include "walsh.hpp"
+#include "walsh2.hpp"
+
+
+struct tfcosine : csnd::Plugin<1, 1> {
+ static constexpr char const *otypes = "k[]";
+ static constexpr char const *itypes = "k[]";
+
+ int init() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ out.init(csound, in.len());
+ return OK;
+ }
+
+ int kperf() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ cosine_transform_data(in, out);
+ return OK;
+ }
+};
+
+struct tfcosineinv : csnd::Plugin<1, 1> {
+ static constexpr char const *otypes = "k[]";
+ static constexpr char const *itypes = "k[]";
+
+ int init() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ out.init(csound, in.len());
+ return OK;
+ }
+
+ int kperf() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ cosine_transform_inverse(in, out);
+ return OK;
+ }
+};
+
+struct tfhaar1 : csnd::Plugin<1, 1> {
+ static constexpr char const *otypes = "k[]";
+ static constexpr char const *itypes = "k[]";
+
+ int init() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ out.init(csound, in.len());
+ return OK;
+ }
+
+ int kperf() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ std::copy(in.begin(), in.end(), out.begin());
+ haar_1d(csound, out);
+ return OK;
+ }
+};
+
+struct tfhaar1inv : csnd::Plugin<1, 1> {
+ static constexpr char const *otypes = "k[]";
+ static constexpr char const *itypes = "k[]";
+
+ int init() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ out.init(csound, in.len());
+ return OK;
+ }
+
+ int kperf() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ std::copy(in.begin(), in.end(), out.begin());
+ haar_1d_inverse(csound, out);
+ return OK;
+ }
+};
+
+struct tfhaar2 : csnd::Plugin<1, 1> {
+ static constexpr char const *otypes = "k[]";
+ static constexpr char const *itypes = "k[]";
+
+ int init() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ out.init(csound, in.len());
+ return OK;
+ }
+
+ int kperf() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ std::copy(in.begin(), in.end(), out.begin());
+ haar(csound, out);
+ return OK;
+ }
+};
+
+struct tfhaar2inv : csnd::Plugin<1, 1> {
+ static constexpr char const *otypes = "k[]";
+ static constexpr char const *itypes = "k[]";
+
+ int init() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ out.init(csound, in.len());
+ return OK;
+ }
+
+ int kperf() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ std::copy(in.begin(), in.end(), out.begin());
+ haarin(csound, out);
+ return OK;
+ }
+};
+
+struct tfwalsh1 : csnd::Plugin<1, 1> {
+ static constexpr char const *otypes = "k[]";
+ static constexpr char const *itypes = "k[]";
+
+ int init() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ out.init(csound, in.len());
+ return OK;
+ }
+
+ int kperf() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ std::copy(in.begin(), in.end(), out.begin());
+ walsh(csound, out);
+ return OK;
+ }
+};
+
+struct tfwalsh2 : csnd::Plugin<1, 1> {
+ static constexpr char const *otypes = "k[]";
+ static constexpr char const *itypes = "k[]";
+
+ int init() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ out.init(csound, in.len());
+ return OK;
+ }
+
+ int kperf() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ std::copy(in.begin(), in.end(), out.begin());
+ fwtanal(out);
+ return OK;
+ }
+};
+
+struct tfwalsh2inv : csnd::Plugin<1, 1> {
+ static constexpr char const *otypes = "k[]";
+ static constexpr char const *itypes = "k[]";
+
+ int init() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ out.init(csound, in.len());
+ return OK;
+ }
+
+ int kperf() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ std::copy(in.begin(), in.end(), out.begin());
+ fwtsynth(out);
+ return OK;
+ }
+};
+
+
+#include <modload.h>
+
+void csnd::on_load(csnd::Csound *csound) {
+ csnd::plugin<tfcosine>(csound, "tfcosine", csnd::thread::ik);
+ csnd::plugin<tfcosineinv>(csound, "tfcosineinv", csnd::thread::ik);
+ csnd::plugin<tfhaar1>(csound, "tfhaar1", csnd::thread::ik);
+ csnd::plugin<tfhaar1inv>(csound, "tfhaar1inv", csnd::thread::ik);
+ csnd::plugin<tfhaar2>(csound, "tfhaar2", csnd::thread::ik);
+ csnd::plugin<tfhaar2inv>(csound, "tfhaar2inv", csnd::thread::ik);
+ csnd::plugin<tfwalsh1>(csound, "tfwalsh1", csnd::thread::ik);
+ csnd::plugin<tfwalsh2>(csound, "tfwalsh2", csnd::thread::ik);
+ csnd::plugin<tfwalsh2inv>(csound, "tfwalsh2inv", csnd::thread::ik);
+}
diff --git a/src/walsh.cpp b/src/walsh.cpp
new file mode 100644
index 0000000..0ae574a
--- /dev/null
+++ b/src/walsh.cpp
@@ -0,0 +1,1081 @@
+/*
+ walsh.cpp
+ Copyright (C) 2004 John Burkardt
+ Copyright (C) 2011 Ken Beauchamp
+ Copyright (C) 2019 Richard Knight
+
+
+ This program is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 3 of the License, or (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with this program; if not, write to the Free Software Foundation,
+ Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+
+ */
+
+#include <cstdlib>
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+#include <ctime>
+#include <cstring>
+#include <plugin.h>
+
+using namespace std;
+
+#include "walsh.hpp"
+
+//****************************************************************************80
+
+void ffwt (csnd::Csound* csound, csnd::Vector<MYFLT> x)
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// FFWT performs an in-place fast Walsh transform.
+//
+// Discussion:
+//
+// This routine performs a fast Walsh transform on an input series X
+// leaving the transformed results in X.
+// X is dimensioned N, which must be a power of 2.
+// The results of this Walsh transform are in sequency order.
+//
+// The output sequence could be normalized by dividing by N.
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 16 March 2011
+//
+// Author:
+//
+// Ken Beauchamp
+//
+// Reference:
+//
+// Ken Beauchamp,
+// Walsh functions and their applications,
+// Academic Press, 1975,
+// ISBN: 0-12-084050-2,
+// LC: QA404.5.B33.
+//
+// Parameters:
+//
+// Input, int N, the number of items in X.
+// N must be a power of 2.
+//
+// Input/output, double X[N], the data to be transformed.
+//
+{
+ double hold;
+ int i;
+ int ii;
+ int j;
+ int j2;
+ int js;
+ int k;
+ int l;
+ int m;
+ int mw;
+ int mw1;
+ int nw;
+ int nz;
+ int nz2;
+ int nzi;
+ int nzn;
+ int *two_power;
+ double z;
+ int n = x.len();
+
+ m = i4_log_2 ( n );
+
+ two_power = (int*) csound->malloc(sizeof(int) * m); // new int[m];
+
+ for ( i = 0; i < m; i++ )
+ {
+ two_power[i] = i4_power ( 2, m - 1 - i );
+ }
+
+ for ( l = 0; l < m; l++ )
+ {
+ nz = i4_power ( 2, l );
+ nzi = 2 * nz;
+ nzn = n / nzi;
+ nz2 = nz / 2;
+ if ( nz2 == 0 )
+ {
+ nz2 = 1;
+ }
+
+ for ( i = 0; i < nzn; i++ )
+ {
+ js = i * nzi;
+ z = 1.0;
+ for ( ii = 0; ii < 2; ii++ )
+ {
+ for ( j = 0; j < nz2; j++ )
+ {
+ js = js + 1;
+ j2 = js + nz;
+ hold = x[js-1] + z * x[j2-1];
+ z = - z;
+ x[j2-1] = x[js-1] + z * x[j2-1];
+ x[js-1] = hold;
+ z = - z;
+ }
+ if ( l == 0 )
+ {
+ break;
+ }
+ z = - 1.0;
+ }
+ }
+ }
+//
+// Bit reversal section.
+//
+ nw = 0;
+ for ( k = 0; k < n; k++ )
+ {
+//
+// Choose correct index and switch elements if not already switched.
+//
+ if ( k < nw )
+ {
+ hold = x[nw];
+ x[nw] = x[k];
+ x[k] = hold;
+ }
+//
+// Bump up series by 1.
+//
+ for ( i = 0; i < m; i++ )
+ {
+ ii = i;
+ if ( nw < two_power[i] )
+ {
+ break;
+ }
+ mw = nw / two_power[i];
+ mw1 = mw / 2;
+ if ( mw <= 2 * mw1 )
+ {
+ break;
+ }
+ nw = nw - two_power[i];
+ }
+ nw = nw + two_power[ii];
+ }
+
+ //delete [] two_power;
+ csound->free(two_power);
+
+ return;
+}
+//****************************************************************************80
+
+void fwt (csnd::Csound* csound, csnd::Vector<MYFLT> x)
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// FWT performs a fast Walsh transform.
+//
+// Discussion:
+//
+// This routine performs a fast Walsh transform on an input series X
+// leaving the transformed results in X.
+// X is dimensioned N, which must be a power of 2.
+// The results of this Walsh transform are in sequency order.
+//
+// The output sequence could be normalized by dividing by N.
+//
+// Note that the program text in the reference included the line
+// y(jd) = abs ( x(j) - x(j2) )
+// which has been corrected to:
+// y(jd) = x(j) - x(j2)
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 16 March 2011
+//
+// Author:
+//
+// Ken Beauchamp
+//
+// Reference:
+//
+// Ken Beauchamp,
+// Walsh functions and their applications,
+// Academic Press, 1975,
+// ISBN: 0-12-084050-2,
+// LC: QA404.5.B33.
+//
+// Parameters:
+//
+// Input, int N, the number of items in X.
+// N must be a power of 2.
+//
+// Input/output, double X[N], the data to be transformed.
+//
+{
+ int i;
+ int j;
+ int j2;
+ int jd;
+ int js;
+ int l;
+ int m;
+ int n2;
+ int nx;
+ int ny;
+ int nz;
+ int nzi;
+ int nzn;
+ double *y;
+ int n = x.len();
+ y = (double*) csound->malloc(sizeof(double) * n);//new double[n];
+
+ n2 = n / 2;
+ m = i4_log_2 ( n );
+
+ for ( l = 1; l <= m; l++ )
+ {
+ ny = 0;
+ nz = i4_power ( 2, l - 1 );
+ nzi = 2 * nz;
+ nzn = n / nzi;
+ for ( i = 1; i <= nzn; i++ )
+ {
+ nx = ny + 1;
+ ny = ny + nz;
+ js = ( i - 1 ) * nzi;
+ jd = js + nzi + 1;
+ for ( j = nx; j <= ny; j++ )
+ {
+ js = js + 1;
+ j2 = j + n2;
+ y[js-1] = x[j-1] + x[j2-1];
+ jd = jd - 1;
+ y[jd-1] = x[j-1] - x[j2-1];
+ }
+ }
+ r8vec_copy (n, y, x );
+ }
+ csound->free(y); //delete [] y;
+
+ return;
+}
+//****************************************************************************80
+
+void haar (csnd::Csound* csound, csnd::Vector<MYFLT> x)
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// HAAR performs a Haar transform.
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 16 March 2011
+//
+// Author:
+//
+// Ken Beauchamp
+//
+// Reference:
+//
+// Ken Beauchamp,
+// Walsh functions and their applications,
+// Academic Press, 1975,
+// ISBN: 0-12-084050-2,
+// LC: QA404.5.B33.
+//
+// Parameters:
+//
+// Input, int N, the number of items in X.
+// N must be a power of 2.
+//
+// Input/output, double X[N], the data to be transformed.
+//
+{
+ int i;
+ int i1;
+ int j;
+ int jj;
+ int k;
+ int l;
+ int l2;
+ int l3;
+ double *y;
+ int n = x.len();
+
+ y = (double*) csound->malloc(sizeof(double) * n); //new double[n];
+
+ k = i4_log_2 ( n );
+
+ for ( i = 1; i <= k; i++ )
+ {
+ l = k + 1 - i;
+ l2 = i4_power ( 2, l - 1 );
+
+ r8vec_copy (2 * l2, x, y );
+
+ for ( j = 1; j <= l2; j++ )
+ {
+ l3 = l2 + j;
+ jj = 2 * j - 1;
+ x[j-1] = y[jj-1] + y[jj];
+ x[l3-1] = y[jj-1] - y[jj];
+ }
+ }
+ csound->free(y); //delete [] y;
+
+ return;
+}
+//****************************************************************************80
+
+void haarin (csnd::Csound* csound, csnd::Vector<MYFLT> x)
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// HAARIN inverts a Haar transform.
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 16 March 2011
+//
+// Author:
+//
+// Ken Beauchamp
+//
+// Reference:
+//
+// Ken Beauchamp,
+// Walsh functions and their applications,
+// Academic Press, 1975,
+// ISBN: 0-12-084050-2,
+// LC: QA404.5.B33.
+//
+// Parameters:
+//
+// Input, int N, the number of items in X.
+// N must be a power of 2.
+//
+// Input/output, double X[N], the data to be transformed.
+//
+{
+ int i;
+ int i1;
+ int j;
+ int jj;
+ int jj1;
+ int k;
+ int l;
+ int lj;
+ double *y;
+ int n = x.len();
+
+ y = (double*) csound->malloc(sizeof(double) * n); //new double[n];
+
+ k = i4_log_2 ( n );
+
+ for ( i = 1; i <= k; i++ )
+ {
+ l = i4_power ( 2, i - 1 );
+ r8vec_copy (2 * l, x, y );
+ for ( j = 1; j <= l; j++ )
+ {
+ lj = l + j;
+ jj = 2 * j;
+ jj1 = jj - 1;
+ x[jj-1] = y[j-1] - y[lj-1];
+ x[jj1-1] = y[j-1] + y[lj-1];
+ }
+ }
+
+ csound->free(y); //delete [] y;
+
+ return;
+}
+//****************************************************************************80
+
+void hnorm (csnd::Vector<MYFLT> x)
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// HNORM computes normalization factors for a forward or inverse Haar transform.
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 16 March 2011
+//
+// Author:
+//
+// Ken Beauchamp
+//
+// Reference:
+//
+// Ken Beauchamp,
+// Walsh functions and their applications,
+// Academic Press, 1975,
+// ISBN: 0-12-084050-2,
+// LC: QA404.5.B33.
+//
+// Parameters:
+//
+// Input, int N, the number of items in X.
+// N must be a power of 2.
+//
+// Input/output, double X[N], the data to be transformed.
+//
+{
+ int i;
+ int ii;
+ int j;
+ int jmax;
+ int jmin;
+ int k;
+ double wlk;
+ int n = x.len();
+
+ k = i4_log_2 ( n );
+
+ x[0] = x[0] / pow ( 2.0, k );
+
+ if ( 1 <= k )
+ {
+ x[1] = x[1] / pow ( 2.0, k );
+ }
+
+ for ( ii = 2; ii <= k; ii++ )
+ {
+ i = ii - 1;
+ wlk = 1.0 / pow ( 2.0, k - i );
+ jmin = i4_power ( 2, i );
+ jmax = i4_power ( 2, ii ) - 1;
+ for ( j = jmin; j <= jmax; j++ )
+ {
+ x[j] = x[j] * wlk;
+ }
+ }
+ return;
+}
+//****************************************************************************80
+
+int i4_log_2 ( int i )
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// I4_LOG_2 returns the integer part of the logarithm base 2 of an I4.
+//
+// Example:
+//
+// I I4_LOG_10
+// ----- --------
+// 0 0
+// 1 0
+// 2 1
+// 3 1
+// 4 2
+// 5 2
+// 7 2
+// 8 3
+// 9 3
+// 1000 9
+// 1024 10
+//
+// Discussion:
+//
+// I4_LOG_2 ( I ) + 1 is the number of binary digits in I.
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 04 January 2004
+//
+// Author:
+//
+// John Burkardt
+//
+// Parameters:
+//
+// Input, int I, the number whose logarithm base 2 is desired.
+//
+// Output, int I4_LOG_2, the integer part of the logarithm base 2 of
+// the absolute value of X.
+//
+{
+ int i_abs;
+ int two_pow;
+ int value;
+
+ if ( i == 0 )
+ {
+ value = 0;
+ }
+ else
+ {
+ value = 0;
+ two_pow = 2;
+
+ i_abs = abs ( i );
+
+ while ( two_pow <= i_abs )
+ {
+ value = value + 1;
+ two_pow = two_pow * 2;
+ }
+ }
+
+ return value;
+}
+//****************************************************************************80
+
+int i4_max ( int i1, int i2 )
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// I4_MAX returns the maximum of two I4's.
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 13 October 1998
+//
+// Author:
+//
+// John Burkardt
+//
+// Parameters:
+//
+// Input, int I1, I2, are two integers to be compared.
+//
+// Output, int I4_MAX, the larger of I1 and I2.
+//
+{
+ int value;
+
+ if ( i2 < i1 )
+ {
+ value = i1;
+ }
+ else
+ {
+ value = i2;
+ }
+ return value;
+}
+//****************************************************************************80
+
+int i4_min ( int i1, int i2 )
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// I4_MIN returns the minimum of two I4's.
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 13 October 1998
+//
+// Author:
+//
+// John Burkardt
+//
+// Parameters:
+//
+// Input, int I1, I2, two integers to be compared.
+//
+// Output, int I4_MIN, the smaller of I1 and I2.
+//
+{
+ int value;
+
+ if ( i1 < i2 )
+ {
+ value = i1;
+ }
+ else
+ {
+ value = i2;
+ }
+ return value;
+}
+//****************************************************************************80
+
+int i4_modp ( int i, int j )
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// I4_MODP returns the nonnegative remainder of I4 division.
+//
+// Discussion:
+//
+// If
+// NREM = I4_MODP ( I, J )
+// NMULT = ( I - NREM ) / J
+// then
+// I = J * NMULT + NREM
+// where NREM is always nonnegative.
+//
+// The MOD function computes a result with the same sign as the
+// quantity being divided. Thus, suppose you had an angle A,
+// and you wanted to ensure that it was between 0 and 360.
+// Then mod(A,360) would do, if A was positive, but if A
+// was negative, your result would be between -360 and 0.
+//
+// On the other hand, I4_MODP(A,360) is between 0 and 360, always.
+//
+// I J MOD I4_MODP I4_MODP Factorization
+//
+// 107 50 7 7 107 = 2 * 50 + 7
+// 107 -50 7 7 107 = -2 * -50 + 7
+// -107 50 -7 43 -107 = -3 * 50 + 43
+// -107 -50 -7 43 -107 = 3 * -50 + 43
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 26 May 1999
+//
+// Author:
+//
+// John Burkardt
+//
+// Parameters:
+//
+// Input, int I, the number to be divided.
+//
+// Input, int J, the number that divides I.
+//
+// Output, int I4_MODP, the nonnegative remainder when I is
+// divided by J.
+//
+{
+ int value;
+
+ if ( j == 0 )
+ {
+ cerr << "\n";
+ cerr << "I4_MODP - Fatal error!\n";
+ cerr << " I4_MODP ( I, J ) called with J = " << j << "\n";
+ exit ( 1 );
+ }
+
+ value = i % j;
+
+ if ( value < 0 )
+ {
+ value = value + abs ( j );
+ }
+
+ return value;
+}
+//****************************************************************************80
+
+int i4_power ( int i, int j )
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// I4_POWER returns the value of I^J.
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 01 April 2004
+//
+// Author:
+//
+// John Burkardt
+//
+// Parameters:
+//
+// Input, int I, J, the base and the power. J should be nonnegative.
+//
+// Output, int I4_POWER, the value of I^J.
+//
+{
+ int k;
+ int value;
+
+ if ( j < 0 )
+ {
+ if ( i == 1 )
+ {
+ value = 1;
+ }
+ else if ( i == 0 )
+ {
+ cerr << "\n";
+ cerr << "I4_POWER - Fatal error!\n";
+ cerr << " I^J requested, with I = 0 and J negative.\n";
+ exit ( 1 );
+ }
+ else
+ {
+ value = 0;
+ }
+ }
+ else if ( j == 0 )
+ {
+ if ( i == 0 )
+ {
+ cerr << "\n";
+ cerr << "I4_POWER - Fatal error!\n";
+ cerr << " I^J requested, with I = 0 and J = 0.\n";
+ exit ( 1 );
+ }
+ else
+ {
+ value = 1;
+ }
+ }
+ else if ( j == 1 )
+ {
+ value = i;
+ }
+ else
+ {
+ value = 1;
+ for ( k = 1; k <= j; k++ )
+ {
+ value = value * i;
+ }
+ }
+ return value;
+}
+//****************************************************************************80
+
+int i4_wrap ( int ival, int ilo, int ihi )
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// I4_WRAP forces an I4 to lie between given limits by wrapping.
+//
+// Example:
+//
+// ILO = 4, IHI = 8
+//
+// I Value
+//
+// -2 8
+// -1 4
+// 0 5
+// 1 6
+// 2 7
+// 3 8
+// 4 4
+// 5 5
+// 6 6
+// 7 7
+// 8 8
+// 9 4
+// 10 5
+// 11 6
+// 12 7
+// 13 8
+// 14 4
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 19 August 2003
+//
+// Author:
+//
+// John Burkardt
+//
+// Parameters:
+//
+// Input, int IVAL, an integer value.
+//
+// Input, int ILO, IHI, the desired bounds for the integer value.
+//
+// Output, int I4_WRAP, a "wrapped" version of IVAL.
+//
+{
+ int jhi;
+ int jlo;
+ int value;
+ int wide;
+
+ jlo = i4_min ( ilo, ihi );
+ jhi = i4_max ( ilo, ihi );
+
+ wide = jhi + 1 - jlo;
+
+ if ( wide == 1 )
+ {
+ value = jlo;
+ }
+ else
+ {
+ value = jlo + i4_modp ( ival - jlo, wide );
+ }
+
+ return value;
+}
+//****************************************************************************80
+
+void r8vec_copy (int n, double a1[], csnd::Vector<MYFLT> a2)
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// R8VEC_COPY copies an R8VEC.
+//
+// Discussion:
+//
+// An R8VEC is a vector of R8's.
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 03 July 2005
+//
+// Author:
+//
+// John Burkardt
+//
+// Parameters:
+//
+// Input, int N, the number of entries in the vectors.
+//
+// Input, double A1[N], the vector to be copied.
+//
+// Output, double A2[N], the copy of A1.
+//
+{
+ int i;
+
+ for ( i = 0; i < n; i++ )
+ {
+ a2[i] = a1[i];
+ }
+ return;
+}
+
+void r8vec_copy ( int n, csnd::Vector<MYFLT> a1, double a2[])
+{
+ int i;
+
+ for ( i = 0; i < n; i++ )
+ {
+ a2[i] = a1[i];
+ }
+ return;
+}
+
+//****************************************************************************80
+
+double *r8vec_copy_new (csnd::Csound* csound, int n, double a1[] )
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// R8VEC_COPY_NEW copies an R8VEC to a "new" R8VEC.
+//
+// Discussion:
+//
+// An R8VEC is a vector of R8's.
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 03 July 2008
+//
+// Author:
+//
+// John Burkardt
+//
+// Parameters:
+//
+// Input, int N, the number of entries in the vectors.
+//
+// Input, double A1[N], the vector to be copied.
+//
+// Output, double R8VEC_COPY_NEW[N], the copy of A1.
+//
+{
+ double *a2;
+ int i;
+
+ a2 = (double*) csound->malloc(sizeof(double) * n); //new double[n];
+
+ for ( i = 0; i < n; i++ )
+ {
+ a2[i] = a1[i];
+ }
+ return a2;
+}
+//****************************************************************************80
+
+//****************************************************************************80
+
+void walsh (csnd::Csound* csound, csnd::Vector<MYFLT> x)
+
+//****************************************************************************80
+//
+// Purpose:
+//
+// WALSH performs a fast Walsh transform.
+//
+// Discussion:
+//
+// This routine performs a fast Wash transform on an input series X
+// leaving the transformed results in X. The array Y is working space.
+// X and Y are dimensioned N, which must be a power of 2.
+// The results of this Walsh transform are in sequency order.
+//
+// The output sequence could be normalized by dividing by N.
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 16 March 2011
+//
+// Author:
+//
+// Ken Beauchamp
+//
+// Reference:
+//
+// Ken Beauchamp,
+// Walsh functions and their applications,
+// Academic Press, 1975,
+// ISBN: 0-12-084050-2,
+// LC: QA404.5.B33.
+//
+// Parameters:
+//
+// Input, int N, the number of items in X.
+// N must be a power of 2.
+//
+// Input/output, double X[N], the data to be transformed.
+//
+{
+ double a;
+ int i;
+ int i1;
+ int is;
+ int j;
+ int j1;
+ int l;
+ int m;
+ int n1;
+ int n2;
+ double w;
+ double *y;
+ double z;
+ int n = x.len();
+
+ n2 = n / 2;
+ y = (double*) csound->malloc(sizeof(double) * n2); //new double[n2];
+ m = i4_log_2 ( n );
+ z = - 1.0;
+
+ for ( j = 1; j <= m; j++ )
+ {
+ n1 = i4_power ( 2, m - j + 1 );
+ j1 = i4_power ( 2, j - 1 );
+ for ( l = 1; l <= j1; l++ )
+ {
+ is = ( l - 1 ) * n1 + 1;
+ i1 = 0;
+ w = z;
+ for ( i = is; i <= is + n1 - 1; i = i + 2 )
+ {
+ a = x[i-1];
+ x[is+i1-1] = a + x[i];
+ i1 = i1 + 1;
+ y[i1-1] = ( x[i] - a ) * w;
+ w = w * z;
+ }
+ for ( i = 1; i <= n1 / 2; i++ )
+ {
+ x[n1/2+is+i-2] = y[i-1];
+ }
+ }
+ }
+
+ csound->free(y);
+
+ return;
+}
diff --git a/src/walsh2.cpp b/src/walsh2.cpp
new file mode 100644
index 0000000..3e15592
--- /dev/null
+++ b/src/walsh2.cpp
@@ -0,0 +1,58 @@
+/*
+ walsh2.cpp
+ Copyright (C) 2018 Mottakin Chowdhury
+ Copyright (C) 2019 Richard Knight
+
+
+ This program is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 3 of the License, or (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with this program; if not, write to the Free Software Foundation,
+ Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+
+ */
+#include <plugin.h>
+
+const int N = 1<<16;
+
+
+void fwtanal(csnd::Vector<MYFLT> io) {
+ int n = io.len();
+ for (int d = 1; d < n; d <<= 1) {
+ for (int i = 0, m = d<<1; i < n; i += m) {
+ for (int j = 0; j < d; j++) { /// Don't forget modulo if required
+ MYFLT x = io[i+j], y = io[i+j+d];
+ io[i+j] = (x+y), io[i+j+d] = (x-y); // xor
+ // io[i+j] = x+y; // and
+ // io[i+j+d] = x+y; // or
+ }
+ }
+ }
+}
+
+
+void fwtsynth(csnd::Vector<MYFLT> io) {
+ int n = io.len();
+ for (int d = 1; d < n; d <<= 1) {
+ for (int i = 0, m = d<<1; i < n; i += m) {
+ for (int j = 0; j < d; j++) { /// Don't forget modulo if required
+ MYFLT x = io[i+j];
+ MYFLT y = io[i+j+d];
+ /// Modular inverse if required here
+ ///////io[i+j] = (x+y)>>1, io[i+j+d] = (x-y)>>1; // xor
+ io[i+j] = (x+y)/2;
+ io[i+j+d] = (x-y)/2;
+ // io[i+j] = x-y; // and
+ // io[i+j+d] = y-x; // or
+ }
+ }
+ }
+} \ No newline at end of file
diff --git a/src/windowing.cpp b/src/windowing.cpp
new file mode 100644
index 0000000..bdc7fdf
--- /dev/null
+++ b/src/windowing.cpp
@@ -0,0 +1,208 @@
+#include <plugin.h>
+
+typedef struct _fft {
+ OPDS h;
+ ARRAYDAT *out;
+ ARRAYDAT *in, *in2;
+ MYFLT *f;
+ MYFLT b;
+ int32_t n;
+ void *setup;
+ AUXCH mem;
+} FFT;
+
+
+static inline void tabinit(CSOUND *csound, ARRAYDAT *p, int size)
+{
+ size_t ss;
+ if (p->dimensions==0) {
+ p->dimensions = 1;
+ p->sizes = (int32_t*)csound->Malloc(csound, sizeof(int32_t));
+ }
+ if (p->data == NULL) {
+ CS_VARIABLE* var = p->arrayType->createVariable(csound, NULL);
+ p->arrayMemberSize = var->memBlockSize;
+ ss = p->arrayMemberSize*size;
+ p->data = (MYFLT*)csound->Calloc(csound, ss);
+ p->allocated = ss;
+ } else if( (ss = p->arrayMemberSize*size) > p->allocated) {
+ p->data = (MYFLT*) csound->ReAlloc(csound, p->data, ss);
+ p->allocated = ss;
+ }
+ if (p->dimensions==1) p->sizes[0] = size;
+ //p->dimensions = 1;
+}
+
+int32_t shiftin_init(CSOUND *csound, FFT *p) {
+ int32_t sizs = CS_KSMPS;
+ if(p->out->sizes[0] < sizs)
+ tabinit(csound, p->out, sizs);
+ p->n = 0;
+ return OK;
+}
+
+int32_t shiftin_perf(CSOUND *csound, FFT *p) {
+ IGN(csound);
+ uint32_t siz = p->out->sizes[0], n = p->n;
+ MYFLT *in = ((MYFLT *) p->in);
+ if (n + CS_KSMPS < siz) {
+ memcpy(p->out->data+n,in,CS_KSMPS*sizeof(MYFLT));
+ }
+ else {
+ int32_t num = siz - n;
+ memcpy(p->out->data+n,in,num*sizeof(MYFLT));
+ memcpy(p->out->data,in+num,(CS_KSMPS-num)*sizeof(MYFLT));
+ }
+ p->n = (n + CS_KSMPS)%siz;
+ return OK;
+}
+
+
+int32_t shiftout_init(CSOUND *csound, FFT *p) {
+ int32_t siz = p->in->sizes[0];
+ p->n = ((int32_t)*((MYFLT *)p->in2) % siz);
+ if (UNLIKELY((uint32_t) siz < CS_KSMPS))
+ return csound->InitError(csound, "%s", Str("input array too small\n"));
+ return OK;
+}
+
+int32_t shiftout_perf(CSOUND *csound, FFT *p) {
+ IGN(csound);
+ uint32_t siz = p->in->sizes[0], n = p->n;
+ MYFLT *out = ((MYFLT *) p->out);
+
+ if (n + CS_KSMPS < siz) {
+ memcpy(out,p->in->data+n,CS_KSMPS*sizeof(MYFLT));
+ }
+ else {
+ int32_t num = siz - n;
+ memcpy(out,p->in->data+n,num*sizeof(MYFLT));
+ memcpy(out+num,p->in->data,(CS_KSMPS-num)*sizeof(MYFLT));
+ }
+ p->n = (n + CS_KSMPS)%siz;
+ return OK;
+}
+
+
+struct windowing : csnd::Plugin<1, 1> {
+ static constexpr char const *otypes = "k[]";
+ static constexpr char const *itypes = "k[]";
+ FWT<MYFLT> fwt;
+
+ int init() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ out.init(csound, in.len());
+ return OK;
+ }
+
+ int kperf() {
+ csnd::Vector<MYFLT> &out = outargs.vector_data<MYFLT>(0);
+ csnd::Vector<MYFLT> &in = inargs.vector_data<MYFLT>(0);
+ std::copy(in.begin(), in.end(), out.begin());
+ fwt.fwtanal(out);
+ return OK;
+ }
+};
+
+int32_t init_window(CSOUND *csound, FFT *p) {
+ int32_t N = p->in->sizes[0];
+ int32_t i,type = (int32_t) *p->f;
+ MYFLT *w;
+ tabinit(csound, p->out, N);
+ if (p->mem.auxp == 0 || p->mem.size < N*sizeof(MYFLT))
+ csound->AuxAlloc(csound, N*sizeof(MYFLT), &p->mem);
+ w = (MYFLT *) p->mem.auxp;
+ switch(type) {
+ case 0:
+ for (i=0; i<N; i++) w[i] = 0.54 - 0.46*cos(i*TWOPI/N);
+ break;
+ case 1:
+ default:
+ for (i=0; i<N; i++) w[i] = 0.5 - 0.5*cos(i*TWOPI/N);
+ }
+ return OK;
+}
+
+int32_t perf_window(CSOUND *csound, FFT *p) {
+ IGN(csound);
+ int32_t i,end = p->out->sizes[0], off = *((MYFLT *)p->in2);
+ MYFLT *in, *out, *w;
+ in = p->in->data;
+ out = p->out->data;
+ w = (MYFLT *) p->mem.auxp;
+ /*while (off < 0) off += end;
+ for (i=0;i<end;i++)
+ out[(i+off)%end] = in[i]*w[i];*/
+ if(off) off = end - off;
+ for(i=0;i<end;i++)
+ out[i] = in[i]*w[(i+off)%end];
+ return OK;
+
+}
+
+
+
+static inline void tabensure2D(CSOUND *csound, ARRAYDAT *p,
+ int32_t rows, int32_t columns)
+{
+ if (p->data==NULL || p->dimensions == 0 ||
+ (p->dimensions==2 && (p->sizes[0] < rows || p->sizes[1] < columns))) {
+ size_t ss;
+ if (p->data == NULL) {
+ CS_VARIABLE* var = p->arrayType->createVariable(csound, NULL);
+ p->arrayMemberSize = var->memBlockSize;
+ }
+ ss = p->arrayMemberSize*rows*columns;
+ if (p->data==NULL) {
+ p->data = (MYFLT*)csound->Calloc(csound, ss);
+ p->dimensions = 2;
+ p->sizes = (int32_t*)csound->Malloc(csound, sizeof(int32_t)*2);
+ }
+ else p->data = (MYFLT*) csound->ReAlloc(csound, p->data, ss);
+ p->sizes[0] = rows; p->sizes[1] = columns;
+ }
+}
+// setrow
+int32_t set_rows_init(CSOUND *csound, FFT *p) {
+ int32_t sizs = p->in->sizes[0];
+ int32_t row = *((MYFLT *)p->in2);
+ tabensure2D(csound, p->out, row+1, sizs);
+ return OK;
+}
+
+int32_t set_rows_perf(CSOUND *csound, FFT *p) {
+ int32_t start = *((MYFLT *)p->in2);
+ if (UNLIKELY(start < 0 || start >= p->out->sizes[0]))
+ return csound->PerfError(csound, &(p->h),
+ Str("Error: index out of range\n"));
+ int32_t bytes = p->in->sizes[0]*sizeof(MYFLT);
+ start *= p->out->sizes[1];
+ memcpy(p->out->data+start,p->in->data,bytes);
+ return OK;
+}
+
+
+// getrow
+int32_t rows_init(CSOUND *csound, FFT *p) {
+ if (p->in->dimensions == 2) {
+ int32_t siz = p->in->sizes[1];
+ tabinit(csound, p->out, siz);
+ return OK;
+ }
+ else
+ return csound->InitError(csound, "%s",
+ Str("in array not 2-dimensional\n"));
+}
+
+int32_t rows_perf(CSOUND *csound, FFT *p) {
+ int32_t start = *((MYFLT *)p->in2);
+ if (LIKELY(start < p->in->sizes[0])) {
+ int32_t bytes = p->in->sizes[1]*sizeof(MYFLT);
+ start *= p->in->sizes[1];
+ memcpy(p->out->data,p->in->data+start,bytes);
+ return OK;
+ }
+ else return csound->PerfError(csound, &(p->h),
+ Str("requested row is out of range\n"));
+} \ No newline at end of file