From 1105cd57a0e820e36b230f533b7c42245ab508f4 Mon Sep 17 00:00:00 2001
From: Chikashi Miyama <chikashimiyama@gmail.com>
Date: Thu, 1 Apr 2021 21:52:56 +0200
Subject: [PATCH 1/4] add convolution

---
 .gitmodules               |  1 +
 convolution/Makefile      | 14 ++++++
 convolution/convolution.c | 17 +++++++
 convolution/convolution.h |  4 ++
 convolution/fft-dif.c     | 97 +++++++++++++++++++++++++++++++++++++++
 convolution/fft.h         | 20 ++++++++
 convolution/ift-dif.c     | 97 +++++++++++++++++++++++++++++++++++++++
 convolution/ift.h         | 20 ++++++++
 convolution/test.c        | 21 +++++++++
 9 files changed, 291 insertions(+)
 create mode 100644 .gitmodules
 create mode 100644 convolution/Makefile
 create mode 100644 convolution/convolution.c
 create mode 100644 convolution/convolution.h
 create mode 100644 convolution/fft-dif.c
 create mode 100644 convolution/fft.h
 create mode 100644 convolution/ift-dif.c
 create mode 100644 convolution/ift.h
 create mode 100644 convolution/test.c

diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000..8b13789
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1 @@
+
diff --git a/convolution/Makefile b/convolution/Makefile
new file mode 100644
index 0000000..fda2928
--- /dev/null
+++ b/convolution/Makefile
@@ -0,0 +1,14 @@
+PROGRAM = convolution_test
+OBJS = convolution.o test.o fft-dif.o ift-dif.o
+
+.SUFFIXES: .c .o
+
+$(PROGRAM): $(OBJS)
+	$(CC) -o $(PROGRAM) $^
+
+.c.o:
+	$(CC) $(CFLAGS) -c $<
+
+.PHONY: clean
+clean:
+	$(RM) $(PROGRAM) $(OBJS)
diff --git a/convolution/convolution.c b/convolution/convolution.c
new file mode 100644
index 0000000..2e06049
--- /dev/null
+++ b/convolution/convolution.c
@@ -0,0 +1,17 @@
+#include "convolution.h"
+#include "fft.h"
+#include "ift.h"
+
+void complexMultiplication(float complex* sample, float complex* ir, float complex* result, int size)
+{
+    for(int i = 0; i < size; ++i)
+        *result++ = *sample++ * *ir++;
+}
+
+void convolve(float complex* sample, float complex* ir, float complex* result, int size)
+{
+	fft(sample, size);
+	fft(ir, size);
+    complexMultiplication(sample, ir, result, size);
+    ift(result, size);
+}
diff --git a/convolution/convolution.h b/convolution/convolution.h
new file mode 100644
index 0000000..f4789f6
--- /dev/null
+++ b/convolution/convolution.h
@@ -0,0 +1,4 @@
+#include <complex.h>
+
+void convolve(float complex* sample, float complex* ir, float complex* result, int size);
+
diff --git a/convolution/fft-dif.c b/convolution/fft-dif.c
new file mode 100644
index 0000000..fc94866
--- /dev/null
+++ b/convolution/fft-dif.c
@@ -0,0 +1,97 @@
+#include "fft.h"
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+static int ctz(size_t N)
+{
+	int ctz1 = 0;
+
+	while( N ) {
+		ctz1++;
+		N >>= 1;
+	}
+
+	return ctz1-1;
+}
+
+static void nop_split(const float complex *x, float complex *X, size_t N)
+{
+	for(size_t n = 0; n < N/2; n++) {
+		X[2*n+0] = x[0/2+n];
+		X[2*n+1] = x[N/2+n];
+	}
+}
+
+static void fft_split(const float complex *x, float complex *X, size_t N, float complex phi)
+{
+	for(size_t n = 0; n < N/2; n++) {
+		X[2*n+0] = (x[0/2+n] + x[N/2+n]);
+		X[2*n+1] = (x[0/2+n] - x[N/2+n]) * cexpf(-2*(float)M_PI*I*phi);
+	}
+}
+
+static size_t revbits(size_t v, int J)
+{
+	size_t r = 0;
+
+	for(int j = 0; j < J; j++) {
+		r |= ( (v>>j)&1 ) << (J-1-j);
+	}
+
+	return r;
+}
+
+static int nop_reverse(int b, float complex *buffers[2], size_t N)
+{
+	int J = ctz(N);
+
+	for(int j = J-2; j >= 0; j--, b++) {
+		size_t delta = N>>j;
+
+		for(size_t n = 0; n < N; n += delta) {
+			nop_split(buffers[b&1]+n, buffers[~b&1]+n, delta);
+		}
+	}
+
+	return b;
+}
+
+static int fft_reverse(int b, float complex *buffers[2], size_t N)
+{
+	int J = ctz(N);
+
+	for(int j = J-1; j >= 0; j--, b++) {
+		size_t delta = N>>j;
+
+		for(size_t n = 0; n < N; n += delta) {
+			float complex phi = (float)revbits(n/delta, j) / (float)(2<<j);
+			fft_split(buffers[b&1]+n, buffers[~b&1]+n, delta, phi);
+		}
+	}
+
+	return b;
+}
+
+int fft(float complex *vector, size_t N)
+{
+	if( !N ) return 0;
+
+	if( N & (N-1) ) return 1;
+
+	float complex *buffers[2] = { vector, malloc(N*sizeof(float complex)) };
+
+	if( !buffers[1] ) return -1;
+
+	int b = 0;
+
+	b = nop_reverse(b, buffers, N);
+	b = fft_reverse(b, buffers, N);
+	b = nop_reverse(b, buffers, N);
+
+	memmove(vector, buffers[b&1], N*sizeof(float complex));
+
+	free( buffers[1] );
+
+	return 0;
+}
diff --git a/convolution/fft.h b/convolution/fft.h
new file mode 100644
index 0000000..02d92bd
--- /dev/null
+++ b/convolution/fft.h
@@ -0,0 +1,20 @@
+#ifndef FFT_H
+#define FFT_H
+
+#include <complex.h>
+#include <stddef.h>
+
+/**
+ * @brief FFT algorithm (forward transform)
+ *
+ * This function computes forward radix-2 fast Fourier transform (FFT).
+ * The output is written in-place over the input.
+ *
+ * @param vector An array of @p N complex values in single-precision floating-point format.
+ * @param N The size of the transform must be a power of two.
+ *
+ * @return Zero for success.
+ */
+int fft(float complex *vector, size_t N);
+
+#endif
diff --git a/convolution/ift-dif.c b/convolution/ift-dif.c
new file mode 100644
index 0000000..39d06db
--- /dev/null
+++ b/convolution/ift-dif.c
@@ -0,0 +1,97 @@
+#include "ift.h"
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+static int ctz(size_t N)
+{
+	int ctz1 = 0;
+
+	while( N ) {
+		ctz1++;
+		N >>= 1;
+	}
+
+	return ctz1-1;
+}
+
+static void nop_split(const float complex *x, float complex *X, size_t N)
+{
+	for(size_t n = 0; n < N/2; n++) {
+		X[2*n+0] = x[0/2+n];
+		X[2*n+1] = x[N/2+n];
+	}
+}
+
+static void fft_split(const float complex *x, float complex *X, size_t N, float complex phi)
+{
+	for(size_t n = 0; n < N/2; n++) {
+		X[2*n+0] = (x[0/2+n] + x[N/2+n])/2;
+		X[2*n+1] = (x[0/2+n] - x[N/2+n])/2 * cexpf(+2*(float)M_PI*I*phi);
+	}
+}
+
+static size_t revbits(size_t v, int J)
+{
+	size_t r = 0;
+
+	for(int j = 0; j < J; j++) {
+		r |= ( (v>>j)&1 ) << (J-1-j);
+	}
+
+	return r;
+}
+
+static int nop_reverse(int b, float complex *buffers[2], size_t N)
+{
+	int J = ctz(N);
+
+	for(int j = J-2; j >= 0; j--, b++) {
+		size_t delta = N>>j;
+
+		for(size_t n = 0; n < N; n += delta) {
+			nop_split(buffers[b&1]+n, buffers[~b&1]+n, delta);
+		}
+	}
+
+	return b;
+}
+
+static int fft_reverse(int b, float complex *buffers[2], size_t N)
+{
+	int J = ctz(N);
+
+	for(int j = J-1; j >= 0; j--, b++) {
+		size_t delta = N>>j;
+
+		for(size_t n = 0; n < N; n += delta) {
+			float complex phi = (float)revbits(n/delta, j) / (float)(2<<j);
+			fft_split(buffers[b&1]+n, buffers[~b&1]+n, delta, phi);
+		}
+	}
+
+	return b;
+}
+
+int ift(float complex *vector, size_t N)
+{
+	if( !N ) return 0;
+
+	if( N & (N-1) ) return 1;
+
+	float complex *buffers[2] = { vector, malloc(N*sizeof(float complex)) };
+
+	if( !buffers[1] ) return -1;
+
+	int b = 0;
+
+	b = nop_reverse(b, buffers, N);
+	b = fft_reverse(b, buffers, N);
+	b = nop_reverse(b, buffers, N);
+
+	memmove(vector, buffers[b&1], N*sizeof(float complex));
+
+	free( buffers[1] );
+
+	return 0;
+}
diff --git a/convolution/ift.h b/convolution/ift.h
new file mode 100644
index 0000000..d00415d
--- /dev/null
+++ b/convolution/ift.h
@@ -0,0 +1,20 @@
+#ifndef IFT_H
+#define IFT_H
+
+#include <complex.h>
+#include <stddef.h>
+
+/**
+ * @brief FFT algorithm (inverse transform)
+ *
+ * This function computes inverse radix-2 fast Fourier transform (FFT).
+ * The output is written in-place over the input.
+ *
+ * @param vector An array of @p N complex values in single-precision floating-point format.
+ * @param N The size of the transform must be a power of two.
+ *
+ * @return Zero for success.
+ */
+int ift(float complex *vector, size_t N);
+
+#endif
diff --git a/convolution/test.c b/convolution/test.c
new file mode 100644
index 0000000..1621bad
--- /dev/null
+++ b/convolution/test.c
@@ -0,0 +1,21 @@
+#include <stdio.h>
+#include "convolution.h"
+
+int main()
+{
+	float complex dirac1[128] = {0.0f};
+	float complex dirac2[128] = {0.0f};
+	float complex result[128] = {0.0f};
+
+	dirac1[2] = 0.8f;
+	dirac2[2] = 0.5f;
+
+	convolve(dirac1, dirac2, result, 128);
+
+	for(int i = 0; i < 10; i++)
+	{
+		printf("index:%d value:%f\n", i, creal(result[i]));
+	}
+
+	return 0;
+}
\ No newline at end of file

From 26a54caaeac7cc4e911362d7e331626ccb300cf4 Mon Sep 17 00:00:00 2001
From: Chikashi Miyama <chikashimiyama@gmail.com>
Date: Thu, 1 Apr 2021 22:17:22 +0200
Subject: [PATCH 2/4] rename dir

---
 Makefile                        |   3 ++-
 {convolution => conv}/Makefile  |   0
 conv/convolution.c              |  23 +++++++++++++++++++++++
 conv/convolution.h              |   5 +++++
 conv/convolution_test           | Bin 0 -> 50048 bytes
 {convolution => conv}/fft-dif.c |   0
 {convolution => conv}/fft.h     |   0
 {convolution => conv}/ift-dif.c |   0
 {convolution => conv}/ift.h     |   0
 conv/test.c                     |  30 ++++++++++++++++++++++++++++++
 convolution/convolution.c       |  17 -----------------
 convolution/convolution.h       |   4 ----
 convolution/test.c              |  21 ---------------------
 13 files changed, 60 insertions(+), 43 deletions(-)
 rename {convolution => conv}/Makefile (100%)
 create mode 100644 conv/convolution.c
 create mode 100644 conv/convolution.h
 create mode 100755 conv/convolution_test
 rename {convolution => conv}/fft-dif.c (100%)
 rename {convolution => conv}/fft.h (100%)
 rename {convolution => conv}/ift-dif.c (100%)
 rename {convolution => conv}/ift.h (100%)
 create mode 100644 conv/test.c
 delete mode 100644 convolution/convolution.c
 delete mode 100644 convolution/convolution.h
 delete mode 100644 convolution/test.c

diff --git a/Makefile b/Makefile
index 3ac4930..f232410 100644
--- a/Makefile
+++ b/Makefile
@@ -2,7 +2,8 @@
 
 lib.name = earplug~
 
-class.sources = earplug~.c
+cflags = -I./conv
+class.sources = earplug~.c conv/convolution.c conv/fft-dif.c conv/ift-dif.c
 
 datafiles = earplug~-help.pd earplug~-meta.pd earplug_data_src.txt \
             parse-to-h.pl README.txt LICENSE.txt
diff --git a/convolution/Makefile b/conv/Makefile
similarity index 100%
rename from convolution/Makefile
rename to conv/Makefile
diff --git a/conv/convolution.c b/conv/convolution.c
new file mode 100644
index 0000000..07707ca
--- /dev/null
+++ b/conv/convolution.c
@@ -0,0 +1,23 @@
+#include "convolution.h"
+#include "fft.h"
+#include "ift.h"
+
+void complexMultiplication(float complex* sample, float complex* ir, int size)
+{
+    while(size--)
+        *sample++ *= *ir++;
+}
+
+void convolve(float complex* sample, float complex* ir, int size)
+{
+	fft(sample, size);
+	fft(ir, size);
+    complexMultiplication(sample, ir, size);
+    ift(sample, size);
+}
+
+void clean(float complex* buffer, int size)
+{
+	while(size--)
+		*buffer++ = 0.f + 0.fI;
+}
\ No newline at end of file
diff --git a/conv/convolution.h b/conv/convolution.h
new file mode 100644
index 0000000..784d740
--- /dev/null
+++ b/conv/convolution.h
@@ -0,0 +1,5 @@
+#include <complex.h>
+
+void convolve(float complex* sample, float complex* ir, int size);
+
+void clean(float complex* buf, int size);
diff --git a/conv/convolution_test b/conv/convolution_test
new file mode 100755
index 0000000000000000000000000000000000000000..3c7547fbf6bd096d8078f5a8515d0598c3e3e86a
GIT binary patch
literal 50048
zcmeHP4{%gPntuby85QV>1}=Iab2{+`U6CLvabd~mOxdF|nAHTs?%)uTjL8|2ND_rB
zU89q^^*N4}qVUSv-MaGX?y9^y>*TJq4tK$uBrpjsgplx05Q%^|L$V|!KtjOG{l0#0
zGLr$7y}hcf+Ur;J>;As({`%|R>#tw`W+v<z``!3;iei|mD9RKhBhq{6isDyHC@IP<
zNF2%IDX>0j|B3y{CrCT_GE5e)Dv=!v2Jm<a?Zt(Yn~47TlXZ$N!xe?npXt>ePks6N
zdR;)4YOm`yOa}Vexmp{N(-O@|XiOG8p1S&mHFXoJRD1b9((EnJTg2sfhYa0b!auU-
zsVJ|iEnkO5srG7hd#iN=Bopov*A@Pi)>PLeCY)-oOSc!&4Ujx(z7vMEg2z)@SHG^Z
z`dM9{YOi^jHs71|dXkf`Nxq`O<7udVv9h|%Q(5h;5%T2c`%k*PV%-SIsWy`B(>i%a
z*s)m)t)54oSpJl@TJkl&;_D(cr}ayTo;|H$kLTH%1ZpH--5#|SSH9L?m&AI@hzrQc
zuLzQ6C#j+_#V^R>%Ga{EU-3D6JXIy0x(#dB)KqzD>(_m&J-xd^tCeiTFT|>HLtH@a
z@syR+mn0Oa{^d7n{>A4)a<aWdPpy=Qcc)|d<93@vLtjxsDf5I(SEFu>aHkZ%Zby}f
znUJE;*cr&<F+^<^LvBRAMb8PJY$V>IC^sSB+^Q&D$Ztd{K%y%?24O%plX0N><f{Vx
z$xiAu0}2z;G^BeQ>ek&`Rk`NgvJF*bkm>phB=Y0b8%6i3%OASYw`=jj_aeV|dLeWs
zB=SLF7#WGmV)_{+Q|?09NOk&U!Su+^W>Z|F<uvYMB&wr1!y>jfL&6vM^VEjA`tr5&
zAD$1tD%Xfc%ALr_{^AxbpQp=&m%Jnb5&?;TL_i`S5s(N-1SA3y0f~S_Kq4R!kO)Ww
zBmxoviGV~vA|Mfv2>dA{P{h@*xUau}t1pjowP8>#>JQla)fQh%Om+49MitgX-`ntj
zeUPgy{M8ovB8jW6!C_`pka!xZ&@-uO6lx<<aT>E6A7(3inWYa|lv#$5on@9w$k53B
z)i5g_Ve?NiOFy%W53@}B=eW5Jan;_1LeVbVR(0T3z8$xQ5N>O?;`Z_u7ND=yF@`<F
z=0k6EQ41GU+EGbgK2aqM7+o-6EP`R<szZ^V-KHpBcC)>gYU$=WNnNO)YYXY<=GuxW
zkpme>`DxMq@mP%ewim2wV{;W9FH;g;JOJoQ8~P*7_<{<;APK5H5^2=rQ66X*3=~0i
zG7=<P@%ABx_P*_=Rc+q+Y;KcSX2L$6h1zf~U_YlnhjZdN+(F1+0XBOY^UuU|$<>{)
zSZq~WbJ8X3z1p9T$Hs&H!4V7}3=Uq3#W2zI*c^e=^h?NdyLsUL+vd^`4`3mTw)^|h
zn0(@DCs#+*1J*4{nxd$0_d@SfTe<2#MCR0B7cz|R-=*cj4$;6LLKay7(>$QJY!7@y
zoq6DGQJ19+@2B?EU%4EMv8Q*U7)(PE9=@%cd|Me6O@e92d0-oj$O9kI_`hz1i;;5R
z+U}1)CXapFsB{AhE`bpEhy-3LYGJc!hc^A6gm`f2Hn_0a5NV!`cA|kVWKcyR`d`43
z2e%PqpESZwq=cSb@FNYwZE*~LgRa7Pc&0rIKRKzn)%gv%brL3E0B%uf!YvZu7G<e!
z-Kn{C5#O!`L<8X#wba}Sq}K~v0j+D4=XEeYeeata;;(L}C%vjI@(XljH#RpHS)aD9
zuohVhpW?6fkD?I|q%YR{s!irW()cg1rZf1Wl?Ru{w$PWf+~3LkUyw$<#MQ0lUeX1B
za2^>{0}k~<#0mQx1|p>Lhgj2o(Y(q0x@hyhxmo1DGQWe|!M26DZ*W=e1+Lo6t9_$O
zUFx7ytv4Il&x<kO$4k-ZZl`&5?uEX;O2h)2n-6#E%?n~-9&nnoxql}QI%4cCGy1sH
zHm*YH=;?E+hxy*q9QGElraqWb2YtN;HS8NS@W7krdZBR!-rP>PRqeL=jvJi84bzqf
z7afI0aRp0eb9F!WeKl^+{av4F^7t;pQI}!FWjN=`+rfO*v?hFGhK3Ay`Y8G5RKIqr
zmvS$}<Cg~>ImgvUPx64RhYI~rTt<G56jZp(ec>tG*Jk9tb`#G#->}=|@2-F6v0z&C
z@!*mu56m8PsDqA{-n7Tm+1GKy0Uk&jd`z`nm*EW9%*N<*a3SGhb8aXh2!9NRoOz*j
zqtQ08(-7K&%UIKE*8`OnTd)@K+GKtg)AohN@hq0)c6)8CW=_tn*cffhlJf(tP&4yB
zy|DBl6#5=a%ujzFo6K8aCDv+<b;g|LCET~w{2mm9xOp=Q?5!CXldF4pY^O_Yar8v(
z>ILy&c;J2UEVx=_&c`!yse2rI2Y8@K*V@b$U(Dd7hbx{bSGBRYs#&HLNkg%>p1@B5
zT>cnq`dNme*vL2#ZZaQ2<w*C`O=6x~=vM<h*z51YK53kug&%ATvuEh7;xwBa!AG%`
zb#mXSajt%74xuSm3(dv8<1n<z+y&Xos#}0eZKU91BUO)lNkep?$eL|t=2S28E!eE+
zdWHXG3$4itJj~r_09!UQ=YcA-#fe)!;uDGjd+tdXe5J}<i}#6B?al3WjGXnoWG+st
z%gpW0jYV(8@{Z%}b+YP^*N$zGfTfvv@wmL%Xz}z;6muR&5%rdHUN+Z;q7|;%P;m)q
zznAs&K(h}*ZQF8a&O<Xe#Ee5I??Aa76$eq#a~jPGFF?D`%LZnWeJIBF)5s@LS_8FZ
z+{=c?&^dxmt>`p7NS)3=GlBxy8{1C3TB*Z!%yC60dcld|F7yadj|du{W-E51#lB-+
zHlEY}=YYxAa+&!b27*?$0+Xt3LDRk#FK0h%r&-fu;OeLRhfZ^G{oU9Qmdra4$A_zJ
z>R#T{=R09g&+?uV{Q6@&)*9<{<Q-@@!TquIGiX!h?5_5{*<wFYNBQ1}?`s3UzCB)p
zV+hv9%Jj>#Xl<;tnN7f!g<BM)<B=86I%)fH_%2@eBK9y5hnQW>oB>p!V9B|l6>4UF
zqZgK@VP4T=c#z5AiRa@~cXIZD(_BCso5dN->f}LNX9f=zc248L`c53Dwu;ro4Wn#h
z+Z$M78ExD!#x}Nc!?+bHr+NUTe6}&%bR^+QS=^Q9;;yWVyRrcV=C{F@j8^7fg4~+X
z#{Bo8ke?C$gcd6EFM@;+)Iz|};Rt3i;UNz$!shEX=ks7z4-clb@!+)II`UfKNr>kS
z!5>Fnw<w=w8(XYRhw$QnL1ER@A^c<hgU|q`kd1{yMa<I0!d)25X0G6B)Ty3MK;&&^
z#=Dqto$+C2Tr!56bIiyrU$Wv;%yNtsA7$>d%ss#?yIHvNbrwEzlv&PN*~&p?`G&0=
zW5wfa<)@G%%)OggqAY8G-PFe{BP@It&+-d1TC<rW?51-peDQV`{@TnkOH3^D`2{TV
z7M59wExFZ%#hA~1AuPz}u^@FX!;ZZ6hOFqnU>?}`Z`9WSt*2o&s{^o{$H_=>1gAad
zRLA(dp~Q@at;|@sTcOOoI}KTJVftmt-0A3>S!l;WIMZz>l~JoPQjL|lwO(%Y(vHzz
zYRBjZR%<_2DGsPBKf@uqM_3_SRrWMEp21@+JoahB-?|LDcnsSOO@TWMJGfzoE03bG
zhIXFUd^iI}=Z_`7OlUr0JGD)v?IjfXDSp-H3>2EPsn^onEV0^aW@c*zOHK^+A`niA
zhIl<`8|q%}`w|<f(_EX_GvNtgMVn@^6G8ZKV$<vp`{cv|L7@XHr9;6sx*xD&sbB}C
z6F@#*(AYhxxrxnfq&kOhETis*+#}H*EMT#p9&oS^&hy|3yop|-2%9#-gYJ<uXV7*q
z4U707;ZeX1w1Qi0Z@A6H_Kd?0Lr5%WLmQM%yt54LSlRZbZsuPOyY>t$YdZ?IjKg9{
z3sB;!BflX>$gGL^og$73&9%;8+CoQg$%7pDSsn|{ZlwaEXvx^}U|K7}YRHw>qD9d$
za>bE%kzzdW3yN%eQy23O!m2QxZyHJ-H7Z7X32zql47SC@TBGO3j0>6ZZzz7UD1O`&
zKYgqiF>{=v$bE^qhnVF63m-)Mj39mvBVG<+Z%1~K6<<QO7io~W51@`^9b-3rE+S|M
z5p)U>G$bPEJPVH@f(8&l)rg=M5J7jc%xCln3L}Er5kWN~f(psg9!K86hJ{g@Pp-CY
z2@C~leYfTw(F&Q#b>e-7Cw5eetwnR4YK+f2HhI2yr)QnT&3!fvfl-)w8TZ+l6ca@v
zvWlQE3bm2cZffjC<L{5DA`w$Xh^bb@RG7#1G;DMk_BmpRq!ulh*lRE2fFF%qh{bRY
z6aG349oqS;M<|;o<Tl7ONjm+qpFfFZB2K;YTKaMlGx9de^Dh72pLczy5cD6413$e+
z-w}bbNl=l<Kj4?Fq<kI9$m<Yr8VgLn2WL7SxSs#;O+kek?goPzMT?(W-PTpsXRMx8
zZCkL2@#LmYncTE~f8+sJnZ&O$MDV8{B3>WBk)_EjK06_P^~rSTCh@C(CymChDAH5-
zRllfhGLticU;VSLXeNdFzBW1aSLSyiPv%z@-^#BlzLQ@ST#a89d_TV;?}=Yg{{N3(
zRS16dB!1zX;8(;>Ch@ELcYqz!Npa~e6#DWK5kC>{o5ay;5}(=(GlEZTL4mz><~#Tl
zVSgo`(zFRaWt!kqA<_O#!QqHc(NOHI|BmyiVhokwQ@yAZe2V6k!l$<5gNF$|Rhz=6
zh$Ru9;yRx?1Ub&9LX-GZ7izwhPpwbzDg6V76iSt?Q>rF&BaJZeVcKMV#pC>nPQG#e
z0FnWI#a_XgV1i}!k`Wy5>Bsas%lZn>SDbWaLx+X~MB!@gxlNO54v=)ei6mBAWiGga
zfQ><a4Fqf?PQb)?pi*5#rM{he4aB(@rl4`Jwrk^FN6~Ejlt#Q9*x#;Byy$u9#EaK2
z)QDFyy5KX{!oPMB|Jrjk{-xF6m5mSmR!&a?$jEVQWR~=sA+W5L$iWPaewhUQTA&}j
zhi%L_Cm*X+0ROC3uq>@cA$kiKKVHU7a6ZAevcHRO(dxT)z9q)q4!01`IqtX*C!l~c
zmZhgrVE$LI#)xjch+NREY7`Q5s{s<xtqw8rA4IoqfvL%KD;pYkqR}lBw9_;lKGvra
zE{$)A^Wq{pFP71HaThDzLubU|{d5Eje}<Fd-lI4aevY({xi2&K8D_x=a0C?tRB;Xo
zAK`Z*>i`u)hNI*~WJJzD#?GUTW$l40-8k6yWa75>IJ4{li&N-$nsu6mqh@L}69XS-
z;S=|<@R!KPZ)Y>%!7&^&SK^pif@9`98_WE03R8PfyrISP1XH^;+74XtvdfIK^@%Ov
z58`W>O%JHq$|HhvUBTRR$`&VZ8sF1<eE)gVeg%hPe)_yo=W$rT+U7Sw<p8I);j!0d
z5t;jsi?LX<;BOy-zv0-M{sI2VF7kKS)ak5k606e*m57=9P=JMdeu<n`bk2}y_dOhr
zqv@<J+Nle!BypNIFq*RRzKh3Yr|c&aeC}pM@ip@~p2X)0zMIdf2oAyL8YxPW_#DWa
zrb_U+Z03Ixt%$k(AMm-mVd@G#hf(K=A-0+uX)j2jbcj#zxhE$+52>Q~Jj#j>Q}h&f
z(q|*#69~iKQW$=Pw41s2GxrDx9T*xahNyxf5-|v#NmTCuvc0Uh2N}gEsNOiU3?K;c
z<UT`aeopau6!D3u96)@YV_9ce_`I1uf<ZLmbCqyE;<FF&c@Xgl4~`=~y@=0h#ODKu
z&r}|$#pfn{hJp${7(I;XCeb|;?N8Z#b2TcM-8U<_F0n|ife%h2K3FLDpi}U{tV{~p
zU3A-r$aM;4SSXmGXy6t-&<Y>ApxZ=Sz(E40EWXCXYw84EPaUM|M%zz?EsEN*WR}SM
z>qXS2FvUBFDbCQC;>2#`P4L6v^F$&O@(IWU37zJzjW5dc`|tbxh|U+qF&WhGT?)X!
zCPA$t&kJFzMjnf4t%g7XzY(!2_##m|!54}3U4<`pPi_jnh@Do6<^0-8_pRdoiMWTv
zy;a<|i+fnyJ8*Xf=j=pgeJYkiJ6CfBY$(xFSYZ`LDt%5R3zfR5bU&3ksDw^R^b;!4
zSLM+Usq_q${*_8qRC*UB<jZ#NaxTN4UoJbsc7`C+#ICsFI_%2gisx~CJjW`JQHmR&
z&u-{S&f{g$qv@mzsz&&cPj&hQ-<0Ik7_Xole@b1x{-K;Q(+ee44doBzcyVkbU+}H0
zxRP`7nf|_g5pqh_ABiTSj5j9!#&xQat^I}Sl-{X|AwLt!Z2T2|_MOx<K1HpQyh{Wm
z0ulj<fJ8tdAQ6xVNCYGT5&?;TL_i`S5s(N-1SA3y0f~S_Kq4R!kO)WwBmxoviGV~v
zA|Mfv2uK7Z0ulj<fJ8tdAQ6xVNCYGT5&?;TL_i`S5s(N-1SA3y0f~S_Kq4R!kO)Ww
zBmxoviGV~vBJf8;pb`IP0R2Y*Z|M0ZJ^y<>e^1Z<Q_r{Rd05Z?XpYN)B?1xwiGV~v
zA|Mfv2uK7Z0ulj<fJ8tdAQ6xVNCYGT5&?;TL_i`S5s(N-1SA3y0f~S_Kq4R!kO)Ww
zBmxoviGV~vA|Mfv2uK7Z0ulj<fJ8tdAQ6xVNCYGT5&?;TL_i`S5s(N-1SA3y0f~S_
zKqBxbgh1vUN-<3*GvDL!)YX@iKIbW|c+T@|L&>@_<;Uxk=KQh^Rb{AYSmRk!SzT5R
z6}7L^daP}zsw=$@s=AC63*M5-D$*OFFD+kR>m?x@f_GhcIh9Q)uPv#nswpMOLQ?tK
zwKXq@$_1KOS6)w*iy_sntE{dkBl>0VDCtVNXKjV2e0_OoLw&iYqP&FMn_sFdTw79E
ztvu)T)+@D@xHU|nPE%@>rB&r6)rwb{vbMIWe7&oos=l(es<O1CzOtrTaaX@kQ-$G6
zrYc=am8nW&K5o5q+iazq9k=Gy`tV&#H&0nz(3sYHtESkBY{vD9lHaKPe+o1rGZ;lB
z{r?KHgsdp5O<MgNUGBO~lNaf-qSxcUS3rM7*?hZJ|A;Od7icn%*X!+<Lq>D_(=K|b
z$?{*P$S<eJO)2slDe}e?`Ijm3Kc>k4mLg|B20Y4*NYjv}Bh5f!NQBZ2NSR3A*UqGP
zuJX)1rQ*SQN^7q6V4>C@t%sI`hn57E9+n8!yuN}bt$wj~Cuy{J&=<0{gyV~MV%buG
Y7HoW}UJ*DRk5FMHzfe+Nj>+QsA0czxV*mgE

literal 0
HcmV?d00001

diff --git a/convolution/fft-dif.c b/conv/fft-dif.c
similarity index 100%
rename from convolution/fft-dif.c
rename to conv/fft-dif.c
diff --git a/convolution/fft.h b/conv/fft.h
similarity index 100%
rename from convolution/fft.h
rename to conv/fft.h
diff --git a/convolution/ift-dif.c b/conv/ift-dif.c
similarity index 100%
rename from convolution/ift-dif.c
rename to conv/ift-dif.c
diff --git a/convolution/ift.h b/conv/ift.h
similarity index 100%
rename from convolution/ift.h
rename to conv/ift.h
diff --git a/conv/test.c b/conv/test.c
new file mode 100644
index 0000000..65f6d6b
--- /dev/null
+++ b/conv/test.c
@@ -0,0 +1,30 @@
+#include <stdio.h>
+#include "convolution.h"
+
+int main()
+{
+	float complex dirac1[128] = {0.0f};
+	float complex dirac2[128] = {0.0f};
+
+	dirac1[2] = 0.8f;
+	dirac2[2] = 0.5f;
+
+	convolve(dirac1, dirac2, 128);
+
+	for(int i = 0; i < 10; i++)
+	{
+		printf("index:%d value:%f\n", i, creal(dirac1[i]));
+	}
+
+	clean(dirac2, 128);
+
+	dirac2[2] = 0.5f;
+	convolve(dirac1, dirac2, 128);
+
+	for(int i = 0; i < 10; i++)
+	{
+		printf("index:%d value:%f\n", i, creal(dirac1[i]));
+	}
+
+	return 0;
+}
\ No newline at end of file
diff --git a/convolution/convolution.c b/convolution/convolution.c
deleted file mode 100644
index 2e06049..0000000
--- a/convolution/convolution.c
+++ /dev/null
@@ -1,17 +0,0 @@
-#include "convolution.h"
-#include "fft.h"
-#include "ift.h"
-
-void complexMultiplication(float complex* sample, float complex* ir, float complex* result, int size)
-{
-    for(int i = 0; i < size; ++i)
-        *result++ = *sample++ * *ir++;
-}
-
-void convolve(float complex* sample, float complex* ir, float complex* result, int size)
-{
-	fft(sample, size);
-	fft(ir, size);
-    complexMultiplication(sample, ir, result, size);
-    ift(result, size);
-}
diff --git a/convolution/convolution.h b/convolution/convolution.h
deleted file mode 100644
index f4789f6..0000000
--- a/convolution/convolution.h
+++ /dev/null
@@ -1,4 +0,0 @@
-#include <complex.h>
-
-void convolve(float complex* sample, float complex* ir, float complex* result, int size);
-
diff --git a/convolution/test.c b/convolution/test.c
deleted file mode 100644
index 1621bad..0000000
--- a/convolution/test.c
+++ /dev/null
@@ -1,21 +0,0 @@
-#include <stdio.h>
-#include "convolution.h"
-
-int main()
-{
-	float complex dirac1[128] = {0.0f};
-	float complex dirac2[128] = {0.0f};
-	float complex result[128] = {0.0f};
-
-	dirac1[2] = 0.8f;
-	dirac2[2] = 0.5f;
-
-	convolve(dirac1, dirac2, result, 128);
-
-	for(int i = 0; i < 10; i++)
-	{
-		printf("index:%d value:%f\n", i, creal(result[i]));
-	}
-
-	return 0;
-}
\ No newline at end of file

From a94a1bd0cd4dbc1a230d472840267f00043ff5a2 Mon Sep 17 00:00:00 2001
From: Chikashi Miyama <chikashimiyama@gmail.com>
Date: Sat, 3 Apr 2021 22:55:01 +0200
Subject: [PATCH 3/4] unit test for multiple blocksize convolution

---
 conv/convolution.c | 46 ++++++++++++++++++++--------
 conv/convolution.h |  4 +--
 conv/test.c        | 75 ++++++++++++++++++++++++++++++++--------------
 3 files changed, 88 insertions(+), 37 deletions(-)

diff --git a/conv/convolution.c b/conv/convolution.c
index 07707ca..e1e5247 100644
--- a/conv/convolution.c
+++ b/conv/convolution.c
@@ -2,22 +2,44 @@
 #include "fft.h"
 #include "ift.h"
 
-void complexMultiplication(float complex* sample, float complex* ir, int size)
+void copyToComplex(float* buffer, int size, float complex* complexBuffer)
 {
-    while(size--)
-        *sample++ *= *ir++;
+	for(int i = 0; i < size; ++i)
+		complexBuffer[i] = buffer[i];
 }
 
-void convolve(float complex* sample, float complex* ir, int size)
+void copyFromComplex(float complex* complexBuffer, int size, float* buffer)
 {
-	fft(sample, size);
-	fft(ir, size);
-    complexMultiplication(sample, ir, size);
-    ift(sample, size);
+	for(int i = 0; i < size; ++i)
+		buffer[i] = creal(complexBuffer[i]);
 }
 
-void clean(float complex* buffer, int size)
+void convolve(float* sample, int size, float* ir, float* result)
 {
-	while(size--)
-		*buffer++ = 0.f + 0.fI;
-}
\ No newline at end of file
+	// fft ir filter
+	float complex complexIr[128] = {0.f};
+	copyToComplex(ir, 128, complexIr);
+	fft(complexIr, 128);
+
+	float complex complexResult[128] = {0.f};
+	int remainingSamples = size;
+
+	do
+	{
+		float complex complexSample[128] = {0.f};
+		int samplesToBeCopied = remainingSamples > 128 ? 128 : remainingSamples; 
+		copyToComplex(sample, samplesToBeCopied, complexSample);
+		fft(complexSample, 128);
+
+		for(int i = 0; i < 128; ++i)
+    	    complexResult[i] = complexSample[i] * complexIr[i];
+
+	    ift(complexResult, 128);
+    	copyFromComplex(complexResult, samplesToBeCopied, result);
+
+		remainingSamples -= 128;
+		sample += 128;
+		result += 128;
+	}while(remainingSamples > 0);
+
+}
diff --git a/conv/convolution.h b/conv/convolution.h
index 784d740..22c8c54 100644
--- a/conv/convolution.h
+++ b/conv/convolution.h
@@ -1,5 +1,3 @@
 #include <complex.h>
 
-void convolve(float complex* sample, float complex* ir, int size);
-
-void clean(float complex* buf, int size);
+void convolve(float* sample, int size, float* ir, float* result);
diff --git a/conv/test.c b/conv/test.c
index 65f6d6b..6465b4e 100644
--- a/conv/test.c
+++ b/conv/test.c
@@ -1,30 +1,61 @@
+#include <assert.h>
+#include <math.h>
 #include <stdio.h>
 #include "convolution.h"
 
-int main()
+void assert_eq(float expected, float actual)
 {
-	float complex dirac1[128] = {0.0f};
-	float complex dirac2[128] = {0.0f};
-
-	dirac1[2] = 0.8f;
-	dirac2[2] = 0.5f;
-
-	convolve(dirac1, dirac2, 128);
-
-	for(int i = 0; i < 10; i++)
-	{
-		printf("index:%d value:%f\n", i, creal(dirac1[i]));
-	}
+	float dif = fabsf(expected -actual);
+	assert(dif < 0.001f);
+}
 
-	clean(dirac2, 128);
-
-	dirac2[2] = 0.5f;
-	convolve(dirac1, dirac2, 128);
-
-	for(int i = 0; i < 10; i++)
-	{
-		printf("index:%d value:%f\n", i, creal(dirac1[i]));
-	}
+int main()
+{
+	float sample64[64] = {0.0f};
+	float sample128[128] = {0.0f};
+	float sample512[512] = {0.0f};
+	float sample2048[2048] = {0.0f};
+
+	float result64[64] = {0.0f};
+	float result128[128] = {0.0f};
+	float result512[512] = {0.0f};
+	float result2048[2048] = {0.0f};
+
+	float ir[128] = {0.0f};
+
+	// the second last sample have dirac
+	sample64[62] = 1.f;
+	sample128[126] = 1.f;
+	sample512[510] = 1.f;
+	sample2048[2046] = 1.f;
+
+	ir[1] = 0.3f; // generate one sample delay + attenuation
+
+	//expect convolution result in buffers with one inpulse (amp = 0.3) at the end of the buffer
+
+	convolve(sample64, 64, ir, result64);
+	for(int i = 0; i < 63; ++i)
+		assert_eq(result64[i], 0.0f);
+	assert_eq(result64[63], 0.3f);
+	puts("block size 64: OK");
+
+	convolve(sample128, 128, ir, result128);
+	for(int i = 0; i < 127; ++i)
+		assert_eq(result128[i], 0.0f);
+	assert_eq(result128[127], 0.3f);
+	puts("block size 128: OK");
+
+	convolve(sample512, 512, ir, result512);
+	for(int i = 0; i < 511; ++i)
+		assert_eq(result512[i], 0.0f);
+	assert_eq(result512[511], 0.3f);
+	puts("block size 512: OK");
+
+	convolve(sample2048, 2048, ir, result2048);
+	for(int i = 0; i < 2047; ++i)
+		assert_eq(result2048[i], 0.0f);
+	assert_eq(result2048[2047], 0.3f);
+	puts("block size 2048: OK");
 
 	return 0;
 }
\ No newline at end of file

From e8743a3a0f0bc9bfd80fe03069994c0a64ea8ee2 Mon Sep 17 00:00:00 2001
From: Chikashi Miyama <chikashimiyama@gmail.com>
Date: Sat, 3 Apr 2021 22:55:23 +0200
Subject: [PATCH 4/4] replace direct convolution with frequency domain one

---
 Makefile   |  3 ++-
 earplug~.c | 75 +++++++++++++++++++++++++++++++-----------------------
 2 files changed, 45 insertions(+), 33 deletions(-)

diff --git a/Makefile b/Makefile
index f232410..37fc0eb 100644
--- a/Makefile
+++ b/Makefile
@@ -3,7 +3,8 @@
 lib.name = earplug~
 
 cflags = -I./conv
-class.sources = earplug~.c conv/convolution.c conv/fft-dif.c conv/ift-dif.c
+common.sources = conv/convolution.c conv/fft-dif.c conv/ift-dif.c
+class.sources = earplug~.c 
 
 datafiles = earplug~-help.pd earplug~-meta.pd earplug_data_src.txt \
             parse-to-h.pl README.txt LICENSE.txt
diff --git a/earplug~.c b/earplug~.c
index 0714d67..9ebdd15 100644
--- a/earplug~.c
+++ b/earplug~.c
@@ -5,10 +5,13 @@
 /* Revised in spring 2009 by Hans-Christoph Steiner to compile in the data file */
 /* Updated in fall 2020 by Dan Wilcox */
 
-#include "earplug~.h"
 #include <math.h>
 #include <string.h>
 #include <errno.h>
+#include <complex.h>
+
+#include "earplug~.h"
+#include "convolution.h"
 
 #define VERSION "0.2.2"
 
@@ -36,6 +39,9 @@ typedef struct _earplug
     t_float azi;
     t_float ele;
      
+    t_float inputBuffer[8192];
+    t_float outputBuffer[2][8192];
+
     t_float crossCoef[8192];
     t_float azimScale[13];
     unsigned int azimOffset[13];
@@ -48,6 +54,12 @@ typedef struct _earplug
     int bufferPin;
 } t_earplug;
 
+void copy(float* source, float* dist, int size)
+{
+    for(int i = 0; i < size; ++i)
+        *dist++ = *source++;
+}
+
 static t_int *earplug_perform(t_int *w)
 {
     t_earplug *x = (t_earplug *)(w[1]);
@@ -143,39 +155,38 @@ static t_int *earplug_perform(t_int *w)
         }
     }
 
-    float inSample;
-    float convSum[2]; /* to accumulate the sum during convolution */
-    int blockScale = 8192 / blocksize;
+    convolve(in, blocksize, x->currentImpulse[0], left_out);
+    convolve(in, blocksize, x->currentImpulse[1], right_out);
 
     /* convolve the interpolated HRIRs (left and right) with the input signal */
-    while (blocksize--)
-    {
-        convSum[0] = 0; 
-        convSum[1] = 0; 
-
-        inSample = *(in++);
-
-        x->convBuffer[x->bufferPin] = inSample;
-        unsigned scaledBlocksize = blocksize * blockScale;
-        unsigned blocksizeDelta = 8191 - scaledBlocksize;
-        for (i = 0; i < 128; i++)
-        { 
-            convSum[0] += (x->previousImpulse[0][i] * x->crossCoef[blocksizeDelta] + 
-                            x->currentImpulse[0][i] * x->crossCoef[scaledBlocksize]) *
-                            x->convBuffer[(x->bufferPin - i) &127];
-                            
-            convSum[1] += (x->previousImpulse[1][i] * x->crossCoef[blocksizeDelta] +
-                            x->currentImpulse[1][i] * x->crossCoef[scaledBlocksize]) *
-                            x->convBuffer[(x->bufferPin - i) &127];
-
-            x->previousImpulse[0][i] = x->currentImpulse[0][i];
-            x->previousImpulse[1][i] = x->currentImpulse[1][i];
-        }   
-        x->bufferPin = (x->bufferPin + 1) & 127;
-
-        *left_out++ = convSum[0];
-        *right_out++ = convSum[1];
-    }
+//    while (blocksize--)
+//    {
+//        convSum[0] = 0;
+//        convSum[1] = 0;
+//
+//        inSample = *(in++);
+//
+//        x->convBuffer[x->bufferPin] = inSample;
+//        unsigned scaledBlocksize = blocksize * blockScale;
+//        unsigned blocksizeDelta = 8191 - scaledBlocksize;
+//        for (i = 0; i < 128; i++)
+//        {
+//            convSum[0] += (x->previousImpulse[0][i] * x->crossCoef[blocksizeDelta] +
+//                            x->currentImpulse[0][i] * x->crossCoef[scaledBlocksize]) *
+//                            x->convBuffer[(x->bufferPin - i) &127];
+//
+//            convSum[1] += (x->previousImpulse[1][i] * x->crossCoef[blocksizeDelta] +
+//                            x->currentImpulse[1][i] * x->crossCoef[scaledBlocksize]) *
+//                            x->convBuffer[(x->bufferPin - i) &127];
+//
+//            x->previousImpulse[0][i] = x->currentImpulse[0][i];
+//            x->previousImpulse[1][i] = x->currentImpulse[1][i];
+//        }
+//        x->bufferPin = (x->bufferPin + 1) & 127;
+//
+//        *left_out++ = convSum[0];
+//        *right_out++ = convSum[1];
+//    }
     return w + 6;
 }