diff --git a/bfps/test/test_fftw.py b/bfps/test/test_fftw.py
new file mode 100644
index 0000000000000000000000000000000000000000..3de54502a47c5cb9dcc6d4298ea43a7705c1b548
--- /dev/null
+++ b/bfps/test/test_fftw.py
@@ -0,0 +1,56 @@
+#! /usr/bin/env python
+
+import os
+import numpy as np
+#import cupy as np
+import h5py
+import sys
+
+import bfps
+from bfps import TEST
+
+import matplotlib.pyplot as plt
+
+def main():
+    niterations = 10
+    c = TEST()
+    c.launch(
+            ['symmetrize_test',
+             '--nx', '32',
+             '--ny', '32',
+             '--nz', '32',
+             '--np', '4',
+             '--ntpp', '1',
+             '--wd', './'] +
+             sys.argv[1:])
+
+    df = h5py.File(c.simname + '.h5', 'r')
+    print(df['kspace/kx'].value)
+    print(df['kspace/ky'].value)
+    print(df['kspace/kz'].value)
+    df = h5py.File(c.simname + '_fields.h5', 'r')
+    field1_complex = df['field1/complex/0'].value
+    field1_real = df['field1/real/0'].value
+
+    np_field1_real = np.fft.irfftn(field1_complex, axes = (0, 1, 2)).transpose(1, 0, 2, 3)
+    print(np.mean(field1_real**2))
+    print(np.mean(np_field1_real**2)*(np_field1_real.size/3)**2)
+    print(np.max(np.abs(field1_real - np_field1_real*(np_field1_real.size/3))))
+
+    np_field1_complex = np.fft.rfftn(field1_real.transpose(1, 0, 2, 3), axes = (0, 1, 2))
+
+    print(np.sum(np.abs(field1_complex)**2))
+    print(np.sum(np.abs(np_field1_complex)**2))
+    print(np.max(np.abs(np_field1_complex - field1_complex)))
+
+    f = plt.figure()
+    a = f.add_subplot(121)
+    a.imshow(np.log(np.abs(np_field1_complex[:, :, 0, 0])), interpolation = 'nearest')
+    a = f.add_subplot(122)
+    a.imshow(np.log(np.abs(field1_complex[:, :, 0, 0])), interpolation = 'nearest')
+    f.savefig('symmetrize_test.pdf')
+    return None
+
+if __name__ == '__main__':
+    main()
+