Re: outstanding changesets

[John W. Eaton]

On 28-Jan-2009, Ben Abbott wrote:

| 
| On Jan 27, 2009, at 11:03 PM, John W. Eaton wrote:
| 
| > On  2-Jan-2009, Ben Abbott wrote:
| >
| > | I have four outstanding changesets. They are below.
| > |
| > *snip*
| > |
| > | (4) freqz produces incorrect output
| > | http://www-old.cae.wisc.edu/pipermail/bug-octave/2008-December/007493.html
| > | changeset-freqz.patch
| >
| > Have we decided that this is the correct fix, or is there more to it?
| >
| > Thanks,
| >
| > jwe
| 
| Regarding freqz, you had pointed out the phase looked strange. I'm not  
| sure I found what you saw, but it appears that the proprietry  
| implementation unwraps the phase in only one direction. This ensures  
| the phase delay is positive (phase slope is negative). I don' t have a  
| nice solution for that yet.
| 
| The phase slope issue is unrelated to the original problem, so there  
| is some merit to committing the current change and following up with a  
| second. Do you have a preference? (anyone else?).

OK, for some random data (attached) here is what I see (attached pdf
file) with the current version of the patch that I have (also
attached).

I'm just doing

  load foo.dat
  freqz (foo)

Do you see the same results for this data file?  Does Matlab produce
the same result?

jwe

foo.dat.gz
Description: Binary data

foo.pdf
Description: Adobe PDF document

diff --git a/scripts/ChangeLog b/scripts/ChangeLog
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -359,6 +359,10 @@
 
        * testfun/test.m: Print "has no tests" message if there are demos
        but no tests instead of printing PASSES 0 out of 0 tests.
+
+2008-12-24  Frederick Umminger <address@hidden>
+
+       * signal/freqz.m: freqz.m: fix for long input
 
 2008-12-23  David Bateman  <address@hidden>
 
diff --git a/scripts/signal/freqz.m b/scripts/signal/freqz.m
--- a/scripts/signal/freqz.m
+++ b/scripts/signal/freqz.m
@@ -127,19 +127,33 @@
     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);
+    if (strcmp (region, "whole"))
+      N= n;
+    else
+      N = 2*n;
+    endif
+
+    f = Fs * (0:n-1)' / N;
+
+    sz = max(size(b,1), size(a,1));
+    pad_sz = N*ceil(sz/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;