# HG changeset patch
# User Mike Miller <mtmiller@ieee.org>
# Date 1375651660 14400
# Node ID 3a23cbde59d5fa1682cc161bf26efdf387529089
# Parent  1d5c0c9b3e990e9bddb3fe083835ae0bc00b63cc
interpft.m: Fix interpolation to preserve spectral symmetry (bug #39566)

* interpft.m: Fix interpolation to preserve spectral symmetry, be compatible
with Matlab. Add test cases.

diff -r 1d5c0c9b3e99 -r 3a23cbde59d5 scripts/general/interpft.m
--- a/scripts/general/interpft.m	Sun Aug 04 14:14:26 2013 -0700
+++ b/scripts/general/interpft.m	Sun Aug 04 17:27:40 2013 -0400
@@ -70,7 +70,7 @@
 
   inc = ceil (m/n);
   y = fft (x) / m;
-  k = floor (m / 2);
+  k = ceil (m / 2);
   sz = size (x);
   sz(1) = n * inc - m;
 
@@ -79,6 +79,18 @@
   z = cat (1, y(idx{:}), zeros (sz));
   idx{1} = k+1:m;
   z = cat (1, z, y(idx{:}));
+
+  ## When m is an even number of rows, the FFT has a single Nyquist bin.
+  ## If zero-padded above, distribute the value of the Nyquist bin evenly
+  ## between the new corresponding positive and negative frequency bins.
+  if (sz(1) > 0 && k == m/2)
+    idx{1} = n * inc - k + 1;
+    tmp = z(idx{:}) / 2;
+    z(idx{:}) = tmp;
+    idx{1} = k + 1;
+    z(idx{:}) = tmp;
+  endif
+
   z = n * ifft (z);
 
   if (inc != 1)
@@ -108,6 +120,18 @@
 %!assert (interpft (y', n), y', 20*eps);
 %!assert (interpft ([y,y],n), [y,y], 20*eps);
 
+%% Test case with complex input from bug #39566
+%!test
+%! x = (1 + j) * [1:4]';
+%! y = ifft ([15 + 15*j; -6; -1.5 - 1.5*j; 0; -1.5 - 1.5*j; -6*j]);
+%! assert (interpft (x, 6), y, 10*eps);
+
+%% Test for correct spectral symmetry with even/odd lengths
+%!assert (max (abs (imag (interpft ([1:8], 20)))), 0, 2*eps);
+%!assert (max (abs (imag (interpft ([1:8], 21)))), 0, 2*eps);
+%!assert (max (abs (imag (interpft ([1:9], 20)))), 0, 2*eps);
+%!assert (max (abs (imag (interpft ([1:9], 21)))), 0, 2*eps);
+
 %% Test input validation
 %!error interpft ()
 %!error interpft (1)