[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: outstanding changesets
From: |
John W. Eaton |
Subject: |
Re: outstanding changesets |
Date: |
Wed, 28 Jan 2009 22:13:20 -0500 |
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;