# HG changeset patch
# User Ben Abbott <bpabbott@mac.com>
# Date 1233765469 18000
# Node ID a891981681758a896e7978d903c508c6924436a6
# Parent  34b7b93f91bcc0a42f17f7d2993a8bf189d8b07b
freqz.m: freqz.m: fix for long input

diff -r 34b7b93f91bc -r a89198168175 scripts/ChangeLog
--- a/scripts/ChangeLog	Wed Feb 04 11:28:24 2009 -0500
+++ b/scripts/ChangeLog	Wed Feb 04 11:37:49 2009 -0500
@@ -1,3 +1,8 @@
+2009-02-04  Frederick Umminger  <Frederick_Umminger@playstation.sony.com>
+
+	* signal/freqz.m: Ensure causal phase response.
+	Handle long input correctly.
+
 2009-02-04  Petr Mikulik  <mikulik@physics.muni.cz>
 
 	* plot/__go_draw_axes__.m: Pass "interpolate 0, 0" to gnuplot
diff -r 34b7b93f91bc -r a89198168175 scripts/signal/freqz.m
--- a/scripts/signal/freqz.m	Wed Feb 04 11:28:24 2009 -0500
+++ b/scripts/signal/freqz.m	Wed Feb 04 11:37:49 2009 -0500
@@ -127,19 +127,38 @@
     k = max (length (b), length (a));
     hb = polyval (postpad (b, k), exp (j*w));
     ha = polyval (postpad (a, k), exp (j*w));
-  elseif (strcmp (region, "whole"))
-    f = Fs * (0:n-1)' / n;
+  else
     ## polyval(fliplr(P),exp(jw)) is O(p n) and fft(x) is O(n log(n)),
     ## where p is the order of the polynomial P.  For small p it
     ## would be faster to use polyval but in practice the overhead for
     ## polyval is much higher and the little bit of time saved isn't
     ## worth the extra code.
-    hb = fft (postpad (b, n));
-    ha = fft (postpad (a, n));
-  else
-    f = Fs/2 * (0:n-1)' / n;
-    hb = fft (postpad (b, 2*n))(1:n);
-    ha = fft (postpad (a, 2*n))(1:n);
+    k = max (length (b), length (a));
+    if (k > n/2 && nargout == 0)
+      ## Ensure a causal phase response.
+      n = n * 2 .^ ceil (log2 (2*k/n));
+    endif
+
+    if (strcmp (region, "whole"))
+      N = n;
+    else
+      N = 2*n;
+    endif
+
+    f = Fs * (0:n-1).' / N;
+
+    pad_sz = N*ceil (k/N);
+    b = postpad (b, pad_sz);
+    a = postpad (a, pad_sz);
+
+    hb = zeros (n, 1);
+    ha = zeros (n, 1);
+
+    for i = 1:N:pad_sz
+      hb = hb + fft (postpad (b(i:i+N-1), N))(1:n);
+      ha = ha + fft (postpad (a(i:i+N-1), N))(1:n);
+    endfor
+
   endif
 
   h = hb ./ ha;