MATLAB for Math of Signals/Systems Book

The directories here contain code associated with "Mathematical Methods and Algorithms for Signal Processing" by Todd K. Moon and Wynn C. Stirling.

DISCLAIMER: We make no warranty, express or implied, that the programs or algorithms presented in this book (or its accompanying media) are free of error, or that they will meet your requirements for any particular applications. They should not be relied upon for solving a problem whose incorrect solution could result in injury to a person or loss of property. Any and all use of the programs and algorithms associated with this book is at your own risk. The authors and publisher disclaim all liability for direct or consequential damages resulting from your use of the programs.

These m-files are user-contributed routines, and are not a product of The MathWorks, Inc. The MathWorks assume no responsibility for any errors that may exist in these routines.

USAGE: You are free to use the programs or any derivative of them for any scientific purpose, but please reference this book.

VERSION: The code listed here has been written with Matlab v. 5 or above. (Most of the algorithms will run with v. 4, but in some cases we make use of 'end', cell arrays, and other version 5 features.

DIRECTORIES:

bookalgs -- contains the algorithms listed or referenced in the book, as well as function which these reference
mkpict -- contains functions and Matlabs scripts that are used to make the Matlab-based plots appearing in the book. These are provided for researchers who want to profit from the form of a plot for their own work, and who want to test the result.

Correspondance between the figures and the Matlab files is provided by the file mkpict/listoffigures.

Many of these functions and scripts are were used as the basis for testing the algorithms referenced in the book, and in that regard should be regarded as part of the testing suite.
testers -- contains functions and scripts to test the algorithms appearing in the book, other than those that appear in the mkpict directory. These can be used as examples of calling sequences.
solutions -- most (but not all) homework problems requiring Matlab code to be written have code provided here. DON'T CHEAT yourself by rushing to the solution. Learn to do it yourself!
misc -- is a treasure trove of functions which were written more or less in parallel with the book. Some categories of functions: number theoretical methods: including fast convolution, CRT for integers and polynomials, Euclidean algorithm, polynomial operations, linear congruences, and many others. interpolation and approximation: including Bernstein polynomials, Neville's algorithm, splines, rational approximation, remez (minimax), Vandermonde matrix equation solver, trigonometric interpolation, pade approximaation, etc. examples related to the book: including neural networks, Viterbi, Pisarenko, Prony, composite mappings, HMM training (with unscaled forward/backward probabilities), etc.

The source of many of these is a part of the book --- which was removed due to lack of space --- about-number theoretical methods in signal processing, and another part --- also removed --- on interpolation and approximation.
data -- example data used in the blind source separation; the image for tomographic reconstruction; wavelet and scaling function data for plots.

TOOLBOXES:

The code written should run with Matlab without additional toolboxes, with the following known exceptions:

bookalgs/eigmakePQ.m requires the use of the sinc function (found in the signal processing toolbox)

bookalgs/testirwls.m makes use of the freqz function for plotting (found in the signal processing toolbox)

misc/kaisfilt.m requires the use of the sinc and the freqz function (found in the signal processing toolbox)

misc/ratintfilt.m makes use of the zplane function for plotting (found in the signal processing toolbox)

mkpict/testeigfil.m and mkpict/testeigfil2.m use the freqz function for plotting (found in the signal processing toolbox)

testers/bliter.m requires the use of the sinc function (found in the signal processing toolbox)

testers/lpfilttest.m makes use of the freqz function for plotting (found in the signal processing toolbox)

testers/testbplp.m makes use of the freqz function for plotting (found in the signal processing toolbox)

Directory: bookalgs

art1.m
% 
% Produce an updated solution x to Ax = b using the algebraic reconstruction
% technique.
% (This function is for demonstration, and does not precompute
% row norms as a more efficient version would.)  No constraints are applied.
%
% function x = art1(A,x,b)
%
% A = input matrix
% x = initial solution
% b = right-hand side
% 
% Output x= updated solution

bayes3.m
% Minimax decision for Gaussian

bidiag.m
% 
% bidiagonalize the real symmetric matrix A: UAV' = B,
% 
% function [B,U,V] = bidiag(A)
%
% A = matrix to be factored
%
% B = bidiagonal
% U = orthogonal transformation (optional output)
% V = orthogonal transformation (optional output)

bliter1.m
% Iterate on bandlimited data using two projections

bss1.m
% Blind source separation example

bssg.m
%
% The nonlinear function in the neural network
%
% function y = bssg(v).

cholesky.m
%
% Compute the Cholesky factorization of B, B = LL'
% (this version does not require additional storage)
%
% function [L] = cholesky(B)
%
% B = matrix to be factored
%
% L = lower triangular factor

compmap2.m
% Run the composite mapping algorithm to make a positive sequence

conjgrad1.m
% function [x,D] = conjgrad1(Q,b)
%
% Solve the equation Qx = b using conjugate gradient, where Q is symmetric
%
% Q = symmetric matrix
% b = right-hand side

conjgrad2.m
%
% Apply the conjugate gradient to minimize a function
%
% function [x,X] = conjgrad2(x,grad,hess)
%
% x = starting point
% grad = the gradient of the function to minimize (a function name)
% hess = the gradient of the function to minimize (a function name)
%
% Output x = update of point
% X = array of points examined (optional)

convnorm.m
%
% Compute the Hamming distance between the branchweights and the input
% This function may be feval'ed for use with the Viterbi algorithm
% (state and nextstate are not used here)
% 
% function d = convnorm(branch,input,state,nextstate)
%

dijkstra.m
% 
% Find the set of shortest paths from vertex a to each of the other vertices
% in the graph represented by the cost matrix.
% If there is no branch from i to j, then cost(i,j) = inf.
%
% function [dist,prevnode] = dijkstra(a,cost)
% 
% a = starting vectex
% cost = cost matrix for graph
%
% dist(i) = shortest distance from vertex a to vertex i.
% prevnode(i) = vertex prior to vertex i on the shortest path to i.

dirtorefl.m
% 
% Convert from direct-form FIR filter coefficients in a
% to lattice form
%
% function kappa = refltodir(a)
%
% a = direct form coefficients
%
% kapp = lattice form coefficients

durbin.m
% 
% solve the nxn Toeplitz system Tx = [r(1)..r(n)]
% given a vector r = (r_0,r_1,\ldots,r_{n+1}),
%
% Since matlab has no zero-based indexing, r(1) = r_0
%
% function [x] = durbin(r)
%
% r = input vector
% x = solution to Tx = r

eigfil.m
%
% Design an eigenfilter
% 
% function h = eigfil(wp,ws,m,alpha)
%
% wp = passband frequency
% ws = stopband frequency
% m = number of coefficients (must be ODD)
% alpha = tradeoff parameter between stopband and passband
%
% h = filter coefficients

eigfilcon.m
%
% Design an eigenfilter with values constrained at some frequencies
%
% function h = eigfilcon(wp,ws,N,alpha, Wo, d)
%
% wp = passband frequency
% ws = stopband frequency
% m = number of coefficients (must be ODD)
% alpha = tradeoff parameter between stopband and passband
% Wo = list of constraint frequencies
% d = desired magnitude values at those frequencies
%
% h = filter coefficients

eigmakePQ.m
%
% Make the P and Q matrices for eigenfiltering
%
% function [P,Q] = eigmakePQ(wp,ws,N)
%
% wp = passband frequency
% ws = stopband frequency
% N = number of coefficients

eigqrshiftstep.m
% 
% Perform the implicit QR shift on T, where the shift
% is obtained by the eigenvalue of the lower-right 2x2 submatrix of T
%
% function [T,Q] = eigqrshiftstep(T,Qin)
% 
% T = tridiagonal matrix
% Qin = (optional) orthogonal input matrix
%
% T = T matrix after QR shift
% Q = (optional) orthogonal matrix

em1.m
% Illustration of example em algorithm computations
%

esprit.m
%
% Compute the frequency parameters using the ESPRIT method
%
%function f = esprit(Ryy,Ryz,p)
%
%
% Ryy = (estimate of) the autocorrelation of the observations
% Ryz = (estimate of) the correlation between y[n] and y[n+1]
% p = number of modes
%
% f = vector of frequencies (in Hz/sample)

et1.m
% et1.m
% Perform tomographic reconstruction using the EM algorithm
% im is a (nr x nc) image scaled so that 0=white, 255=black
% nr = number of rows;   nc = number of columns

fblp.m
%
% Forward-backward linear predictor
%
% function [a,sigma] = fblp(x,n)
%
% Given a sequence of data in the columns vector x, 
% return the coefficients of an nth order forward-backward linear predictor in a
%
% x = data sequence
% n = length of filter
%
% a = filter coefficients
% sigma = standard deviation of noise

fordyn.m
%
% Do forward dynamic programing
%
% function [pathlist,cost] = fordyn(G,W)
%
% G = list of next vertices
% W = list of costs to next vertices
%
% pathlist = list of paths through graphs
% cost = cost of the paths

gamble.m
%
% Return the bets for a race with win probabilities p and subfair odds o,
%
% function [b0,B] = gamble(p,o)
%
% p = probability of win
% o = subfair odds

gaussseid.m
% 
% Produce an updated solution x to Ax = b using Gauss-Seidel iteration
%
% function x = gausseid(A,x,b)
%
% A = input matrix
% x = initial solution
% b = right-hand side
% 
% Output x= updated solution

golubkahanstep.m
% 
% Given a bidiagonal matrix B with NO zeros on the diagonal or
% superdiagonal, create a new B <-- U'BV, using an implicit QR shift on
% T = B'B
%
% function [B,U,V] = golubkahanstep(B,Uin,Vin)
%
% B = bidiagonal matrix
% Uin, Vin = last estimate of U and V
%
% B = new bidiagonal matrix
% U, V = new estimate of U and V

gramschmidt1.m
% 
% Compute the Gram-Schmidt orthogonalization of the 
% columns of A, assuming nonzero columns of A
%
% function [Q,R] = gramschmidt1(A)
%
% A = input matrix to be factored (assuming nonzero columns)
%
% Q = orthogonal matrix
% R = upper triangular matrix such that A = QR

hmmApiupn.m
% 
% update the HMM probabilities A and pi using the normalized forward and
% backward probabilities alphahat and betahat
%
% function [A,pi] = hmmapiupn(y,alphahat,betahat,HMM)
%
% y = input sequence
% alphahat = normalized alphas
% betahat = normalized betas
% HMM = current model parameters
%
% A = updated transition probability matrix
% pi = updated initial probability matrix

hmmabn.m
% 
% compute the normalized forward and backward probabilities for the model HMM
% and the output probabilities and the normalization factor
%
% function [alphahat,betahat,f,c] = hmmabn(y,HMM)
%
% y = input sequence y[1],y[2], ..., y[T]
% HMM = HMM data
%
% alphahat = [ alphahat(:,1) alphahat(:,2) ... alphahat(:,T)]
%            (alphahat(:,t) = alphahat(y_t^T,:))
% betahat = [betahat(:,2) ... betahat(:,T-1) betahat(:,T) betahat(:,T+1)]
%            (betahat(:,t) = betahat(y_{t+1}^T|:))
% f = initial probability types
% c = normalizing factors

hmmdiscf.m
%
% Compute the pmf value for a discrete distribution
%
% function f = hmmdiscf(f)
% 
% y = output value
% f = output distribution
% s = state

hmmdiscfupn.m
% 
% Update the discrete HMM output distribution f using the normalized forward
% and backward probabilities alphahat and betahat
% 
% function f = hmmdiscfupn(y,alphahat,betahat,c,HMM)
%
% y = output sequence
% alphahat = normalized alphas
% betahat = normalized betas
% c = normalization constants
% HMM = current model parameters
%
% f = updated output distribution

hmmf.m
% 
% Determine the likelihood of the output y for the model HMM
% This function acts as a clearinghouse for different probability types
%
% function p = hmmf(y,f,s)
%
% y = output value
% f = probability distribution
% s = state
%
% p = likelihood

hmmfupdaten.m
% 
% Provide an update to the state output distributions for the HMM model
% using the normalized probabilities alphahat and betahat
%
% function f = hmmfupdaten(y,alphahat,betahat,HMM)
%
% y = output sequence
% alphahat = normalized alphas
% betahat = normalized betas
% c = normalization factors
% HMM = current model paramters
%
% f = updated output distribution

hmmgausf.m
% 
% Compute the pmf value for a Gaussian distribution
%
% function f = hmmgausf(y,f,s)
%
% y = output value
% f = output distribution
% s = state

hmmgausfupn.m
% 
% Update the Gaussian output distribution f of the HMM using the normalized
% probabilities alphahat and betahat
%
% function f = hmmgausfupn(y,alphahat,betahat,c,HMM)
%
% y = output sequence
% alphahat = normalized alphas
% betahat = normalized betas
% c = normalization constants
% HMM = current model parameters
%
% f = updated output distribution

hmmgendat.m
% 
% Generate T outputs of a Hidden Markov model HMM
%
% function [y,ss] = hmmgendat(T,HMM)
%
% T = number of samples to produce
% HMM = HMM model parameters
%
% y = output sequence
% s = (optional) state sequence

hmmgendisc.m
% 
% Generate T outputs of a HMM with a discrete output distribution
%
% function y = hmmgendisc(T,HMM)
% 
% T = number of outputs to produce
% HMM = model parameters
%
% y = output sequence
% ss = (optional) state sequence (for testing purposes)

hmmgengaus.m
% 
% Generate T outputs of a HMM with a Gaussian output distribution
%
% function y = hmmgengaus(T,HMM)
%
% T = number of outputs to produce
% HMM = model parameters
%
% y = output sequence
% ss = (optional) state sequence (for testing purposes)

hmminitvit.m
% 
% Initialize the Viterbi algorithm stuff for HMM sequence identification
%
% function hmminitvit(inHMM,inpathlen)
% 
% inHMM = a structure containing the initial probabilities, state transition
%         probabilities, and output probabilities
% inpathlen = length of window used in VA

hmmlpyseqn.m
%
% Find the log likelihood of the sequence y[1],y[2],...,y[T],
% using the parameters in HMM.% That is, compute log P(y|HMM),
%
% function lpy = hmmlpyseqn(y,HMM)
%
% y = input sequence = y[1],y[2],\ldots,y[T]
% HMM = current estimate of HMM parameters
%
% lpy = log P(y|HMM)

hmmlpyseqv.m
% 
% Find the log likelihood of the sequence y[1],y[2],...,y[T]
% i.e., compute log P(y|HMM),
% using the the Viterbi algorithm
%
% function lpy = hmmlpyseqv(y,HMM)
%
% y = input sequence
% HMM = HMM parameters
%
% lpy = log likelihood value

hmmnorm.m
% 
% Compute the branch norm for the HMM using the Viterbi approach
%
% function d = hmmnorm(branchweight,y,state,nextstate)
%
% branchweight= log transition probability 
% y = output
% state = current state in trellis
% nextstate = next state in trellis
%
% d = branch norm (log-likelihood)

hmmnotes.m
% Notes on data structures and functions for the HMM
%

hmmtest2vb.m
% Test the HMM using both Viterbi and EM-algorithm based training methods

hmmupdaten.m
%
% Compute updated HMM model from observations
%
% function HMM = hmmupdaten(y,HMM) 
%
% y = output sequence
% HMM = current model parameters
%
% hmmo = updated model parameters

hmmupdatev.m
% 
% Compute updated HMM model from observations y using Viterbi methods
% Assumes only a single observation sequence.
%
% function HMM = hmmupdatev(y,HMM) 
%
% y = sequence of observations
% HMM = old HMM (to be updated)
%
% hmmo = updated HMM

hmmupfv.m
% 
% Compute an update to the distribution f based upon the data y
% and the (assumed) state assignment in statelist
%
% function fnew = hmmupfv(y,statelist,n,f)
%
% y = sequence of observations
% statelist = state assignments
% n = number of states
% f = distribution (cell) to update
%
% fnew = updated distribution

houseleft.m
%
% Apply the Householder transformation based on v to A on the left
%
% function A = houseleft(A,v)
%
% A = an mxn matrix
% v = a household vector
%
% B = H_v A

houseright.m
%
% Apply the householder transformation based on v to A on the right
%
% function A = houseright(A,v)
%
% A = an mxn matrix
% v = a household vector
%
% B = H_v A

ifs3a.m
% Plot the logistic map and the orbit of a point
%

initcluster.m
%
% 
% Choose an initial cluster at random
% 
% function Y = initcluster(X,m)
%
% X = input data: each column is a training data vector
% m = number of clusters
% Y = initial cluster: each column is a point

initpoisson.m
% 
% Initialize the global variables for the poisson generator
% 
% function initpoisson

initvit1.m
% 
% Initialize the data structures and pointers for the Viterbi algorithm
% 
% function initvit1(intrellis, inbranchweight, inpathlen, innormfunc)
%
% intrellis: a description of the successor nodes 
%    e.g. [1 3; 3 4; 1 2; 3 4]
% inbranchweight: weights of branches used for comparison, saved as
%    cells in branchweight{state_number, branch_number}
%    branchweight may be a vector
%    e.g.  branchweight{1,1} = 0; branchweight{1,2} = 6;
%          branchweight{2,1} = 3; branchweight{2,2} = 3;
%          branchweight{3,1} = 6; branchweight{3,2} = 0;
%          branchweight{4,1} = 3; branchweight{4,2} = 3;
% inpathlen: length of window over which to compute
% normfun: the norm function used to compute the branch cost

invwavetrans.m
%
% Compute the inverse discrete wavelet transform 
%
% function c = invwavetrans(C,ap,coeff)
%
% C = input data (whose inverse transform is to be found)
% ap = index for start of coefficients for the jth level
% coeff = wavelet coefficients
%
% c = inverse transformed data

invwavetransper.m
%
% Compute the periodized inverse discrete wavelet transform 
%
% function c = invwavetransper(C,coeff,J)
%
% C = input data
% coeff = wavelet coefficients
% J = (optional) number of levels of inverse transform to compute
%    If length(C) is not a power of 2, J must be specified.
%
% c = inverse discrete wavelet transform of C

irwls.m
% 
% Computes the minimum solution c to ||x-Ac||_p using
% iteratively reweighted least squares
%
% function c = irwls(A,x)
%
% A = system matrix
% x = rhs of equation
% p = L_p norm
%
% c = solution vector

jacobi.m
% 
% Produce an updated solution x to Ax = b using Jacobi iteration
%
% function x = jacobi(A,x,b)
%
% A = input matrix
% x = initial solution
% b = right-hand side
% 
% Output x= updated solution

kalex1.m
% Kalman filter example 1
%

kalman1.m
% 
% Computes the Kalman filter esimate xhat(t+1|t+1)
% for the system x(t+1) = Ax(t) + w
%                y(t) = Cx(t) + v
% where cov(w) = Q  and cov(v) = R, 
% The prior estimate is x0, and the prior covariance is P0.
%

karf.m
%
% Evaluate the potential function f(x,c)
% for karmarkers algorithm
%
% function f = karf(x,c)
% 
% x = value of x
% c = constraint vector
%
% f = potential function

karmarker.m
% 
% Implement a Karmarker-type algorithm for linear programming
% to solve a problem in "Karmarker standard form"
%  min       c'x
% subject to Ax=0,  sum(x)=1, x >=0
%
% function x = karmarker(A,c)
%
% A,c = system matrices
%
% x = solution

kissarma.m
% 
% Determine the ARMA parameters a and b of order p based upon the data in y.
%
% function [a,b] = kissarma(y,p)
%
% y = sequence 
% p = order of AR part
%
% a = AR coefficients
% b = MA coefficients

levinson.m
% 
% Given a vector r = (r_0,r_1,\ldots,r_{n-1}),
% and a vector b = (b_1,b_2,\ldots,b_n)
% solve the nxn Toeplitz system Tx = b
%
% function [y] = levinson(r,b)
%
% Since Matlab has no zero-based indexing, r(1) = r_0
%
% r = vector of coefficients for Toeplitz matrix
% b = right-hand side
% 
% y = solution to Tx = b

lgb.m
% 
% Find m clusters on the data X
%
% function [Y,d] = lgb(X,m)
%
%
% X = input data: each column is a training data vector
% m = number of clusters to find
%
% Y = set of clusters: each column is a cluster centroid
% d = minimum total distortion

lms.m
% 
% Given a (real) scalar input signal x and a desired scalar signal d,
% compute an LMS update of the weight vector h.
% This function must be initialized by lmsinit
%
% function [h,eap] = lms(x,d)
%
% x = input signal (scalar)
% d = desired signal (scalar)
%
% h = updated LMS filter coefficient vector
% eap = (optional) a-priori error

lmsinit.m
%
% Initialize the LMS filter
% 
% function lmsinit(m,mu)
%
% m = dimension of vector
% mu = lms stepsize

logistic.m
% 
% Compute the logistic function y = lambda*x*(1-x)
%
% function y = logistic(x,lambda)
%
% x = input value (may be a vector)
% lambda = factor of the function

lpfilt.m
% 
% Design an optimal linear-phase filter using linear programming
%
% function [h,delta] = lpfilt(fp,fs,n)
%
% fp = pass-band frequency  (0.5=Fs/2)
% fs = stop-band frequency
% n = number of coefficients (assumed odd here)
%
% h = filter coefficients
% delta1 = pass-band ripple

lsfilt.m
%
% Determine a least-squares filter h with m coefficients 
%
% function [h,X] = lsfilt(f,d,m,type)
%
% f = input data
% d = desired output data
% m = order of filter
% type = data matrix type
%     type=1: "covariance" method    2: "autocorrelation" method
%          3: prewindowing           4: postwindowing
%
% h = least-squares filter
% X = (optional) data matrix

makehankel.m
% 
% form a hankel matrix from the input data y
%
% [H] = makehankel(y,m,n)
%
% y = input data  = [y1 y2 ...] (a series of vectors in a _row_)
% m = number of block rows in H
% n = number of block columns in H
%
% H = Hankel matrix formed from y

makehouse.m
%
% Make the Householder vector v such that Hx has zeros in 
% all but the first component
%
% function v = makehouse(x)
%
% x = vector to be transformed
%
% v = Householder vector

massey.m
%
% Return the shortest binary (GF(2)) LFSR consistent with the data sequence y
%
% function [c] = massey(y)
%
% y = input sequence 
%
% c = LFSR connections, c = 1 + c(2)D + c(3)D^2 + ... c(L+1)D^L
%     (Note: opposite from usual Matlab order)

maxeig.m
%
% Compute the largest eigenvalue and associated eigenvector of 
% a matrix A using the power method
%
% function [lambda,x] = maxeig(A)
%
% A = matrix whose eigenvalue is sought
%
% lambda = largest eigenvalue
% x = corresponding eigenvector

mineig.m
% 
% Compute the smallest eigenvalue and associated eigenvector of 
% a matrix A using the power method
% function [lambda,x] = mineig(A)
%
% A = matrix whose eigenvalue is sought
%
% lambda = minimum eigenvalue
% x = corresponding eigenvector

musicfun.m
%
% Compute the "MUSIC spectrum" at a frequency f.
%
% function pf = musicfunc(f,p,V)
%
% f = frequency (may be an array of frequencies)
% p = order of system
% V = eigenvectors of autocorrelation matrix
%
% pf = plotting value for spectrum

neweig.m
% 
% Compute the eigenvalues and eigenvector of a real symmetric matrix A
%
% function [T,Q] = neweig(A)
%
% A = matrix whose eigendecomposition is sought
%
% T = diagonal matrix of eigenvalues
% Q = (optional) matrix of eigenvectors

newlu.m
%
% Compute the lu factorization of A
%
% function [lu,indx] = newlu(A)
%
% A = matrix to be factored
%
% lu = matrix containg L and U factors
% indx = index of pivot permutations

newsvd.m
% 
% Compute the singular value decomposition of the mxn matrix A, as A= u s v'.
% We assume here that m>n
% 
% [u,s,v] = newsvd(A)
% or
% s = newsvd(A)
%
% A = matrix to be factored
% 
% Output:
% s = singular values
% u,v = (optional) orthogonal matrices

nn1.m
%
% Compute the output of a neural network with weights in w
%
% function [y,V,Y] = nn1(xn,w)
% 
% xn = input
% w = cell array of weights
%
% y = output layer output
% V = (optional) internal activity
% Y = (optional) neuron output
%    The optional arguments V and Y are used for training to store output for
%    each layer:
%    Y{1} = input, Y{2} = first hidden layer, etc.
%    V{1} = first hidden layer, etc.

nnrandw.m
% 
% Generate an initial set of weights for a neural network at random,
% based upon the list in m
%
% function w = nnrandw(m)
% m = list of weights 
%   m(1) = number of inputs, m(2) = first hidden layer, etc.
%
% w = random weights

nntrain1.m
%
% Train a neural network using the input/output training data [x,d]
%
% function w = nntrain(x,d,m,ninput,mu)
%
% x = [x(1) x(2) ... x(N)] = input training data   
% d = [d(1) d(2) ... d(N)] = output training data
% nlayer = number of layers
% m = number of neurons on each layer, 
%     m(1) = input layer, ... m(nlayer+1) = ouput layer
% mu = steepest descent step size
% alpha = (optional) momentum constant
% maxiter = (optional) maximum number of iterations (w = no maximum)
% w = (optional) starting weights
%
% w = new weights
% err = (optional) total squared error from training

permutedata.m
% 
% Randomly permute the columns of the data x.
%
% function xp = permutedata(x,type)
%
% x = data to permute
% type=type of permutation
%   type=1: Choose a random starting point, and go sequentially
%   type=2: random selection without replacement (not really a permutation)
%
% xp = permuted x

pisarenko.m
%
% Compute the the modal frequencies using Pisarenko's method,
% then find the amplitudes
%
% function [f,P] = pisarenko(Ryy)
%
% Ryy = autocorrelation function from observed data
%
% f = vector of frequencies
% P = vector of amplitudes

pivottableau.m
% 
% Perform pivoting on an augmented tableau until 
% there are no negative entries on the last row
%
% function [tableau,basicptr] = pivottableau(intableau,inbasicptr)
%
% intableau = input tableau tableau,
% inbasicptr = a list of the basic variables, such as [1 3 4]
%
% tableau = pivoted tableau 
% basicptr = new list of basic variables

poisson.m
% 
% Generate a sample of a random variable x with mean lambda
% (Following Numerical Recipes in C, 2nd ed., p. 294)
% This function should be initialized by initpoisson.m
%
% function x = poisson(lambda)
%
% lambda = Poisson mean
%
% x = Poisson random variable

ptls1.m
% 
% Compute the Partial Total Least Squares solution of Ax = b
% where the first k columns of A are not modified
%
% function [x,Ahat,bhat] = ptls1(A,b,k)
%
% A = system matrix
% b = right-hand side
% k = number of columns of A not modified
%
% x = ptls solution to Ax=b
% Ahat = modified A matrix
% bhat = modified b matrix

ptls2.m
% 
% Find the partial total least-squares solution to Ax = b,
% where k1 rows and k2 columns of A are unmodified
% 
% function [x] = ptls2(A,b,k1,k2)
%
% A = system matrix
% b = right-hand side
% k1 = number of rows of A not modified
% k2 = number of columns of A not modified
%
% x = PTLS solution to Ax=b

qf.m
% 
% Compute the Q function:
%
% function p = qf(x)
%   p = 1/sqrt(2pi)int_x^infty exp(-t^2/2)dt

qfinv.m
% 
% Compute the inverse of the q function
%
% function x = qfinv(q)

qrgivens.m
% 
% Compute the QR factorization of a matrix A without column pivoting
% using Givens rotations
%
% function [R,thetac,thetas] = qrgivens(A)
% 
% A = mxn matrix (assumed to have full column rank)
%
% R = upper triangular matrix
% thetac = matrix of c values 
% thetas = matrix of s values (these must be converted to Q)

qrhouse.m
% 
% Compute the QR factorization of a matrix A without column pivoting
% using Householder transformations
% 
% function [V,R] = qrhouse(A)
%
% 
% A = mxn matrix (assumed to have full column rank)
%
% V = [v1 v2 ... vn] = matrix of Householder vectors 
%           (must be converted to Q ) using qrmakeq.
%           
% R = upper triangular matrix

qrmakeq.m
% 
% Convert the V matrix returned by qrhouse into a Q matrix
% 
% function Q = qrmakeq(V)
%
% V = [v1 v2 .... vr], Householder vectors produced by qrhouse
%
% Q = [Q1 Q2 ... Qr], the orthogonal matrix in the QR factorization

qrmakeqgiv.m
% 
% Given thetac and thetas containing rotation parameters from Givens rotations,
% (produced using qrqrgivens), compute Q
% function Q = qrmakeqgiv(thetac,thetas)
% 
% thetac = cosine component of Givens rotation
% thetas = sin component of Givens rotation
%
% Q = orthogonal matrix formed by Givens rotations

qrqtb.m
% 
% Given a matrix V containing Householder vectors as columns
% (produced using qrhouse), compute Q^H B, where Q is formed (implicitly)
% from the Householder vectors in V.
%
% function [B] = qrqtb(B,V)
% 
% B = matrix to be operated on
% V = matrix of Household vectors (as columns)
%
% output: B = Q^H B

qrqtbgiv.m
% 
% Given thetac and thetas containing rotation parameters from Givens rotations,
% (produced using qrqrgivens), compute Q^H B, where Q is formed (implicitly)
% from the rotations in Theta.
%
% function [B] = qrqtbgiv(B,thetac,thetas)
%
% B = matrix to be opererated on
% thetac = cosine component of rotations from Givens rotations
% thetas = sine component of rotations from Givens rotations
%
% Output: B <-- Q^H B

qrtheta.m
% 
% Given x and y, compute cos(theta) and sin(theta) for a Givens rotation
%
% function [c,s] = qrtheta(x,y)
%
% (x, y) = point to determine rotation
%
% c = cos(theta),   s=sin(theta)

reducefree.m
% 
% Perform elimination on the free variables in a linear programming problem
%
% function [A,b,c,value,savefree,nfree] = reducefree(A,b,c,freevars)
% 
% A,b,c = parameters from linear programming problem
% freevars = list of free variables
%
% A,b,c = new parameters for linear programming problem (with free variables
%         eliminated)
% value = value of linear program
% savefree = tableau information for restoring free variables
% nfree = number of free variables found

refltodir.m
% 
% Convert from a set of reflection coefficients kappa(1)...kappa(m)
% to FIR filter coefficients in direct form
%
% function a = refltodir(kappa)
%
% kappa = vector of reflection coefficients
%
% a = output filter coefficients = [1 a(1) a(2) ... a(m)]

restorefree.m
% 
% Restore the free variables by back substitution
%
% function x = restorefree(inx,savefree,freevars)
% 
% inx = linear programming solution (without free variables)
% savefree = information from tableau for substitution
% freevars = list of free variables
% 
% x = linear programming solution (including free variables)

rls.m
%
% Given a scalar input signal x and a desired scalar signal d,
% compute an RLS update of the weight vector h.
%
% This function must be initialized by rlsinit
%
% function [h,eap] = rls(x,d)
%
% x = input signal
% d = desired signal
%
% h = updated filter weight vector
% eap = (optional) a-priori estimation error

rlsinit.m
%
% Initialize the RLS filter
%
% function rlsinit(m,delta)

simplex1.m
% 
% Find the solution of a linear programming problem in standard form
%  minimize c'x
%  subject to Ax=b
%             x >= 0
%
% function [x,value,w] = simplex1(A,b,c,freevars)
% 
% A, b, c: system problem
% freevars = (optional) list of free variables in problem
%
% x = solution
% value = value of solution
% w = (optiona)l solution of the dual problem.
%    If w is used as a return value, then the dual problem is also solved.
%    (In this implementation, the dual problem cannot be solved when free
%    variables are employed.)

sor.m
% 
% Produce an updated solution x to Ax = b successive over-relaxation
% A must be Hermitian positive definite
%
% function x = sor(A,x,b)
%
% A = input matrix
% x = initial solution
% b = right-hand side
% omega = relaxation parameter
% 
% Output x= updated solution

sysidsvd.m
% 
% given a sequence of impulse responses in h
% identify a system (A,B,C)
%
% function [A,B,C] = sysidsvd(h,tol)
% h = impulse responses (a cell array)
% tol = (optional) tolerance used in rank determination
% ord = system order (overrides tolerance)
%
% (A,B,C) = estimated system matrix parameters

testet.m
% testet.m
% Test the emission tomography code
% This script loads an image file, plots it, then calls the 
% code to test the tomographic reconstruction

testirwls.m
% Test the irwls routine for a filter design problem
% After [Burrus 1994, p. 2934]

testnn10.m
% Test the neural network on a pattern recognition problem
%

tls.m
% 
% determine the total least-squares solution of Ax=b
%
% function x = tls(A,b)
%
% A = system matrix
% b = right-hand side
%
% x = TLS solution to Ax=b

tohankel.m
% 
% Determine the matrix nearest to A which is Hankel and has rank r
% using the composite mapping algorithm
%
% function A = tohankel(A,r)
%
% A = input/output matrix
% r = desired rank
%
% Ouptut A = Hankel matrix with desired properties
% d = norm of difference between matrices

tohanktoep.m
% 
% Determine the matrix nearest to A which is the stack
%  [ Hankel
%    Toeplitz ] 
% with rank r using the composite mapping algorithm
%
% function A = tohantoep(A,r)
%
% A = input matrix
% r = desired rank
%
% output: A = matrix with desired properties
% normlist = (optional) vector of errors

tokarmarker.m
% 
% Given a linear programming problem in standard form
%  min       c'x
% subject to Ax = b,  x >= 0
%
% return new values of A and c in "Karmarker standard form"
%  min       c'x
% subject to Ax=0,  sum(x)=n, x >= 0
%
% function [Anew,cnew] = tokarmarker(A,b,c)
%
% (A,b,c) = matrices in standard form
%
% Anew, cnew = matrices in Karmarker standard form

tostoch.m
% 
% Determine the matrix nearest to A which is stochastic using
% the composite mapping algorithm
%
% function A = tostoch(A)
% A = input matrix
%
% Output: A = nearest stochastic A

tridiag.m
% 
% tridiagonalize the real symmetric matrix A
% 
% function [T,Q] = tridiag(A)
%
% A = matrix whose tridiagonal form is sought
%
% T = tridiagonal matrix
% Q = (optional) orthogonal transformation

vitbestcost.m
% 
% Returns the best cost so far in the Viterbi algorithm
%
% function c = vitbestcost

viterbi1.m
% 
% Run the Viterbi algorithm on the input for one branch
% Before calling this function, you must call initvit1
%
% function [p] = viterbi1(r)
%
% r = input value (scalar or vector)
%
% p = 0 if number of branches in window is not enough
% p = statenum on path if enough branches in window
%
% Call vitflush to clear out path when finished

vitflush.m
% 
% Flush out the rest of the viterbi paths
%
% function [plist] = vitflush(termnode)
%
% termnode = 0 or list of allowable terminal nodes
%    If termnode==0, then choose path with lowest cost
%    Otherwise, choose path with best cost among termnode
%
% plist = list of paths from trellis

warp.m
% 
% find the dynamic warping between A and B (which may not be of the
% same length)
%
% function [path] = warp(A,B)
%
% A = cells of the vectors, A{1}, A{2}, ..., A{M}
% B = cells of the vectors, B{1}, B{2}, ..., B{N}
%
% path = Kx2 array of (i,j) correspondence

warshall.m
% 
% Find the transitive closure of the graph represented by the adjacency
% matrix A
%
% function Anew = warshall(A)
%
% A = adjacency matrix of graph
% 
% Anew = adjacency matrix for transitive closure of graph

wavecoeff.m
% Coefficients for Daubechies wavelets
%

wavetest.m
% Test the wavelet transform in matrix notation

wavetesto.m
% Test the wavelet transform in wavelet notation (alternate indexing)

wavetrans.m
%
% Compute the (nonperiodized) discrete wavelet traqnsform 
%
% function [C,ap] = wavetrans(c,coeff,J)
%
% c = data to be transformed
% coeff = wavelet transform coefficients
% J = number of levels of transform
%    If J is specified, then J levels of the transform are computed.  
%    Otherwise, the largest possible number of levels are used.
%
% C = transformed data
%     The output is stacked in C in wavelet-coefficient first order,
%     C = [d1 d2 ... dJ cJ]
% ap indexes the start of the coefficients for the jth level,
%     ap(j+1) indexes the start of the coefficients for the jth level,
%     except that ap(1) indicates the number of original datapoints
%
% This function simply stacks up the data for a call to the function wave

wavetransper.m
%
% Compute the periodized discrete wavelet transform of the data
%
% function [C] = wavetransper(c,coeff,J)
%
% c = data to be transformed
% coeff = transform coefficients
% J indicates how many levels of the transform to compute.
%    If length(c) is not a power of 2, J must be specified.
%
% C = output vector
%    The output is stacked in C in wavelet-coefficient first order,
%    C = [d1 d2 ... dJ cJ]

wftest.m
% Test the Wiener filter equalizer for a first-order signal and first-order
% channel

zerorow.m
% 
% Zero a row by a series of Givens rotations
%
% function [B,U] = zerorow(B,f,U)
%
% B = matrix to have row zeroed
% f = vector of row indices that are zero on the diagonal
% U = (optional) rotation matrix
%
% B = modified matrix
% U = (optional) rotation matrix

Directory: misc

H.m
% 
% Compute the binary entropy function
% 
% function h = H(p)
%
% p = crossover probability
%
% h = binary entropy

b2n.m
% convert an m-bit binary sequence b to an integer

bernapprox.m
% Plot Bernstein polynomials

bernpoly.m
% 
% compute the Bernstein polynomial g_{nk}(t)
%
% function g = bernpoly(n,k,t)
%
% n = degree
% k = order
% t = location
% 
% g = value

bsplineval1.m
% 
% When f(x) = sum_i c_i^k B_i^k(x)
% where B_i^k is the kth order spline across t_i,t_{i+1},\ldots,t_{i+k+1}
% This function evaluates f(x) for a given x.
% (See Kincaid-Cheney, p. 396)
%
% function s = bsplineval1(c,t,x,k)
%
% c = set of coefficients
% t = knot set
% x = value.  x must fall in range of knots, t_i <= x < t_{i+1}
% k = order
%
% s = spline value

chebinterp.m
% data for  Chebyshev interpolation example

chi2.m
% 
% Compute the pdf for a chi-squared random variable
%
% function f = chi2(x,n)
%
% x = value
% n = order of chi-squared
%
% f = pdf value

compmap1.m
% Test the composite mapping algorithm on the positive sequence

compmap4.m
% test the positive semi-definite mapping

computeremez.m
% Test the remez algorithm

crtgamma.m
% 
% Compute the gammas for the CRT
%
% function gamma = crtgamma(m)
%
% m = set of modulos
%
% gamm = set of gammas

crtgammapoly.m
% 
% Compute the gammas for the CRT for polynomials
%
% function gamma = crtgammapoly(m)
%
% m = set of modulos (polynomials)
%
% gamma = set of gammas (polynomials)

crypttest.m
% test the cryptographic example

d2b.m
% convert n to an m-bit binary representation

diagstack.m
% Stack matrices diagonally:
% D = [X 0
%      0 Y];
%
% function D = diagstack(X,Y)
% X, Y = input matrices
%
% D = diagonal stack

discapprox.m
% find a discrete approximating polynomial

divdiff.m
% 
% Compute the upper row of a divided difference table
%
% function c = divdiff(ts,fs)
%
% ts = sample locations
% fs = function values
%
% c = set of divided differences

eigfilcon0.m
%
% Find eigenilter constrained so that response is 0 at some frequencies
%
% function h = eigfilcon(wp,ws,N,alpha, Wo)
%
% wp = passband frequency
% ws = stopband frequency
% m = number of coefficients (must be ODD)
% alpha = tradeoff parameter between stopband and passband
% Wo = list of constraint frequencies at which response is 0
%
% h = filter coefficients

eigfilcon0new.m
%
% Find eigenilter constrained so that response is 0 at some frequencies
%
% function h = eigfilcon(wp,ws,N,alpha, Wo)
%
% wp = passband frequency
% ws = stopband frequency
% m = number of coefficients (must be ODD)
% alpha = tradeoff parameter between stopband and passband
% Wo = list of constraint frequencies at which response is 0
%
% h = filter coefficients

elem.m
%
% Return an elementary matrix E_{rs} of size mxn

fact.m
% compute the factorial

findprim5.m
% Find a primitive polynomial in GF(5)

fromcrt.m
% 
% Given a sequence [y1,y2,\ldots,y2] that is a representation 
% of an integer x in CRT form, convert back to x.
%
% function x = fromcrt(y,m)
%
% y = sequence
% m = [m1,m2,...,mr]
% gammain (optional) = set of gamma factors.  
%   If not passed in, gamma is computed
%
% x = integer representation
% gamma = gamma values

fromcrtpoly.m
% 
% Compute the representation of the polynomial f using the Chinese Remainder
% Theorem (CRT) using the moduli m = [m1,m2,...,mr].  It is assumed 
% (without checking) that the moduli are relatively prime.  
% Optionally, the gammas may be passed in (speeding computation), and
% are returned as optional return values.
%
% function [f,gamma] = fromcrtpoly(y,m,gammain)
% function [f] = fromcrtpoly(y,m)
% function [f,gamma] = fromcrtpoly(y,m)
% function [f] = fromcrtpoly(y,m,gammain)
%
% y = list of polynomials (cell array)
% m = list of moduli (cell array)
% gammain = (optional)list of gammas (cell array)
%
% f = reconstructed polynomial
% gamma = gamma

fromhankel.m
%
% Pull sequential data out of a Hankel matrix X 
%
% X = Hankel matrix
% d = (optional) dimension of data
%
% x = Sequential data

fromhankel2.m
% Pull sequential data out of a Hankel matrix X  (cell array version)
%
% X = Hankel matrix
% d = (optional) dimension of data
%
% x = Sequential data (in a cell array)

gam1.m
% 
% Determine the optimum bet on a subfair track
%
% function [B,b0,b,l] = gam1(p,o)
%
% p = probability of win
% o = subfair odds
%
% B = other bets
% b0 = amount witheld
% b = bet

gcdint1.m
% 
% Compute (only) the GCD (a,b) using the Euclidean algorithm
%
% function g = gcdint1(b,c)
%
% b,c = integers
%
% g = GCD(b,c)

gcdint2.m
% 
% Compute the GCD g = (b,c) using the Euclidean algorithm
% and return s,t such that bs+ct = g
%
% function [g,s,t] = gcdint2(b,c)
%
% b,c = integers
% g = GCD(b,c)
% s,t = integers

gcdpoly.m
% 
% Compute the GCD g = (b,c) using the Euclidean algorithm
% and return s,t such that bs+ct = g, where b and c are polynomials
% with real coefficients
%
% function [g,s,t] = gcdpoly(b,c)
%
% b,c = polynomials
% thresh = (optional) threhold argument used to truncate small remainders
%
% g = GCD(b,c)
% s,t = polynomials

genardat.m
% 
% Generate N points of AR data with a = [a(1) a(2), \ldots, a(n)]'
% and input variance sigma^2
%
% function x = genardat(a,sigma,N)
%
% a = AR parameters
% sigma = standard deviation of input noise
% N = number of points
%
% x = AR process

greedyperm.m
% 
% Using a greedy algorithm, determine a permutation P such that Px=z
% as closely as possible.
%
% function P = greedyperm(x,z)

greedyperm2.m
% 
% Using a greedy algorithm, determine a permutation P such that Px=z
% as closely as possible.
% This algorithm is more complex than greedyperm
%
% function P = greedyperm2(x,z)

haar1.m
% Do the computations for the Haar transform, working toward the
% wavelet lifting transform
%
% s = input data
% nlevel = number of levels
%
% h = Haar transform

haarinv1.m
% Do the computations for the inverse Haar transform, working toward the
% wavelet lifting transform

hmmApiup.m
% 
% Update the A and pi probabilities in the HMM using the forward and
% backward probabilities alpha and beta
%
% function [A,pi] = hmmapiup(y,alpha,beta,HMM)
%
% y = input sequence
% alpha = forward probability
% beta = backward probability
% f = distribution
% HMM = model parameters
%
% A = updated state transition probability
% pi = updated initial state probability

hmmab.m
% 
% Compute the forward and backward probabilities for the model HMM
% and the output probabilities

hmmdiscfup.m
% 
% Update the output probability distribution f of the HMM using the forward
% and backward probabilities alpha and beta
%
% function f = hmmdiscfup(y,alpha,beta,HMM)
%
% y = input sequence
% alpha = forward probabilities
% beta = backward probabilities
% HMM = current model parameters
%
% f = updated distribution

hmmfupdate.m
% 
% Provide an update to the state output distributions for the HMM model
%
% function f = hmmfupdate(y,alpha,beta,HMM)
% 
% y = input sequence
% alpha = forward probabilities
% beta = backward probabilities
% HMM = current model parameters
% 
% f = updated distribution

hmmgausfup.m
% 
% Update the Gaussian output distribution f of the HMM using the
% probabilities alpha and beta
%
% function f = hmmgausfup(y,alpha,beta,HMM)
% 
% y = input sequence
% alpha = forward probabilities
% beta = backward probabilities
% HMM = current model parameters
% 
% f = updated distribution

hmmlpyseq.m
%
% Find the log likelihood of the sequence y[1],y[2],...,y[T], 
% i.e., compute log P(y|HMM)
%
% function lpy = hmmlpyseq(y,HMM)
%
% y = input sequence
% HMM = current model parameters

hmmupdate.m
% 
% Compute updated HMM model from observations (nonnormalized)
%
% function HMM = hmmupdate(y,HMM) 
%
% y = output sequence
% HMM = current model parameters
%
% hmmo = updated model

ifs2.m
% test some ifs stuff

ifs2ex.m
% find an affine transformation Ax + b that transforms from
% {x00,x10,x20,x30} to {x01,x11,x21,x31}

initvit2.m
% 
% Initialize the data structures and pointers for the Viterbi algorithm
%
% function initvit2(intrellis, inbranchweight, inpathlen, innormfunc)
%
% intrellis: a description of the successor nodes using a list, e.g.
%          intrellis{1} = [1 3]; intrellis{2} = [3 4];
%          intrellis{3} = [1 2]; intrellis{4} = [3 4];
% inbranchweight: weights of branches used for comparison, saved as
%    cells in branchweight{state_number, branch_number}
%    branchweight may be a vector
%    e.g.  branchweight{1,1} = 0; branchweight{1,2} = 6;
%          branchweight{2,1} = 3; branchweight{2,2} = 3;
%          branchweight{3,1} = 6; branchweight{3,2} = 0;
%          branchweight{4,1} = 3; branchweight{4,2} = 3;
% inpathlen: length of window over which to compute
% normfun: the norm function used to compute the branch cost

interplane.m
% 
% find the intersecting point for the planes
%  m1'(x-x1) = 0   and  m2'(x-x2)=0.
%
% function x = interplane(m1,x1,m2,x2)
%
% m1 = normal to plane
% x1 = point on plane
% m2 = normal to plane
% x2 = point on plane
%
% x = intersecting point

invdiff.m
% 
% Compute the inverse differences for a rational interpolation function,
% returning the vector of inverse differences that are necessary for
% interpolation
%
% function phis = invdiff(ts,fs)
%
% ts = vector of independent variable
% fs = vector of dependent variable
%
% phis = inverse differences

jordanex.m
% example of Jordan forms

kaisfilt.m
% test the design of a Kaiser filter

kronsum.m
%
% Kronecker sum of A and B.  A and B are assumed square.
%
% function C = kronsum(A,B).

lagrangepoly.m
% 
% Lagrange interpolator

latexform.m
% 
% Display a matrix X in latex form
%
%function latexform(fid,X,[nohfill])
% 
% fid = output file id (use 1 for terminal display)
% X = matrix of vector to display
% nohfill = 1 if no hfill is wanted

lfsr.m
% 
% Produce m outputs of an lfsr with coefficient c and initial values y0
%
% function y = lfsr(c,y0,m)
% y0 = [y_0,...,y_{p-2},y_{p-1}]
% c = [c(1),c(2),...,c(p)]
% 
% y_j = sum_{i=1}^p y_{j-i} c(i)

lfsrfind.m
% 
% Find a good lfsr c to match Ac=b
%
% function c = lfsrfind(A,b)

lfsrfind2.m
% 
% Find a good lfsr c to match Ac=b
% where A and b are formed by the lfsr
% In this case, feed the error back around
%
% function c = lfsrfind2(y,m)

lsdata.m
% Make least-squares data matrices

makemarkov.m
% 
% Return the sequence of n impulse response samples into the cell array y
% y{1},y{2}, ... y{n}
%
% function y = makemarkov(A,B,C,n)
%
% (A,B,C) = system
% n = number of samples
%
% y = cell array of impulse responses

makeperm.m
% Return all permutations of length n

maketoeplitz.m
% 
% Form a toeplitz matrix from the input data y
%
% function [H] = maketoeplitz(y,m,n)
%
% y = input data = [y1 y2 ...] (a series of vectors in a _row_)
% m = number of block rows in H
% n = number of block columns in H

marv.m
% 
% Prony: given a sequence of (supposedly) pure sinusoidal data, 
% determine the frequency using model methods
%
% function f = marv(x,fs)
% x = data sequence
% fs = sampling rate
%
% f = frequencies found

masseyinit.m
% Initialize the iteratively called massey's algorithm
%
% function masseyinit()

masseyit.m
% 
% Compute the lfsr connection polynomial using Massey's algorithm
% accepting new data at each iteration.
% masseyinit should be called before calling this the first time
%
% y = new data point
%
% c = updated connection polynomial

miniapprox1.m
% minimax approximation example

modaldata1.m
% data for a modal analysis problem

myplaysnd.m
% 
% Modified and simplified from playsnd, to make the sample rate stuff work

n2b.m
% convert n to an m-bit binary representation

neville.m
% 
% Neville's algorithm for computing a value for an interpolating polynomial
% y = NEVILLE(x,X,Y) takes the (xi,yi) coordinate pairs in the
% vectors X and Y and computes the value of the unique
% interpolating polynomial that passes through the list of points
% at the given value of x. 
% 
% function y = neville(x, X, Y)
%
% x = interpolated point
% X = X data
% Y = Y data
%
% y = interpolated value

nntrain2.m
% 
% train a neural network using the input/output training data [x,d]
% with sequential selection of the data
%
% function w = nntrain(x,d,m,ninput,mu)
%
% x = [x(1) x(2) ... x(N)]   
% d= [d(1) d(2) ... d(N)]
% nlayer = number of layers
% m = number of neurons on each layer, 
%     m(1) = input layer, ... m(nlayer+1) = ouput layer
% mu = steepest descent step size
% alpha = (optional) momentum constant
% maxiter = (optional) maximum number of iterations (w = no maximum)
% w = (optional) starting weights
%
% err = (optional) total squared error from training

pade1.m
% Pade example

padefunct.m
% 
% Find the Pade approximation from the Maclaurin series coefficients
%
% function [A,B] = padefunct(c,m,k)
%
% c = Maclaurin series coefficients (need m+k+1)
% m = degree of numerator polynomial
% k = degree of denominator polynomial
%
% A = coefficients of numerator polynomial (in Matlab order)
% B = coefficients of denominator polynomial (in Matlab order)

permer.m
% 
% function permlist = permer(n1,p,perm,permnew,permlist)

pisexamp.m
% Example for Pisarenko Harmonic Decomposition

plotbernapprox.m
% plot the Benstein polynomial approximation to $f(t) = e^t$

plotbernpoly.m
% plot the Benstein polynomial

plotfplane.m
% plot a function and a linear approximating surface

plotplane.m
% determine points in the plane m'(x-x0) = 0 for plotting purposes
%

polyadd.m
%
% Add the polynomials p=a+b
% 
% function p = polyadd(a,b)
%
% a,b = polynomial
%
% p = polynomial sum.

polydiv.m
% 
% Divide a(x)/b(x), and return quotient and remainder in q and r
% Coefficients are assumed to be in Matlab standard order (highest order first)
%
%
% function [q,r] = polydiv(a,b)
%
% a = numerator
% b = denominator
%
% q = quotient
% r = remainder

polydivgfp.m
% 
% Divide a(x)/b(x), and return quotient and remainder in q and r
% using arithmetic in GF(p)
% Coefficients are assumed to be in Matlab standard order (highest order first)
%
% function [q,r] = polydivgfp(a,b)
%
% a = numerator
% b = denominator
%
% q = quotient
% r = remainder

polymult.m
% 
% Multipoly the polynomials p=a*b
%
% function p = polymult(a,b)
%
% a,b = polynomials
%
% p = product

polysub.m
% 
% Subtract the polynomials p=a-b
%
% a,b = polynomials
%
% p = difference

psdarma.m
%
% Plot the psd of an arma model
%
% function [w,h] = psdarma(b,a)
%
% b = numerator coefficients
% a = denominator coefficients
%
% w = frequency values
% h = absolute value of response

ratinterp.m
% 
% Compute the rational function interpolation
% from the data in ts and fs.
% Polynomial coefficients returned in Matlab order (largest to smallest)

ratinterp1.m
% 
% Compute a single interpolated point f(t) given the interpolating data
% and the inverse differences
%
% function f = ratinterp1(t,ts,fs,phis)
%
% t = point at which to evaluate
% ts = vector of independent data
% fs = vector of depdendent data
% phis = inverse differences
%
% f = interpolated value

ratintfilt.m
% Try some data for a rationally-interpolated filter

res.m
% 
% Computes <a^n>_m
%
% function d = res(a,n,m)
%
% a = value
% n = exponent
% m = modulo
%
% d = remainder(a^n,m0

schurcohn.m
% 
% Returns 1 if p is a Schur polynomial (all roots inside unit circle)
%
% function stable = schurcohn(p)
% 
% p = polynomial coefficients
%
% stable = 1 if stable polynomial

simppivot.m
% 
% Pivot a linear programming tableau about the p,q entry
% 
% function tableau = simppivot(intableau,p,q)
%
% intableau = tableau
% (p,q) = point about which to pivot
%
% tableau = pivoted tableau

solvlincong.m
% 
% Ddetermine the solution to the linear congruence
% a x equiv b (mod m), if it exists
%
% function x = solvlincong(a,m,b)

sreal.m
% sysreal.m
% data for the system identification example in the SVD stuff

sreal1.m
% SVD realization

sugiyama.m
% 
% Compute the GCD g = (b,c) using the Euclidean algorithm
% and return s,t such that bs+ct = g, where b and c are polynomials
% with real coefficients
%
% thresh = (optional) threshold argument used to truncate small remainders

sysidsvd2.m
% 
% given a sequence of impulse responses in h (a cell array)
% identify a system (A,B,C)
% This uses the tohankel method of finding a nearest hankel matrix
% of desired rank
%
% function [A,B,C] = sysidsvd(h,order)
%
% h = impulse response sequence (cell array)
% order = desired order of system
%
% (A,B,C) = system

taylorf.mm
(* example of a taylor series *)

tocrt.m
% 
% Compute the representation of the scalar x using the
% using the Chinese Remainder Theorem (CRT) with
% moduli m = [m1,m2,...,mr].  It is assumed (without checking)
% that the moduli are relatively prime
%
% function y = tocrt(x,m)
%
% x = number to convert
% m = set of moduli
%
% y = CRT representation of x

tocrtpoly.m
% 
% Compute the representation of the polynomial f using the
% using the Chinese Remainder Theorem (CRT) with
% moduli m = {m1,m2,...,mr}.  It is assumed (without checking)
% that the moduli are relatively prime.
% m is passed in as a cell array containing polynomial vectors
% and y is returned as a cell array containing polynomial vectors
%
% function y = tocrt(f,m)
%
% f = polynomial
% m = set of modulo polynomials
%
% y = CRT form of f

tohankelbig.m
% 
% Determine the matrix nearest to A which is (block) Hankel and has rank r
% using the composite mapping algorithm
%
% function A = tohankelbig(A,r)
%
% A = input matrix
% r = desired ranke
% d = (optional) block size (default=1)
%
% A = nearest rank r Hankel matrix
% diff = norm of difference between matrices

triginterp.m
% demonstrate trigonometric interpolation

vandsolve1.m
%
% Solves the equation Vx = fs, where V is the Vandermonde
% matrix determined from ts.
%
% function a = vandsolve1(ts,fs)
%
% ts = abscissa values
% fs = ordinate values
%
% a = solution

vitnop.m
% 
% Compute the norm of the difference between inputs
% This function may be feval'ed for use with the Viterbi algorithm
% In this case, the norm is simply taken as the branch number
%
% function d = vitnop(branch,input)
%

vitsqnorm.m
% 
% Compute the square norm of the difference between inputs
% This function may be feval'ed for use with the Viterbi algorithm
% (state and nextstate are not used here)
%
% function d = vitsqnorm(branch,input,state,nextstate)

wino3by3.m
% 
% Convolve the 3-sequence a with the 3-sequence b 
% a and b are both assumed to be column vectors
% using Winograd convolution
%
% function c = wino3by3(a,b)

winotest.m
% Set up data for a Winograd convolution algorithm

winotest2.m
% Set up data for a Winograd convolution algorithm

Directory: mkpict

attract1.m
% a plot showing an attractor

attract2.m
% a plot showing an attractor

bayes1.m
% Bayes decision tests

bayes2.m
% Bayes decision tests for Gaussian

bayes4.m
% show the decision regions for a 3-way test

binchan.m
% 
% Data for Bayesian detection on the binary channel

binchanex.m
% data for binary channel

chebyplot.m
% Plot Chebyshev polynomials

chi2plot.m
% 
% Plot the chi-squared r.v.

compmap3.m
% make figure comppos1

condhilb.m
% Plot the condition of the Hilbert matrix

drawtrellis.m
% 
% Draw a trellis in LaTeX picture mode
%
% function drawtrellis(fid,numbranch,r,p)
%
% fid = output file id
% numbranch = number of branches to draw
% r = path cost values
% p = flag
%
% Other values are contained in global variables.  See the file

drawtrelpiece.m
%
% Draw a piece of a trellis in LaTeX picture mode
%
% fname = file name
% trellis = trellis description
% branchweight = weights of branches

drawvit.m
% Program to draw the paths for the Viterbi algorithm using a LaTeX picture

duality1.m
% Make a plot illustrating duality

eigdir.m
% make a contour plot of eigenstuff

eigdir2.m
% make a contour plot of eigenstuff

eigdirex.m
% make a contour plot of eigenstuff

eigdist.m
% show the asymptotic equal distribution of eigenvalues

ellipse.m
% Plot contours of an ellipse with large eigenvalue disparity
% and the results of steepest descent

ellipsecg.m
% Plot contours of an ellipse with large eigenvalue disparity
% and the results of conjugate gradient.

entplot.m
% plot the binary entropy function

expmod.m
% Test Cadzow's results on the sinusoidal modeling

fourser.m
% example Fourier series

hilb1.m
% Program to generate the data for the hilbert approximation to
% the exponential function

ifs3.m
% Plot the logistic map and the orbit of a point

ifs3b.m
% Plot the logistic map and the orbit of a point
% do not specify lambda and x0 here: it is done by an upper script

ifs4.m
% Demonstrate the logistic map

ifsex3.m
% find an affine transformation Ax + b that transforms from
% {x00,x10,x20} to {x01,x11,x21}

ifsfig1.m
% Make side-by-side figures

legendreplot.m
% Plot legendre polynomials

makeim.m
% make a test image for tomography example

matcond.m
% Make an ill-conditioned matrix of sinusoids.

matcond2.m
% Set up an ill-conditioned matrix of sinusoids

min1.m
% make the contour plot for wftest

min2.m
% make the contour plot for wftest

moveiter.m
% test the solution of a moving RHS in the equation Ax=b

newt1.m
% Demonstrate newton's stuff

newt2.m
% Demonstrate newton's stuff on Rosenbrocks function

oddeven.m
% data for odd/even game

orthog.mma
(* sample file for orthogonalization *)

patrec1.m
% generate some simple pattern recognition example data

plotI0.m
% Plot the Bessel function

plotJsurf.m
% plot a quadratice error surface

plotbpsk.m
% Plot the probability of error for BPSK

plotgauss.m
% Plot the Gaussian function

plotgauss2.m
% Plot approximations to the central limit theorem

plotgauss3.m
% plot a Gauss surface plot

plotwavelet.m
% plot the wavelet data

roc1.m
% plot the roc for a gaussian r.v.

roc2.m
% plot the roc for a a xi^2

roc3.m
% plot the roc for a gaussian r.v. and its conjugate

rosenbrock.m
% Plot the Rosenbrock function contours

rosengrad.m
% 
% compute the gradient of the rosenbrock function for test purposes
% function grad = rosengrad(x)

saddle1.m
% make a saddle plot

scatter.m
% create a scatter plot to demonstrate principal component

scatterex.m
% create a scatter plot to demonstrate principal component

sigmoid.m
% plot the sigmoid function

steeperr.m
% Plot errors of the steepest descent

steeperrplot.m
% Make plots of error for steepest descent

steepest1.m
% Demonstrate steepest descent on Rosenbrocks function

sugitest.m
% test the Sugiyama algorithm

surf1.m
% make a surface plot

test2regress.m
% Test the formulas for regression in two dimensions
% input: x and y vectors

test2regress2.m
% Test the formulas for regression in two dimensions
% input: x and y vectors

testeigfil.m
% Test the eigenfilter stuff

testeigfil2.m
% Test the eigenfilter stuff

testeigfil3.m
% test the eigenfilter stuff

testexlms.m
% Test the lms in a system identification setting
% Assume a Gaussian input

testlms.m
% test the lms in an equalizer setting
% Assume a binary +/- 1 input.

testmusic.m
% Test the music algorithm

testnn1.m
% test the neural network stuff

testnn2.m
% test the neural network stuff
% (run testnn1.m first to get the network trained)
%
% does some plots after the initial training is finished

testnn3.m
% test the neural network stuff
% try different values of mu and alpha
% run testnn1 first to get the training data

testrls.m
% test the rls in an equalizer setting
% Assume a binary +/- 1 input.

testrls2.m
% test the rls in a system identification setting
% Assume a binary +/- 1 input.

testrls2ex.m
% test the rls in a system identification setting
% Assume a binary +/- 1 input.

testrlsex.m
% test the rls in an equalizer setting
% Assume a binary +/- 1 input.

testrot.m
% test the procrustes rotation

testtls.m
% Test tls stuff

vq1.m
% Generate random Gaussian data, determine a codebook for it, and plot

wftestcont.m
% make the contour plot for wftest

Directory: solutions

ator2.m
% 
% Given the coefficients from a 2nd-order AR model
% y[t+2] + a1 y[t+1] + a2 y[t] = f[t+2],
% where f has variance sigmaf2, compute sigma_y^2, r[1], and r[2].
% 
% function [sigma2,r1,r2] = ator2(a1,a2,sigmaf2)
%
% a1, a2 -- AR model coefficients
% sigmaf2 -- input noise variance
%
% sigma2 -- output noise variance
% r1, r2 -- covariance values

backdyn.m
%
% Backward dynamic programming
%
% function [pathlist,cost] = backdyn(H,W)
%
% H = graph
% W = costs
%
% pathlist = list of paths
% cost = cost of paths

backsub.m
% 
% solve Ux = b, where U is upper triangular
%
% function x = backsub(U,b)
% U = upper triangular matrix
% b = right and side
%
% x = solution

bayesest1.m
% Example of non-Gaussian Bayes estimate

correst.m
% 
% Estimate the autocorrelation function
% the returned values are offset (by Matlab requirements) so that
% r(1) = r[0], etc.
% Only correlations for positive lags are returned.  For other values,
% use the fact that r[k] = conj(r[-k])
%
% function r = correst(x)
%
% x = data sequence
%
% r = estimated correlations

fcm.m
%
% Find k clusters on the data X using fuzzy clustering
% function [Y,d] = fcm(X,m)
%
% X = input data: each column is a training data vector
% m = number of clusters to find
%
% Y = set of clusters: each column is a cluster centroid
% U = membership functions

findperm.m
% 
% Determine a permutation P such that Px = z (as close as possible)
% using a composite mapping
%
% function P = findperm(x,z,maxiter)
%
% x = input data (to be permuted)
% z = permuted data
% maxiter = optional argument on number of iterations
%
% P is an index listing, so the permutation is obtained by x(P) = z.

forbacksub.m
% 
% Solve Ax = b, where A has been factored as PA = LU
%
% function [x] = forbacksub(b,LUin,indx)

forbacksubround.m
% 
% Solve Ax = b, where A has been factored as PA = LU
% using digits places after the decimal point.
%
% function [x] = forbacksubround(b,LUin,indx,digits)
%
% b = right hand side
% LUin = LUfactorization (from newlu)
% indx = pivot index list (from newlu)
% digits = number of digits to retain
%
% x = solution

forsub.m
% 
% Solve Lx = b, where L is lower triangular
% %
% function x = forsub(L,b)
%
% L = lower triangular matrix
% b = right-hand side
%
% x = solution

getulu.m
% 
% Return the L and U matrix from the LU factorization computed by newlu
%
% function [L,U] = getulu(luin,indx)
%
% luin = lu matrix from newlu
% index = pivot index from newlu
%
% L = lower triangle
% U = upper triangle

gramschmidt2.m
% 
% Modified Gram-Schmidt: Compute the Gram-Schmidt orthogonalization of the 
% columns of A, A = QR
%
% function [Q,R] = gramschmidt2(A)
%
% A = matrix to be factored
%
% Q = orthogonal matrix
% R = upper triangular matrix

gramschmidtW.m
% 
% Compute the Gram-Schmidt orthogonalization of the 
% columns of A with the inner product <x,y> = x'*W*y
% W should be symmetric
% 
% function [Q,R] = gramschmidtW(A,W)
%
% A = matrix to be factored
% W = weighting matrix
%
% Q = orthogonal matrix
% R = upper triangular matrix

gs.m
% 
% Gram-Schmidt using symbolic toolbox
%
% function q = gs(P,a,b,w,t)
%
% P = list of functions in P(1), P(2), ...
% a = lower limit of integration
% b = upper limit of integration
% w = weight function
% t = variable of integration
%
% q = array of orthogonal functions

gs.mma
(* A Gram-Schmidt procedure *)

ifs1.m
% test some ifs stuff

mgs.m
% 
% Compute the Gram-Schmidt orthogonalization of the 
% columns of A, assuming nonzero columns of A
% using the modified Gram-Schmidt algorithm.  A = QR
% 
% function [Q,R] = mgs(A)
%
% A = matrix to be factored
%
% Q = orthogonal matrix
% R = upper triangular

newluround.m
% 
% Compute the lu factorization of A
% controlling the pivoting and the rounding
% dopivot = 1 if piviting desired
% digits = number of decimal places to retain in rounding (digits=3 for 3
% dec. digits)
%
% function [lu,indx] = newluround(A,dopivot,digits)
%
% A = matrix to factor
% dopivot = flag if pivoting desired
% digits = number of digits to retain
%
% lu = matrix containg L and U factors
% indx = index of pivot permutations

newprony.m
% 
% Prony's method: Given a sequence of supposedly sinusoidal data with p
% modes, determine the  vector a of characteristic equation coefficients and
% modes --- the roots of the characteristic polynomial 
%
% function [a,r] = newprony(x,p)
% 
% x = sinusoidal data vector
% p = number of modes
%
% a = characteristic polynomial 
% r = roots of the characteristic polynomial

plotellipse.m
% 
% Determine the points to plot an ellispe in two dimensions, 
% described by (x-x0)'*A*(x-x0) = c, where A is symmetric
%
%  function [x] = plotellipse(A,x0,c)
%
% A = symmetric matrix describing ellipse
% x0 = center point
% c = constant
%
% x = 2 x n list of data points

plotlfsrautoc.m
% plot the autocorrelation of the output of an LFSR

rdigits.m
% 
% round the input to digits places
%
% function x = rdigits(y,digits)
%
% y = input value
% digits = number of digits
%
% x = rounded value

speceig.m
% Set up a matrix with specified eigenvalue strucure

speceig1.m
% Construct a matrix with given eigenspace structure

speceig2.m
% Construct a matrix with given eigenspace structure

wrls.m
% 
% Given a scalar input signal x and a desired scalar signal d,
% compute an RLS update of the weight vector h.
% eap is an optional return parameter, the a-priori estimation error
% This function must be initialized by wrlsinit
%
% function [h,eap] = wrls(x,d)
%
% x = scalar input
% d = desired value
%
% h = filter coefficient vector
% eap = a priori estimation error

wrlsinit.m
% 
% Initialize the weighted RLS filter
%
% function rlsinit(m,delta)
%
% m = dimension of vector
% delta = a small positive constant used for initial correlation inverse
% lambdain = value of lambda to use for decay factor

Directory: testers

arttest.m
% test the ART algorithm

bliter.m
%
% test reconstruction of bandlimited function

conjgradtest.m
% Test the conjugate gradient algorithm

dijtest.m
% Test Sisjstras and Warshall's algorithms

givens1.m
% test a Givens matrix

hmmtest.m
% Test data for the HMM

hmmtest2.m
% Test data for the HMM

hmmtest21.m

hmmtest22.m

hmmtest23.m

hmmtest24.m
% axis([1 4 -25 -5])

hmmtest2v.m
% Test the HMM using both Viterbi and EM-algorithm based training methods
%  the file hmmtest2vb.m contains identical results, but without the
%  plotting instructions at the end.

hmmvit.m
% test the HMM stuff using the VA

houseW.m
% test the weighted Householder idea

lpfilttest.m
% test the linear programming filter design

testbackdyn.m
% test the forward dynamic programming algorithm

testeigcomp.m
% Test the eigencomputation routine

testfblp.m
% Test the fblp AR estimator

testfindperm.m
% Test the findperm program

testfordyn.m
% Test the forward dynamic programming algorithm

testkarmarker.m
% testkarmarker.m: Test the Karmarker linear programming solution

testmarv.m
% Test the Prony's method in marv.m

testmassey.m
% a sequence with a linear complexity of 5

testmassey2.m
% Test the routines to find lfsrs for data compression

testpolydiv.m
% test the polynomial division stuff

testprim5.m
% generate elements in GF(5^2) using p(x) = x^2 + x+2

testprony.m
% Test the prony method

vit2test.m
% Test and plot the results of the viterbi algorithm

vittest.m
% 
% Test and plot the results of the viterbi
%
% function vittest()

vittest1.m
% demonstrate the VA

vittest2.m
%
% Test and plot the results of the viterbi
%
%  function vittest()

vittest3.m
%
% Ttest and plot the results of the viterbi
% 
% function vittest()

vittest4.m
% Test and plot the results of the viterbi

vittest5.m
% 
% Test and plot the results of the viterbi
% 
% function vittest()

wavetest2.m
% Test some computations for wavelet transforms (matrix-based)

wftest2.m
% test the Wiener filter equalizer stuff for a first-order signal and
% first-order channel

youyouhu / ecen671-book-matlab Goto Github PK

ecen671-book-matlab's Introduction

MATLAB for Math of Signals/Systems Book

DIRECTORIES:

TOOLBOXES:

ecen671-book-matlab's People

Contributors

Forkers

Recommend Projects

React

Vue.js

Typescript

TensorFlow

Django

Laravel

D3

Recommend Topics

javascript

web

server

Machine learning

Visualization

Game

Recommend Org

Facebook

Microsoft

Google

Alibaba

D3

Tencent