Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[8792] trunk/octave-forge/extra/control-devel
|
From: <par...@us...> - 2011-10-19 18:16:33
|
Revision: 8792
http://octave.svn.sourceforge.net/octave/?rev=8792&view=rev
Author: paramaniac
Date: 2011-10-19 18:16:26 +0000 (Wed, 19 Oct 2011)
Log Message:
-----------
control-devel: add code for fitfrd (3)
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/devel/makefile_ident.m
trunk/octave-forge/extra/control-devel/inst/fitfrd.m
Added Paths:
-----------
trunk/octave-forge/extra/control-devel/devel/test_fitfrd.m
trunk/octave-forge/extra/control-devel/src/slsb10yd.cc
Modified: trunk/octave-forge/extra/control-devel/devel/makefile_ident.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/makefile_ident.m 2011-10-19 17:01:28 UTC (rev 8791)
+++ trunk/octave-forge/extra/control-devel/devel/makefile_ident.m 2011-10-19 18:16:26 UTC (rev 8792)
@@ -24,7 +24,8 @@
MB03UD.f MB01SD.f
## fit state-space model to frequency response data
-mkoctfile SB10YD.f DG01MD.f AB04MD.f SB10ZP.f AB07ND.f \
+mkoctfile slsb10yd.cc \
+ SB10YD.f DG01MD.f AB04MD.f SB10ZP.f AB07ND.f \
MC01PD.f TD04AD.f TD03AY.f TB01PD.f TB01XD.f \
AB07MD.f TB01UD.f TB01ID.f MB01PD.f MB03OY.f \
MB01QD.f
Added: trunk/octave-forge/extra/control-devel/devel/test_fitfrd.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/test_fitfrd.m (rev 0)
+++ trunk/octave-forge/extra/control-devel/devel/test_fitfrd.m 2011-10-19 18:16:26 UTC (rev 8792)
@@ -0,0 +1,6 @@
+sys = ss (-1, 1, 1, 0);
+n = 1;
+
+ret0 = fitfrd (sys, n)
+
+ret1 = fitfrd (sys, n, 1)
\ No newline at end of file
Modified: trunk/octave-forge/extra/control-devel/inst/fitfrd.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/fitfrd.m 2011-10-19 17:01:28 UTC (rev 8791)
+++ trunk/octave-forge/extra/control-devel/inst/fitfrd.m 2011-10-19 18:16:26 UTC (rev 8792)
@@ -16,9 +16,8 @@
## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
-## @deftypefn{Function File} {[@var{K}, @var{N}, @var{gamma}, @var{info}] =} ncfsyn (@var{G}, @var{W1}, @var{W2}, @var{factor})
-## Normalized Coprime Factor (NCF) H-infinity synthesis.
-## Compute positive feedback controller using the McFarlane/Glover Loop Shaping Design Procedure.
+## @deftypefn{Function File} {[@var{sys}, @var{n}] =} fitfrd (@var{dat}, @var{n})
+## Fit frequency response data with a stable, minimum-phase state-space system.
##
## @strong{Inputs}
## @table @var
@@ -68,8 +67,25 @@
## Created: October 2011
## Version: 0.1
-function [retsys, n] = fitfrd (sys, n, code)
+function [sys, n] = fitfrd (dat, n, flag = 0)
+ if (nargin == 0 || nargin > 3)
+ print_usage ();
+ endif
+ if (! isa (dat, "frd"))
+ dat = frd (dat);
+ endif
+
+ if (! issample (n, 0) || n != round (n))
+ error ("fitfrd: second argument must be an integer >= 0");
+ endif
+ [H, w, tsam] = frdata (dat, "vector");
+ dt = isdt (dat);
+
+ [a, b, c, d, n] = slsb10yd (real (H), imag (H), w, n, dt, flag);
+
+ sys = ss (a, b, c, d, tsam);
+
endfunction
\ No newline at end of file
Added: trunk/octave-forge/extra/control-devel/src/slsb10yd.cc
===================================================================
--- trunk/octave-forge/extra/control-devel/src/slsb10yd.cc (rev 0)
+++ trunk/octave-forge/extra/control-devel/src/slsb10yd.cc 2011-10-19 18:16:26 UTC (rev 8792)
@@ -0,0 +1,159 @@
+/*
+
+Copyright (C) 2011 Lukas F. Reichlin
+
+This file is part of LTI Syncope.
+
+LTI Syncope 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 3 of the License, or
+(at your option) any later version.
+
+LTI Syncope 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 LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
+
+Fit FRD with SS model.
+Uses SLICOT SB10YD by courtesy of NICONET e.V.
+<http://www.slicot.org>
+
+Author: Lukas Reichlin <luk...@gm...>
+Created: October 2011
+Version: 0.1
+
+*/
+
+#include <octave/oct.h>
+#include <f77-fcn.h>
+#include <complex>
+#include "common.cc"
+
+extern "C"
+{
+ int F77_FUNC (sb10yd, SB10YD)
+ (int& DISCFL, int& FLAG,
+ int& LENDAT,
+ double* RFRDAT, double* IFRDAT,
+ double* OMEGA,
+ int& N,
+ double* A, int& LDA,
+ double* B,
+ double* C,
+ double* D,
+ double& TOL,
+ int* IWORK, double* DWORK, int& LDWORK,
+ Complex* ZWORK, int& LZWORK,
+ int& INFO);
+}
+
+DEFUN_DLD (slsb10yd, args, nargout,
+ "-*- texinfo -*-\n\
+Slicot SB10YD Release 5.0\n\
+No argument checking.\n\
+For internal use only.")
+{
+ int nargin = args.length ();
+ octave_value_list retval;
+
+ if (nargin != 6)
+ {
+ print_usage ();
+ }
+ else
+ {
+ // arguments in
+ Matrix rfrdat = args(0).matrix_value ();
+ Matrix ifrdat = args(1).matrix_value ();
+ Matrix omega = args(2).matrix_value ();
+ int n = args(3).int_value ();
+ int discfl = args(4).int_value ();
+ int flag = args(5).int_value ();
+
+ int lendat = omega.rows (); // number of frequencies
+
+ int lda = max (1, n);
+
+ // arguments out
+ Matrix a (lda, n);
+ Matrix b (n, 1);
+ Matrix c (1, n);
+ Matrix d (1, 1);
+
+ // workspace
+ int liwork = max (2, 2*n + 1);
+ int ldwork;
+ int lzwork;
+ int hnpts = 2048;
+
+ int lw1 = 2*lendat + 4*hnpts;
+ int lw2 = lendat + 6*hnpts;
+ int mn = min (2*lendat, 2*n+1);
+ int lw3;
+ int lw4;
+
+ if (n > 0)
+ {
+ lw3 = 2*lendat*(2*n+1) + max (2*lendat, 2*n+1)
+ + max (mn + 6*n + 4, 2*mn + 1);
+ lzwork = lendat*(2*n+3);
+ }
+ else
+ {
+ lw3 = 4*lendat + 5;
+ lzwork = lendat;
+ }
+
+ if (flag == 1)
+ lw4 = max (n*n + 5*n, 6*n + 1 + min (1, n));
+ else
+ lw4 = 0;
+
+ ldwork = max (2, lw1, lw2, lw3, lw4);
+
+ OCTAVE_LOCAL_BUFFER (int, iwork, liwork);
+ OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
+ OCTAVE_LOCAL_BUFFER (Complex, zwork, lzwork);
+
+ // tolerance
+ double tol = 0;
+
+ // error indicator
+ int info;
+
+
+ // SLICOT routine SB10YD
+ F77_XFCN (sb10yd, SB10YD,
+ (discfl, flag,
+ lendat,
+ rfrdat.fortran_vec (), ifrdat.fortran_vec (),
+ omega.fortran_vec (),
+ n,
+ a.fortran_vec (), lda,
+ b.fortran_vec (),
+ c.fortran_vec (),
+ d.fortran_vec (),
+ tol,
+ iwork, dwork, ldwork,
+ zwork, lzwork,
+ info));
+
+ if (f77_exception_encountered)
+ error ("fitfrd: slsb10yd: exception in SLICOT subroutine SB10YD");
+
+ if (info != 0)
+ error ("fitfrd: slsb10yd: SB10YD returned info = %d", info);
+
+ // return values
+ retval(0) = a;
+ retval(1) = b;
+ retval(2) = c;
+ retval(3) = d;
+ retval(4) = octave_value (n);
+ }
+
+ return retval;
+}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
