Mercurial > hg > octave-terminal
changeset 561:e79ff1f4df3c
[project @ 1994-07-25 22:32:08 by jwe]
Initial revision
author | jwe |
---|---|
date | Mon, 25 Jul 1994 22:32:08 +0000 |
parents | 34c713db72b6 |
children | 50bf8d1b024e |
files | scripts/general/postpad.m scripts/polynomial/compan.m scripts/polynomial/polyderiv.m scripts/polynomial/polyreduce.m scripts/polynomial/polyval.m scripts/polynomial/polyvalm.m scripts/signal/filter.m |
diffstat | 7 files changed, 368 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/scripts/general/postpad.m @@ -0,0 +1,42 @@ +function y = postpad(x,l,c) +#postpad(x,l) +#Appends zeros to the vector x until it is of length l. +#postpad(x,l,c) appends the constant c instead of zero. +# +#If length(x) > l, elements from the end of x are removed +#until a vector of length l is obtained. + +# Author: +# Tony Richardson +# amr@mpl.ucsd.edu +# June 1994 + + if(nargin == 2) + c = 0; + elseif(nargin<2 || nargin>3) + error("usage: postpad(x,l) or postpad(x,l,c)"); + endif + + if(is_matrix(x)) + error("first argument must be a vector"); + elseif(!is_scalar(l)) + error("second argument must be a scaler"); + endif + + if(l<0) + error("second argument must be non-negative"); + endif + + lx = length(x); + + if(lx >= l) + y = x(1:l); + else + if(rows(x)>1) + y = [ x; c*ones(l-lx,1) ]; + else + y = [ x c*ones(1,l-lx) ]; + endif + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/compan.m @@ -0,0 +1,48 @@ +function A = compan(c) +#compan (c) +#Compute the companion matrix corresponding to polynomial vector c. +# +#In octave a polynomial is represented by it's coefficients (arranged +#in descending order). For example a vector c of length n+1 corresponds +#to the following nth order polynomial +# +# p(x) = c(1) x^n + ... + c(n) x + c(n+1). +# +#The corresponding companion matrix is +# _ _ +# | -c(2)/c(1) -c(3)/c(1) ... -c(n)/c(1) -c(n+1)/c(1) | +# | 1 0 ... 0 0 | +# | 0 1 ... 0 0 | +# A = | . . . . . | +# | . . . . . | +# | . . . . . | +# |_ 0 0 ... 1 0 _| +# +#The eigenvalues of the companion matrix are equal to the roots of the +#polynomial. +# +#SEE ALSO: poly, roots, residue, conv, deconv, polyval, polyderiv, polyinteg + +# Author: +# Tony Richardson +# amr@mpl.ucsd.edu +# June 1994 + + if(nargin != 1) + error("usage: compan(vector)"); + endif + + if(!is_vector(c)) + error("compan: expecting a vector argument."); + endif + + # Ensure that c is a row vector. + if(rows(c) > 1) + c = c.'; + endif + + n = length(c); + A = diag(ones(n-2,1),-1); + A(1,:) = -c(2:n)/c(1); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/polyderiv.m @@ -0,0 +1,33 @@ +function p = polyderiv(p) +#polyderiv(c) +#Returns the coefficients of the derivative of the polynomial whose +#coefficients are given by vector c. +# +#SEE ALSO: poly, polyinteg, polyreduce, roots, conv, deconv, residue, +# filter, polyval, polyvalm + +# Author: +# Tony Richardson +# amr@mpl.ucsd.edu +# June 1994 + + if(nargin != 1) + error("usage: polyderiv(vector)"); + endif + + if(is_matrix(p)) + error("argument must be a vector"); + endif + + lp = length(p); + if(lp == 1) + p = 0; + return; + elseif (lp == 0) + p = []; + return; + end + + p = p(1:(lp-1)) .* [(lp-1):-1:1]; + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/polyreduce.m @@ -0,0 +1,29 @@ +function p = polyreduce(p) +#polyreduce(c) +#Reduces a polynomial coefficient vector to a minimum number of terms, +#i.e. it strips off any leading zeros. +# +#SEE ALSO: poly, roots, conv, deconv, residue, filter, polyval, polyvalm, +# polyderiv, polyinteg + +# Author: +# Tony Richardson +# amr@mpl.ucsd.edu +# June 1994 + + index = find(p==0); + + index = find(index == 1:length(index)); + + if (length(index) == 0) + return; + endif + + if(length(p)>1) + p = p(index(length(index))+1:length(p)); + endif + + if(length(p)==0) + p = 0; + endif +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/polyval.m @@ -0,0 +1,41 @@ +function y = polyval(c,x) +#Evaluate a polynomial. +# +#In octave, a polynomial is represented by it's coefficients (arranged +#in descending order). For example a vector c of length n+1 corresponds +#to the following nth order polynomial +# +# p(x) = c(1) x^n + ... + c(n) x + c(n+1). +# +#polyval(c,x) will evaluate the polynomial at the specified value of x. +# +#If x is a vector or matrix, the polynomial is evaluated at each of the +#elements of x. +# +#SEE ALSO: polyvalm, poly, roots, conv, deconv, residue, filter, +# polyderiv, polyinteg + +# Author: +# Tony Richardson +# amr@mpl.ucsd.edu +# June 1994 + + if(nargin != 2) + error("usage: polyval(c,x)"); + endif + + if(is_matrix(c)) + error("poly: first argument must be a vector."); + endif + + if(length(c) == 0) + y = c; + return; + endif + + n = length(c); + y = c(1)*ones(rows(x),columns(x)); + for index = 2:n + y = c(index) + x .* y; + endfor +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/polyvalm.m @@ -0,0 +1,40 @@ +function y = polyvalm(c,x) +#Evaluate a polynomial in the matrix sense. +# +#In octave, a polynomial is represented by it's coefficients (arranged +#in descending order). For example a vector c of length n+1 corresponds +#to the following nth order polynomial +# +# p(x) = c(1) x^n + ... + c(n) x + c(n+1). +# +#polyvalm(c,X) will evaluate the polynomial in the matrix sense, i.e. matrix +#multiplication is used instead of element by element multiplication as is +#used in polyval. +# +#X must be a square matrix. +# +#SEE ALSO: polyval, poly, roots, conv, deconv, residue, filter, +# polyderiv, polyinteg + +# Author: +# Tony Richardson +# amr@mpl.ucsd.edu +# June 1994 + + if(nargin != 2) + error("usage: polyvalm(c,x)"); + endif + + if(is_matrix(c)) + error("poly: first argument must be a vector."); + endif + + if(!is_square(x)) + error("poly: second argument must be a square matrix."); + endif + + [v, d] = eig(x); + + y = v * diag(polyval(c,diag(d))) * v'; + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/signal/filter.m @@ -0,0 +1,135 @@ +function [y, w] = filter(b,a,x,w) +#Filter a vector. +#y = filter(b,a,x) returns the solution to the following linear, +#time-invariant difference equation: +# +# N M +# sum a(k+1) y(n-k) + sum b(k+1) x(n-k) = 0 for 1<=n<=length(x) +# k=0 k=0 +# +#where N=length(a)-1 and M=length(b)-1. An equivalent form of this +#equation is: +# +# N M +# y(n) = sum c(k+1) y(n-k) + sum d(k+1) x(n-k) for 1<=n<=length(x) +# k=1 k=0 +# +#where c = a/a(1) and d = b/a(1). +# +#In terms of the z-transform, y is the result of passing the discrete- +#time signal x through a system characterized by the following rational +#system function: +# +# M +# sum d(k+1) z^(-k) +# k=0 +# H(z) = ---------------------- +# N +# 1 + sum c(k+1) z(-k) +# k=1 +# +#[y, sf] = filter(b,a,x,si) sets the initial state of the system, si, +#and returns the final state, sf. The state vector is a column vector +#whose length is equal to the length of the longest coefficient vector +#minus one. If si is not set, the initial state vector is set to all +#zeros. +# +#The particular algorithm employed is known as a transposed Direct +#Form II implementation. +# +#SEE ALSO: poly, roots, conv, deconv, residue, polyval, polyderiv, polyinteg + +# Author: +# Tony Richardson +# amr@mpl.ucsd.edu +# June 1994 + + if(nargin < 3 || nargin > 4) + error("usage: [y, sf] = filter(b,a,x[,si])"); + endif + + if(is_matrix(a) || is_matrix(b) || is_matrix(x)) + error("Argument must be a vector."); + endif + + N = length(a); + M = length(b); + L = length(x); + + MN = max([N M]); + lw = MN - 1; + + # Ensure that all vectors have the assumed dimension. + if(columns(a) > 1) + a = reshape(a,N,1); + endif + if(columns(b) > 1) + b = reshape(b,M,1); + endif + + if(nargin == 3) + # Set initial state to zero. + w = zeros(lw,1); + else + if(is_matrix(w) || length(w) != lw) + error("state vector has the wrong dimensions."); + endif + if(columns(w) > 1) + w = reshape(w,lw,1); + endif + endif + + # Allocate space for result. + y = zeros(1,L); + + # It's convenient to pad the coefficient vectors to the same length. + b = postpad(b,MN); + + norm = a(1); + if (norm == 0.) + error("First element in second argument must be non-zero."); + endif + + if (norm != 1.) + b = b/norm; + endif + + # Distinguish between IIR and FIR cases. The IIR code can easily be made + # to work for both cases, but the FIR code is slightly faster when it can + # be used. + + if (N > 1) + # IIR filter. + a = postpad(a,MN); + if (norm != 1.) + a = a/norm; + endif + for index = 1:L + y(index) = w(1) + b(1)*x(index); + # Update state vector + if(lw > 1) + w(1:(lw-1)) = w(2:lw) - a(2:lw)*y(index) + b(2:lw)*x(index); + w(lw) = b(MN)*x(index) - a(MN) * y(index); + else + w(1) = b(MN)*x(index) - a(MN) * y(index); + endif + endfor + else + # FIR filter. + if(lw > 0) + for index = 1:L + y(index) = w(1) + b(1)*x(index); + if ( lw > 1) + # Update state vector + w(1:lw-1) = w(2:lw) + b(2:lw)*x(index); + w(lw) = b(MN)*x(index); + else + w(1) = b(2)*x(index); + endif + endfor + else + # Handle special case where there is no delay separately. + y = b(1)*x; + endif + endif +endfunction