initial

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