diff --git a/unit_tests/avwutils/fslcomplex/feedsRun b/unit_tests/avwutils/fslcomplex/feedsRun
new file mode 100755
index 0000000000000000000000000000000000000000..dd219a522f00bf149c1646a33def9a8f4663d749
--- /dev/null
+++ b/unit_tests/avwutils/fslcomplex/feedsRun
@@ -0,0 +1,114 @@
+#!/usr/bin/env python
+"""Test that fslcomplex produces correct outputs. Specifically tests that the
+output files have orientation info that matches that of the input files.
+"""
+
+import numpy      as np
+import nibabel    as nib
+import subprocess as sp
+import os.path    as op
+import               shlex
+import               sys
+import               traceback
+
+from fsl.utils.tempdir import tempdir
+from fsl.data.image    import addExt
+
+
+def sprun(cmd):
+    print(f'RUN {cmd}')
+    sp.run(shlex.split(cmd), check=True)
+
+
+# radio=[True|False] -> whether or not
+# to add a flip on the affine X axis
+def create_test_input_data(radio):
+
+    real  = np.random.randint(0, 100, (10, 10, 10)).astype(np.float32)
+    imag  = np.random.randint(0, 100, (10, 10, 10)).astype(np.float32)
+    comp  = real + imag * 1j
+    abso  = np.abs(comp)
+    phase = np.arctan2(comp.imag, comp.real)
+
+    affine        = np.diag([3, 3, 3, 0])
+    affine[:3, 3] = [10, 20, 30]
+
+    if radio:
+        affine[0, 0] *= -1
+
+    real  = nib.Nifti1Image(real,  affine)
+    imag  = nib.Nifti1Image(imag,  affine)
+    comp  = nib.Nifti1Image(comp,  affine)
+    abso  = nib.Nifti1Image(abso,  affine)
+    phase = nib.Nifti1Image(phase, affine)
+
+    real .set_qform(*real .get_sform(coded=True))
+    imag .set_qform(*comp .get_sform(coded=True))
+    comp .set_qform(*comp .get_sform(coded=True))
+    abso .set_qform(*abso .get_sform(coded=True))
+    phase.set_qform(*phase.get_sform(coded=True))
+
+    real .to_filename('real_in.nii.gz')
+    imag .to_filename('imag_in.nii.gz')
+    comp .to_filename('complex_in.nii.gz')
+    abso .to_filename('abs_in.nii.gz')
+    phase.to_filename('phase_in.nii.gz')
+
+
+def compare_images(file1, file2):
+    img1  = nib.load(addExt(file1))
+    img2  = nib.load(addExt(file2))
+    data1 = np.asanyarray(img1.dataobj)
+    data2 = np.asanyarray(img2.dataobj)
+    assert img1.header['sform_code'] == img2.header['sform_code']
+    assert img1.header['qform_code'] == img2.header['qform_code']
+    assert np.all(np.isclose(img1.get_sform(),      img2.get_sform()))
+    assert np.all(np.isclose(img1.get_qform(),      img2.get_qform()))
+    assert np.all(np.isclose(img1.header['pixdim'], img2.header['pixdim']))
+    assert np.all(np.isclose(data1,                 data2))
+
+
+def test_fslcomplex_call(args, infiles, outfiles, radio):
+    with tempdir():
+        create_test_input_data(radio)
+        sprun(f'fslcomplex {args} {infiles} {outfiles}')
+        for outfile in outfiles.split(' '):
+            prefix = outfile.split('_')[0]
+            infile = f'{prefix}_in'
+            print(f'Comparing {outfile} against {infile} [prefix {prefix}]')
+            compare_images(infile, outfile)
+
+
+if __name__ == '__main__':
+    tests = [
+        ('-realabs',       'complex_in',      'abs_out',           False),
+        ('-realabs',       'complex_in',      'abs_out',           True),
+        ('-realphase',     'complex_in',      'phase_out',         False),
+        ('-realphase',     'complex_in',      'phase_out',         True),
+        ('-realpolar',     'complex_in',      'abs_out phase_out', False),
+        ('-realpolar',     'complex_in',      'abs_out phase_out', True),
+        ('-realcartesian', 'complex_in',      'real_out imag_out', False),
+        ('-realcartesian', 'complex_in',      'real_out imag_out', True),
+        ('-complex',       'real_in imag_in', 'complex_out',       False),
+        ('-complex',       'real_in imag_in', 'complex_out',       True),
+        ('-complexpolar',  'abs_in phase_in', 'complex_out',       False),
+        ('-complexpolar',  'abs_in phase_in', 'complex_out',       True),
+        # ('-complexsplit',  [], [], False), # ??
+        # ('-complexmerge',  [], [], False), # ??
+        # ('-copyonly',      [], [], False), # ??
+        # ('-complexsplit',  [], [], True), # ??
+        # ('-complexmerge',  [], [], True), # ??
+        # ('-copyonly',      [], [], True), # ??
+    ]
+
+    result = 0
+    for args, infiles, outfiles, radio in tests:
+        try:
+            test_fslcomplex_call(args, infiles, outfiles, radio)
+            print(f'\nTest fslcomplex {args} {infiles} {outfiles} PASSED')
+        except Exception as e:
+            print(f'\nTest fslcomplex {args} {infiles} {outfiles} FAILED')
+            traceback.print_exc()
+            result = 1
+
+    sys.exit(result)