Mercurial > hg > octave-lojdl
changeset 5377:bd4ee620c38a
[project @ 2005-06-02 15:42:39 by jwe]
author | jwe |
---|---|
date | Thu, 02 Jun 2005 15:42:39 +0000 |
parents | e0b390a01639 |
children | b2a5596a3f7b |
files | scripts/ChangeLog scripts/polynomial/polygcd.m scripts/signal/freqz.m |
diffstat | 3 files changed, 126 insertions(+), 42 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,10 @@ +2005-06-02 Paul Kienzle <pkienzle@users.sf.net> + + * signal/freqz.m: Use correct calculations when given a vector of + frequencies. Improve accuracy of returned frequency vector. + Improve speed for medium length filters (at a slight cost for slow + filters). Add test cases. + 2005-05-27 "Dmitri A. Sergatskov" <dasergatskov@gmail.com> * plot/loglog.m: Fix set commands.
--- a/scripts/polynomial/polygcd.m +++ b/scripts/polynomial/polygcd.m @@ -1,3 +1,68 @@ +## Copyright (C) 2000 Paul Kienzle +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; if not, write to the Free Software +## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{q}]} polygcd (@var{b}, @var{a}, @var{tol}) +## +## Find greatest common divisor of two polynomials. This is equivalent +## to the polynomial found by multiplying together all the common roots. +## Together with deconv, you can reduce a ratio of two polynomials. +## Tolerance defaults to +## @example +## sqrt(eps). +## @end example +## Note that this is an unstable +## algorithm, so don't try it on large polynomials. +## +## Example +## @example +## polygcd(poly(1:8),poly(3:12)) - poly(3:8) +## deconv(poly(1:8),polygcd(poly(1:8),poly(3:12))) - poly(1:2) +## @end example +## @end deftypefn +## +## @seealso{poly, polyinteg, polyderiv, polyreduce, roots, conv, deconv, +## residue, filter, polyval, and polyvalm} + +function x = polygcd(b,a,tol) + if (nargin<2 || nargin>3) + usage("x=polygcd(b,a [,tol])"); + endif + if (nargin<3), tol=sqrt(eps); endif + if (length(a)==1 || length(b)==1) + if a==0, x=b; + elseif b==0, x=a; + else x=1; + endif + return; + endif + a = a./a(1); + while (1) + [d, r] = deconv(b, a); + nz = find(abs(r)>tol); + if isempty(nz) + x = a; + return; + else + r = r(nz(1):length(r)); + endif + b = a; + a = r./r(1); + endwhile +endfunction ## Copyright (C) 2000 Paul Kienzle ## ## This program is free software; you can redistribute it and/or modify
--- a/scripts/signal/freqz.m +++ b/scripts/signal/freqz.m @@ -14,8 +14,8 @@ ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, write to the Free -## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA -## 02110-1301, USA. +## Software Foundation, 59 Temple Place - Suite 330, Boston, MA +## 02111-1307, USA. ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{h}, @var{w}] =} freqz (@var{b}, @var{a}, @var{n}, "whole") @@ -70,7 +70,7 @@ ## Author: jwe ??? -function [h_r, w_r] = freqz (b, a, n, region, Fs) +function [h_r, f_r] = freqz (b, a, n, region, Fs) if (nargin < 1 || nargin > 5) usage ("[h, w] = freqz (b, a, n [, \"whole\"] [, Fs])"); @@ -90,6 +90,9 @@ endif endif + if (isempty (b)) + b = 1; + endif if (isempty (a)) a = 1; endif @@ -111,56 +114,65 @@ endif endif - la = length (a); - a = reshape (a, 1, la); - lb = length (b); - b = reshape (b, 1, lb); - k = max ([la, lb]); + a = a(:).'; + b = b(:).'; - if (! isscalar (n)) - if (nargin == 4) ## Fs was specified - w = 2*pi*n/Fs; - else - w = n; + if (! isscalar (n)) ## Explicit frequency vector given + w = f = n; + if (nargin == 4) ## Sampling rate Fs was specified + w = 2*pi*f/Fs; endif - n = length (n); - extent = 0; + hb = polyval (fliplr(b), exp(-j*w)); + ha = polyval (fliplr(a), exp(-j*w)); elseif (strcmp (region, "whole")) - w = 2 * pi * (0:n-1) / n; - extent = n; + f = Fs * (0:n-1) / n; + ## polyval(fliplr(P),exp(-jw)) is O(p n) and fft(x) is O(n log(n)), + ## where p is the order of the 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 - w = pi * (0:n-1) / n; - extent = 2 * n; + f = Fs/2 * (0:n-1) / n; + hb = fft (postpad (b, 2*n))(1:n); + ha = fft (postpad (a, 2*n))(1:n); endif - if (length (b) == 1) - if (length (a) == 1) - hb = b * ones (1, n); - else - hb = b; - endif - elseif (extent >= k) - hb = fft (postpad (b, extent)); - hb = hb(1:n); - else - hb = polyval (postpad (b, k), exp (j*w)); - endif - if (length (a) == 1) - ha = a; - elseif (extent >= k) - ha = fft (postpad (a, extent)); - ha = ha(1:n); - else - ha = polyval (postpad (a, k), exp (j*w)); - endif h = hb ./ ha; - w = Fs * w / (2*pi); if (nargout != 0), # return values and don't plot h_r = h; - w_r = w; + f_r = f; else # plot and don't return values - freqz_plot (w, h); + freqz_plot (f, h); end endfunction + +%!test # correct values and fft-polyval consistency +%! # butterworth filter, order 2, cutoff pi/2 radians +%! b = [0.292893218813452 0.585786437626905 0.292893218813452]; +%! a = [1 0 0.171572875253810]; +%! [h,w] = freqz(b,a,32); +%! assert(h(1),1,10*eps); +%! assert(abs(h(17)).^2,0.5,10*eps); +%! assert(h,freqz(b,a,w),10*eps); # fft should be consistent with polyval + +%!test # whole-half consistency +%! b = [1 1 1]/3; # 3-sample average +%! [h,w] = freqz(b,1,32,'whole'); +%! assert(h(2:16),conj(h(32:-1:18)),20*eps); +%! [h2,w2] = freqz(b,1,16,'half'); +%! assert(h(1:16),h2,20*eps); +%! assert(w(1:16),w2,20*eps); + +%!test # Sampling frequency properly interpreted +%! b = [1 1 1]/3; +%! [h,f] = freqz(b,1,16,320); +%! assert(f,[0:15]*10,10*eps); +%! [h2,f2] = freqz(b,1,[0:15]*10,320); +%! assert(f2,[0:15]*10,10*eps); +%! assert(h,h2,20*eps); +%! [h3,f3] = freqz(b,1,32,'whole',320); +%! assert(f3,[0:31]*10,10*eps);