# HG changeset patch
# User jwe
# Date 1021580897 0
# Node ID 2ca2d23a49a739db1058b51d9290225a5a734a7e
# Parent  fab8337340a10f131cd88930b86b3a00ae37eb12
[project @ 2002-05-16 20:27:59 by jwe]

diff --git a/scripts/signal/freqz_plot.m b/scripts/signal/freqz_plot.m
new file mode 100644
--- /dev/null
+++ b/scripts/signal/freqz_plot.m
@@ -0,0 +1,100 @@
+## Copyright (C) 2002 John W. Eaton
+##
+## This file is part of Octave.
+##
+## Octave 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, or (at your option)
+## any later version.
+##
+## Octave 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 Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
+## 02111-1307, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} freqz_plot (@var{w}, @var{h})
+## Plot the pass band, stop band and phase response of @var{h}.
+## @end deftypefn
+
+## Author: Paul Kienzle <pkienzle@users.sf.net>
+
+function freqz_plot(w,h)
+
+    n = length (w);
+
+    ## ## exclude zero-frequency
+    ## h = h (2 : length (h));
+    ## w = w (2 : length (w));
+    ## n = n-1;
+
+    mag = 20 * log10 (abs (h));
+    phase = unwrap (arg (h));
+    maxmag = max (mag);
+
+    unwind_protect
+
+      ## Protect graph state.
+
+      if (gnuplot_has_multiplot)
+      	subplot (311);
+      	gset lmargin 10;
+	axis ("labely");
+	xlabel ("");
+      endif
+      grid ("on");
+      axis ([ w(1), w(n), maxmag-3, maxmag ]);
+      plot (w, mag, ";Pass band (dB);");
+
+      if (gnuplot_has_multiplot)
+      	subplot (312);
+	axis ("labely");
+	title ("");
+	xlabel ("");
+      	gset tmargin 0;
+      else
+	input ("press any key for the next plot: ");
+      endif
+      grid ("on");
+      if (maxmag - min (mag) > 100)
+      	axis ([ w(1), w(n), maxmag-100, maxmag ]);
+      else
+      	axis ("autoy");
+      endif
+      plot (w, mag, ";Stop band (dB);");
+      
+      if (gnuplot_has_multiplot)
+      	subplot (313);
+	axis ("label");
+	title ("");
+      else
+	input ("press any key for the next plot: ");
+      endif
+      grid ("on");
+      axis ("autoy");
+      xlabel ("Frequency");
+      axis ([ w(1), w(n) ]);
+      plot (w, phase/(2*pi), ";Phase (radians/2pi);");
+      
+    unwind_protect_cleanup
+
+      ## Restore graph state.
+
+      ## XXX FIXME XXX -- if automatic_replot is non-zero, this will
+      ## mess up the graph, however if we don't do it here then the user
+      ## will have to do it themselves.
+
+      grid ("off");
+      axis ("auto", "label");
+      gset lmargin;
+      gset tmargin;
+      oneplot ();
+
+    end_unwind_protect
+
+endfunction
diff --git a/scripts/signal/unwrap.m b/scripts/signal/unwrap.m
new file mode 100644
--- /dev/null
+++ b/scripts/signal/unwrap.m
@@ -0,0 +1,99 @@
+## Copyright (C) 2000  Bill Lash
+##
+## This file is part of Octave.
+##
+## Octave 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, or (at your option)
+## any later version.
+##
+## Octave 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 Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
+## 02111-1307, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{b} =} unwrap (@var{a}, @var{tol}, @var{dim})
+## 
+## Unwrap radian phases by adding multiples of 2*pi as appropriate to
+## remove jumps greater than @var{tol}.  @var{tol} defaults to pi.
+##
+## Unwrap will unwrap along the columns of @var{a} unless the row
+## dimension of @var{a} is 1 or @var{dim} is given with a
+## value of 1, when it will unwrap along the row(s).
+## @end deftypefn
+
+## Author: Bill Lash <lash@tellabs.com>
+
+function retval = unwrap (a, tol, dim)
+        
+  if (nargin < 1 || nargin > 3)
+    usage ("unwrap (a [, tol [, dim]])")
+  endif
+
+  if (nargin < 3)
+    if (rows (a) == 1)
+      ## Row vector, go along the row.
+      dim = 1;
+    else
+      ## Otherwise go along the columns.
+      dim = 2;
+    endif
+  endif
+
+  if (nargin < 2)
+    tol = pi;
+  endif
+
+  ## If TOL is not provided but dim is, handle it.
+  if (tol == [])
+    tol = pi;
+  endif
+
+  ## Don't let anyone use a negative value for TOL.
+  tol = abs (tol);
+  
+  rng = 2*pi;
+  
+  ## Put data in a form so that we can unwrap each column.
+  if (dim == 1)
+    ra = a.';
+  else
+    ra = a;
+  endif
+  n = columns (ra);
+  m = rows (ra);
+
+  ## Handle case where we are trying to unwrap a scalar, or only have
+  ## one sample in the specified dimension.
+  if (m == 1)       
+    retval = a;     
+    return;         
+  endif
+
+  ## Take first order difference to see so that wraps will show up
+  ## as large values, and the sign will show direction.
+  d = [ zeros(1,n); ra(1:m-1,:)-ra(2:m,:) ];
+
+  ## Find only the peaks, and multiply them by the range so that there
+  ## are kronecker deltas at each wrap point multiplied by the range
+  ## value.
+  p =  rng * (((d > tol) > 0) - ((d < -tol) > 0));
+
+  ## Now need to "integrate" this so that the deltas become steps.
+  r = cumsum (p);
+
+  ## Now add the "steps" to the original data and put output in the
+  ## same shape as originally.
+  if (dim == 1)
+    retval = (ra + r).';
+  else
+    retval = (ra + r);
+  endif
+
+endfunction