Mercurial > hg > octave-terminal
changeset 559:4e826edfbc56
[project @ 1994-07-25 22:18:28 by jwe]
Initial revision
author | jwe |
---|---|
date | Mon, 25 Jul 1994 22:19:05 +0000 |
parents | faf108b99d21 |
children | 34c713db72b6 |
files | scripts/general/prepad.m scripts/image/colormap.m scripts/image/gray.m scripts/image/gray2ind.m scripts/image/image.m scripts/image/imagesc.m scripts/image/imshow.m scripts/image/ind2gray.m scripts/image/ind2rgb.m scripts/image/loadimage.m scripts/image/ntsc2rgb.m scripts/image/ocean.m scripts/image/octtopnm.c scripts/image/rgb2ind.m scripts/image/rgb2ntsc.m scripts/image/saveimage.m scripts/polynomial/polyinteg.m scripts/polynomial/residue.m scripts/polynomial/roots-amr.m scripts/set/complement.m scripts/set/create_set.m scripts/set/intersection.m scripts/set/union.m scripts/signal/sinc.m |
diffstat | 24 files changed, 1280 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/scripts/general/prepad.m @@ -0,0 +1,43 @@ +function y = prepad(x,l,c) +#prepad(x,l) +#Prepends zeros to the vector x until it is of length l. +#prepad(x,l,c) prepends the constant c instead of zero. +# +#If length(x) > l, elements from the beginning 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: prepad(x,l) or prepad(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(lx-l+1:lx); + else + if(rows(x)>1) + y = [ c*ones(l-lx,1); x ]; + else + y = [ c*ones(1,l-lx) x ]; + endif + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/image/colormap.m @@ -0,0 +1,41 @@ +function cmap = colormap(map) +#Set the current colormap. +# +#colormap(map) sets the current colormap to map. map should be an n row +#by 3 column matrix. The columns contain red, green, and blue intensities +#respectively. All entries should be between 0 and 1 inclusive. The new +#colormap is returned. +# +#colormap("default") restores the default colormap (a gray scale colormap +#with 64 entries). The default colormap is returned. +# +#colormap with no arguments returns the current colormap. + +#Author: +# Tony Richardson +# amr@mpl.ucsd.edu +# July 1994 + + global CURRENT_COLOR_MAP + + cmap_name = "CURRENT_COLOR_MAP"; + + if(nargin == 1) + if(isstr(map)) + if(strcmp(map,"default")) + CURRENT_COLOR_MAP = gray; + else + error("invalid argument"); + endif + else + # Set the new color map + CURRENT_COLOR_MAP = map; + endif + elseif(exist(cmap_name) == 0) + # If global color map doesn't exist, create the default map. + CURRENT_COLOR_MAP = gray; + endif + + # Return current color map. + cmap = CURRENT_COLOR_MAP; +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/image/gray.m @@ -0,0 +1,12 @@ +function map = gray(number) + + if(nargin == 0) + number = 64; + endif + + gr = [0:(number-1)]'; + + + map = [ gr gr gr ]/(number-1); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/image/gray2ind.m @@ -0,0 +1,11 @@ +function [X, map] = gray2ind(I,n) + + if(nargin == 1) + n = 64; + endif + + map = gray(n); + + X = round(I*(n-1)) + 1; + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/image/image.m @@ -0,0 +1,53 @@ +function image(x, zoom) +#Display an octave image matrix. +# +#image(x) displays a matrix as a color image. The elements of x are indices +#into the current colormap and should have values between 1 and the length +#of the colormap. +# +#image(x,zoom) changes the zoom factor. The default value is 4. +# +#SEE ALSO: imshow, imagesc, colormap. + +#Author: +# Tony Richardson +# amr@mpl.ucsd.edu +# July 1994 +# +#Modified: +# Tony Richardson +# amr@mpl.ucsd.edu +# July 1994 +# (Modifications based on suggestions from John Eaton.) + + global IMAGEDIR + + if (nargin == 0) + # Load Bobbie Jo Richardson (Born 3/16/94) + x = loadimage([IMAGEDIR,"/default.img"]); + zoom = 2; + elseif(nargin == 1) + zoom = 4; + elseif(nargin > 2) + error("usage: image (matrix, [zoom])"); + endif + + # Generate random file name + rnd_str = num2str(fix(rand*10000)); + ppm_name = ["image.", rnd_str, ".ppm" ]; + + saveimage(ppm_name,x,"ppm"); + + # Start the viewer + # Try xv, then xloadimage. + + xv = sprintf("xv -expand %f %s",zoom,ppm_name); + xloadimage = sprintf("xloadimage -zoom %f %s",zoom*100, ppm_name); + rm = sprintf("rm -f %s",ppm_name); + + command = sprintf("( %s || %s && %s ) > /dev/null 2>&1 &", ... + xv, xloadimage, rm); + + shell_cmd(command); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/image/imagesc.m @@ -0,0 +1,47 @@ +function x = imagesc(x, zoom) +#Scale and display a matrix as an image. +# +#imagesc(x) displays a scaled version of the matrix x. The matrix is +#scaled so that its entries are indices into the current colormap. +#The scaled matrix is returned. +# +#imagesc(x,zoom) sets the magnification, the default value is 4. +# +#SEE ALSO: image, imshow + +#Author: +# Tony Richardson +# amr@mpl.ucsd.edu +# July 1994 +# +#Modified: +# Tony Richardson +# amr@mpl.ucsd.edu +# July 1994 +# (Modifications based on suggestions from John Eaton.) + + if (nargin < 1 || nargin > 2) + error("usage: image (matrix, [zoom])"); + endif + + if (nargin == 1) + zoom = 4; + endif + + [ high, wide ] = size(x); + + maxval = max(max(x)); + minval = min(min(x)); + + # Rescale matrix so that all values are in the range 0 to + # length(colormap) inclusive + if (maxval == minval) + x = ones(high, wide); + else + # Rescale values to between 1 and length(colormap) inclusive. + x = fix((x - minval)/(maxval - minval) * (length(colormap)-1)) + 1; + endif + + image(x,zoom); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/image/imshow.m @@ -0,0 +1,33 @@ +function imshow(a1,a2,a3) +#Display images. +# +#imshow(X) displays an indexed image using the current colormap. +# +#imshow(X,map) displays an indexed image using the specified colormap. +# +#imshow(I,n) displays a gray scale intensity image. +# +#imshow(R,G,B) displays an RGB image. +# +#SEE ALSO: image, imagesc, colormap, gray2ind, rgb2ind. + +#Author: +# Tony Richardson +# amr@mpl.ucsd.edu +# July 1994 + + if (nargin == 0) + error("usage: imshow requires at least one argument."); + elseif(nargin == 2) + if(length(a2)==1) + [a1 a2] = gray2ind(a1,a2); + endif + colormap(a2); + elseif(nargin == 3) + [a1 a2] = rgb2ind(a1,a2,a3); + colormap(a2); + endif + + image(a1); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/image/ind2gray.m @@ -0,0 +1,35 @@ +function Y = ind2gray(X,map) +#Convert an octave indexed image to a gray scale intensity image. +# +#Y = ind2gray(X) converts an indexed image to a gray scale intensity +#image. The current colormap is used to determine the intensities. +#The intensity values lie between 0 and 1 inclusive. +# +#Y = ind2gray(X,map) uses the specified colormap instead of the current +#one in the conversion process. +# +#SEE ALSO: gray2ind, rgb2ntsc, image, colormap + + if (nargin == 1) + map = colormap; + endif + + # Convert colormap to intensity values. + yiq = rgb2ntsc(map); + y = yiq(:,1); + + # We need Fortran indexing capability, but be sure to save the user's + # preference. + pref = do_fortran_indexing; + do_fortran_indexing = "true"; + + # Replace indices in the input matrix with indexed values in the output + # matrix. + [rows, cols] = size(X); + Y = y(X(:)); + Y = reshape(Y,rows,cols); + + # Restore the user's preference. + do_fortran_indexing = pref; + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/image/ind2rgb.m @@ -0,0 +1,29 @@ +function [R G B] = ind2rgb(X,map) +#Convert an indexed image to red, green, and blue color components. +# +#[R G B] = ind2rgb(X) uses the current colormap for the conversion. +# +#[R G B] = ind2rgb(X,map) uses the specified colormap. +# +#SEE ALSO: rgb2ind, image, imshow, ind2gray, gray2ind. + + if(nargin == 1) + map = colormap; + endif + + [hi wi] = size(X); + + pref = do_fortran_indexing; + do_fortran_indexing = "true"; + + R = map(X(:),1); + G = map(X(:),2); + B = map(X(:),3); + + R = reshape(R,hi,wi); + G = reshape(G,hi,wi); + B = reshape(B,hi,wi); + + do_fortran_indexing = pref; + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/image/loadimage.m @@ -0,0 +1,12 @@ +function [X, map] = loadimage(filename) +#Load an image file. +# +#[X, map] = loadimage(img_file) loads an image and it's associated color +#map from file img_file. The image must be in stored in octave's image +#format. +# +#SEE ALSO: saveimage, load, save + + eval(['load ', filename]); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/image/ntsc2rgb.m @@ -0,0 +1,7 @@ +function rgb = ntsc2rgb(yiq) + + trans = [ 1. 1. 1.; 0.95617 -0.27269 -1.10374; 0.62143 -0.64681 1.70062 ]; + + rgb = yiq * trans; + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/image/ocean.m @@ -0,0 +1,18 @@ +function map = ocean(number) + + if(nargin == 0) + number = 64; + endif + + cutin = fix(number/3); + + + dr = (number-1)/(cutin); + r = prepad([0:dr:(number-1)],number)'; + dg = (number-1)/(2*cutin); + g = prepad([0:dg:(number-1)],number)'; + b = [0:(number-1)]'; + + map = [ r g b ]/(number-1); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/image/octtopnm.c @@ -0,0 +1,302 @@ +# /* +cc -s -o octtopnm octtopnm.c +exit +*/ + +#include <stdio.h> +#include <string.h> +#include <malloc.h> + +/* usage: octtopnm [-a] octfile */ + +static void usage(message) +char *message; +{ + if(message != NULL) { + fprintf(stderr,"octtopnm: %s\n",message); + } + fprintf(stderr,"usage: octtopnm [-a] octavefile\n"); + exit(1); +} + +static void fatal(message) +char *message; +{ + if(message != NULL) { + fprintf(stderr,"octtopnm: %s\n",message); + } + exit(1); +} + +int main(argc, argv) +int argc; +char **argv; +{ + int rawbits = 1, row, col, index; + int cmap_rows, cmap_cols, img_rows, img_cols; + int gray, pbm, pgm, ppm; + unsigned char **rgb, byte; + unsigned short **img; + char *oct_file_name; + FILE *oct_file; + char cmap_name[4], cmap_type[7], img_name[2], img_type[7]; + double mat_val; + int option; + extern char *optarg; + extern int optind; + + if(argc == 1) { + usage(NULL); + } + + while((option = getopt(argc,argv,"ha")) != EOF) { + switch(option) { + case 'h': + /* help */ + usage(NULL); + break; + case 'a': + rawbits = 0; + break; + case '?': + default: + usage("unrecognized option"); + } + } + + if(optind+1 != argc) { + usage("input file name missing"); + } + + oct_file_name = argv[optind]; + if((oct_file = fopen(oct_file_name,"r")) == NULL) { + fatal("unable to open input file"); + } + + if(fscanf(oct_file,"# name: %s\n",cmap_name) != 1 || + strcmp(cmap_name,"map") != 0) { + fatal("not a valid octave image file"); + } + + if(fscanf(oct_file,"# type: %s\n",cmap_type) != 1 || + strcmp(cmap_type,"matrix") != 0) { + fatal("not a valid octave image file"); + } + + if(fscanf(oct_file,"# rows: %d\n",&cmap_rows) != 1) { + fatal("error reading octave image file"); + } + + if(fscanf(oct_file,"# columns: %d\n",&cmap_cols) != 1) { + fatal("error reading octave image file"); + } + + if(cmap_cols != 3) { + fatal("invalid color map in octave image file"); + } + + if((rgb = (unsigned char **) + malloc(cmap_rows*sizeof(unsigned char *))) == NULL) { + fatal("out of memory"); + } + + if((rgb[0] = (unsigned char *) + malloc(cmap_rows*cmap_cols*sizeof(unsigned char))) == NULL) { + fatal("out of memory"); + } + + for(row=1; row<cmap_rows; row++) { + rgb[row] = rgb[row-1]+3; + } + + gray = 1; + for(row=0; row<cmap_rows; row++) { + for(col=0; col<cmap_cols; col++) { + if(fscanf(oct_file,"%lf",&mat_val) != 1) { + fatal("error reading color map entries"); + } + if(mat_val < 0) mat_val = 0.; + if(mat_val > 1) mat_val = 1.; + rgb[row][col] = mat_val*255; + } + if(gray) { + if(rgb[row][0] != rgb[row][1] || rgb[row][0] != rgb[row][2]) { + /* It's a color image. */ + gray = 0; + } + } + } + + if(fscanf(oct_file,"\n# name: %s\n",img_name) != 1 || + strcmp(img_name,"X") != 0) { + fatal("not a valid octave image file"); + } + + if(fscanf(oct_file,"# type: %s\n",img_type) != 1 || + strcmp(img_type,"matrix") != 0) { + fatal("not a valid octave image file"); + } + + if(fscanf(oct_file,"# rows: %d\n",&img_rows) != 1) { + fatal("error reading octave image file"); + } + + if(fscanf(oct_file,"# columns: %d\n",&img_cols) != 1) { + fatal("error reading octave image file"); + } + + if((img = (unsigned short **) + malloc(img_rows*sizeof(unsigned short *))) == NULL) { + fatal("out of memory"); + } + + if((img[0] = (unsigned short *) + malloc(img_rows*img_cols*sizeof(unsigned short))) == NULL) { + fatal("out of memory"); + } + + for(row=1; row<img_rows; row++) { + img[row] = img[row-1]+img_cols; + } + + for(row=0; row<img_rows; row++) { + for(col=0; col<img_cols; col++) { + if(fscanf(oct_file,"%lf",&mat_val) != 1) { + fatal("error reading color map entries"); + } + if(mat_val < 1) mat_val = 1.; + if(mat_val > cmap_rows) mat_val = cmap_rows; + img[row][col] = mat_val; + } + } + + pbm = pgm = ppm = 0; + + if(cmap_rows == 2 && gray && + ((rgb[0][0] == 0 && rgb[1][0] == 255) || + (rgb[0][0] == 255 && rgb[1][0] == 0))) { + /* Create a bitmap only if there are two colormap entries and they are + black and white. */ + pbm = 1; + } + else if(gray) { + /* If not a bitmap, create a gray scale image if the entries within + each row of the color map are equal. */ + pgm = 1; + } + else { + /* Otherwise create a full color image. */ + ppm = 1; + } + + if(rawbits) { + if(pbm) { + printf("P4\n"); + printf("%d %d\n",img_cols,img_rows); + index = 0; + for(row=0; row<img_rows; row++) { + for(col=0; col<img_cols; col++) { + if(index == 7) { + byte = (byte << 1) + !rgb[img[row][col]-1][0]; + fwrite(&byte,sizeof(unsigned char),1,stdout); + byte = 0; + index = 0; + } + else { + byte = (byte << 1) + !rgb[img[row][col]-1][0]; + index++; + } + } + } + if(index != 0) { + printf("\n"); + } + } + else if(pgm) { + printf("P5\n"); + printf("%d %d\n",img_cols,img_rows); + printf("255\n"); + for(row=0; row<img_rows; row++) { + for(col=0; col<img_cols; col++) { + fwrite(rgb[img[row][col]-1],sizeof(unsigned char),1,stdout); + } + } + } + else { + printf("P6\n"); + printf("%d %d\n",img_cols,img_rows); + printf("255\n"); + for(row=0; row<img_rows; row++) { + for(col=0; col<img_cols; col++) { + fwrite(rgb[img[row][col]-1],sizeof(unsigned char),3,stdout); + } + } + } + } + else { + if(pbm) { + printf("P1\n"); + printf("%d %d\n",img_cols,img_rows); + index = 0; + for(row=0; row<img_rows; row++) { + for(col=0; col<img_cols; col++) { + if(index == 30) { + printf("%d\n",!rgb[img[row][col]-1][0]); + index = 0; + } + else { + printf("%d ",!rgb[img[row][col]-1][0]); + index++; + } + } + } + if(index != 0) { + printf("\n"); + } + } + else if(pgm) { + printf("P2\n"); + printf("%d %d\n",img_cols,img_rows); + printf("255\n"); + index = 0; + for(row=0; row<img_rows; row++) { + for(col=0; col<img_cols; col++) { + if(index == 12) { + printf("%d\n",rgb[img[row][col]-1][0]); + index = 0; + } + else { + printf("%d ",rgb[img[row][col]-1][0]); + index++; + } + } + } + if(index != 0) { + printf("\n"); + } + } + else { + printf("P3\n"); + printf("%d %d\n",img_cols,img_rows); + printf("255\n"); + index = 0; + for(row=0; row<img_rows; row++) { + for(col=0; col<img_cols; col++) { + if(index == 4) { + printf("%d %d %d\n",rgb[img[row][col]-1][0], + rgb[img[row][col]-1][1],rgb[img[row][col]-1][2]); + index = 0; + } + else { + printf("%d %d %d ",rgb[img[row][col]-1][0], + rgb[img[row][col]-1][1],rgb[img[row][col]-1][2]); + index++; + } + } + } + if(index != 0) { + printf("\n"); + } + } + } +}
new file mode 100644 --- /dev/null +++ b/scripts/image/rgb2ind.m @@ -0,0 +1,31 @@ +function [X map] = rgb2ind(R,G,B) +#Convert and RGB image to an octave indexed image. +# +#[X map] = rgb2ind(R,G,B) +# +#SEE ALSO: ind2rgb, rgb2ntsc. +# +#Bugs: The color map may have duplicate entries. + + if(nargin != 3) + error("usage: [X, map] = rgb2ind(R,G,B)"); + endif + + [hi wi] = size(R); + + X = zeros(hi,wi); + + map = zeros(hi*wi,3); + + pref = do_fortran_indexing; + do_fortran_indexing = "true"; + + map(:,1) = R(:); + map(:,2) = G(:); + map(:,3) = B(:); + + X(:) = 1:(hi*wi); + + do_fortran_indexing = pref; + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/image/rgb2ntsc.m @@ -0,0 +1,6 @@ +function yiq = rgb2ntsc(rgb) + + trans = [ 0.299 0.596 0.211; 0.587 -0.274 -0.523; 0.114 -0.322 0.312 ]; + yiq = rgb * trans; + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/image/saveimage.m @@ -0,0 +1,75 @@ +function saveimage(filename,X,img_form,map) +#Save a matrix to disk in image format. +# +#saveimage(filename,x) saves matrix x to file filename in octave's image +#format. The current colormap is saved in the file also. +# +#saveimage(filename,x,"img") saves the image in the default format and +#is the same as saveimage(filename,x). +# +#saveimage(filename,x,"ppm") saves the image in ppm format instead of +#the default octave image format. +# +#saveimage(filename,x,"ps") saves the image in PostScript format instead +#of the default octave image format. (Note: images saved in PostScript format +#can not be read back into octave with loadimage.) +# +#saveimage(filename,x,format,map) saves the image along with the specified +#colormap in the specified format. +# +#Note: If the colormap contains only two entries and these entries are black +#and white, the bitmap ppm and PostScript formats are used. If the image is +#a gray scale image (the entries within each row of the colormap are equal) +#the gray scale ppm and PostScript image formats are used, otherwise the full +#color formats are used. +# +#SEE ALSO: loadimage, save, load, colormap + + global IMAGEDIR + + if(nargin < 2) + error("usage: saveimage(filename,matrix,[format, [colormap]])"); + elseif(nargin == 2) + if(!isstr(filename)) + error("File name must be a string."); + endif + map = colormap; + img_form = "img"; + elseif(nargin == 3) + if(!isstr(img_form)) + error("Image format specification must be a string"); + endif + map = colormap; + endif + + if (strcmp(img_form,"img") == 1) + oct_file = filename; + elseif (strcmp(img_form,"ppm") == 1) + # We need a temporary octave image file name. + oct_file = sprintf("image.%s.img",num2str(fix(rand*10000))); + ppm_file = filename; + elseif (strcmp(img_form,"ps") == 1) + # We need a temporary octave image file name. + oct_file = sprintf("image.%s.img",num2str(fix(rand*10000))); + ps_file = filename; + endif + + # Save image in octave image file format + eval(['save ', oct_file, ' map X']); + + # Convert to another format if requested. + if (strcmp(img_form,"ppm") == 1) + octtopnm = sprintf([IMAGEDIR,"/octtopnm %s > %s"],oct_file,filename); + rm = sprintf("rm -f %s",oct_file); + shell_cmd(octtopnm); + shell_cmd(rm); + elseif (strcmp(img_form,"ps") == 1) + octtopnm = sprintf([IMAGEDIR,"/octtopnm %s"],oct_file); + ppmtops = sprintf("pnmtops > %s 2> /dev/null", filename); + octtops = [ octtopnm, " | ", ppmtops ]; + rm = sprintf("rm -f %s",oct_file); + shell_cmd(octtops); + shell_cmd(rm); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/polyinteg.m @@ -0,0 +1,38 @@ +function p = polyinteg(p) +#polyinteg(c) +#Returns the coefficients of the integral the polynomial whose coefficients +#are represented by the vector c. +# +#The constant of integration is zero. +# +#SEE ALSO: poly, polyderiv, polyreduce, roots, conv, deconv, residue, +# filter, polyval, polyvalm + +# Author: +# Tony Richardson +# amr@mpl.ucsd.edu +# June 1994 + + if(nargin != 1) + error("usage: polyinteg(vector)"); + endif + + if(is_matrix(p)) + error("argument must be a vector"); + endif + + lp = length(p); + + if(lp == 0) + p = []; + return; + end + + if(rows(p) > 1) + # Convert to column vector + p = p.'; + endif + + p = [ p 0 ] ./ [lp:-1:1 1]; + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/residue.m @@ -0,0 +1,268 @@ +function [r, p, k, e] = residue(b,a,toler) +#[r p k e] = residue(b,a) +#If b and a are vectors of polynomial coefficients, then residue +#calculates the partial fraction expansion corresponding to the +#ratio of the two polynomials. The vector r contains the residue +#terms, p contains the pole values, k contains the coefficients of +#a direct polynomial term (if it exists) and e is a vector containing +#the powers of the denominators in the partial fraction terms. +#Assuming b and a represent polynomials P(s) and Q(s) we have: +# +# P(s) M r(m) N +# ---- = # ------------- + # k(n)*s^(N-n) +# Q(s) m=1 (s-p(m))^e(m) n=1 +# +#(# represents summation) where M is the number of poles (the length of +#the r, p, and e vectors) and N is the length of the k vector. +# +#[r p k e] = residue(b,a,tol) +#This form of the function call may be used to set a tolerance value. +#The default value is 0.001. The tolerance value is used to determine +#whether poles with small imaginary components are declared real. It is +#also used to determine if two poles are distinct. If the ratio of the +#imaginary part of a pole to the real part is less than tol, the imaginary +#part is discarded. If two poles are farther apart than tol they are +#distinct. +# +#Example: +# b = [1 1 1]; +# a = [1 -5 8 -4]; +# +# [r p k e] = residue(b,a) +# +# returns +# +# r = [-2 7 3]; p = [2 2 1]; k = []; e = [1 2 1]; +# +# which implies the following partial fraction expansion +# +# s^2 + s + 1 -2 7 3 +# ------------------- = ----- + ------- + ----- +# s^3 - 5s^2 + 8s - 4 (s-2) (s-2)^2 (s-1) +# +#SEE ALSO: poly, roots, conv, deconv, polyval, polyderiv, polyinteg + +# Author: +# Tony Richardson +# amr@mpl.ucsd.edu +# June 1994 + +# Here's the method used to find the residues. +# The partial fraction expansion can be written as: +# +# +# P(s) D M(k) A(k,m) +# ---- = # # ------------- +# Q(s) k=1 m=1 (s - pr(k))^m +# +# (# is used to represent a summation) where D is the number of +# distinct roots, pr(k) is the kth distinct root, M(k) is the +# multiplicity of the root, and A(k,m) is the residue cooresponding +# to the kth distinct root with multiplicity m. For example, +# +# s^2 A(1,1) A(2,1) A(2,2) +# ------------------- = ------ + ------ + ------- +# s^3 + 4s^2 + 5s + 2 (s+2) (s+1) (s+1)^2 +# +# In this case there are two distinct roots (D=2 and pr = [-2 -1]), +# the first root has multiplicity one and the second multiplicity +# two (M = [1 2]) The residues are actually stored in vector format as +# r = [ A(1,1) A(2,1) A(2,2) ]. +# +# We then multiply both sides by Q(s). Continuing the example: +# +# s^2 = r(1)*(s+1)^2 + r(2)*(s+1)*(s+2) + r(3)*(s+2) +# +# or +# +# s^2 = r(1)*(s^2+2s+1) + r(2)*(s^2+3s+2) +r(3)*(s+2) +# +# The coefficients of the polynomials on the right are stored +# in a row vector called rhs, while the coefficients of the +# polynomial on the left is stored in a row vector called lhs. +# If the multiplicity of any root is greater than one we'll +# also need derivatives of this equation of order up to the +# maximum multiplicity minus one. The derivative coefficients +# are stored in successive rows of lhs and rhs. +# +# For our example lhs and rhs would be: +# +# | 1 0 0 | +# lhs = | | +# | 0 2 0 | +# +# | 1 2 1 1 3 2 0 1 2 | +# rhs = | | +# | 0 2 2 0 2 3 0 0 1 | +# +# We then form a vector B and a matrix A obtained by evaluating the +# polynomials in lhs and rhs at the pole values. If a pole has a +# multiplicity greater than one we also evaluate the derivative +# polynomials (successive rows) at the pole value. +# +# For our example we would have +# +# | 4| | 1 0 0 | | r(1) | +# | 1| = | 0 0 1 | * | r(2) | +# |-2| | 0 1 1 | | r(3) | +# +# We then solve for the residues using matrix division. + + if(nargin < 2 || nargin > 3) + error("usage: residue(b,a[,toler])"); + endif + + if (nargin == 2) + # Set the default tolerance level + toler = .001; + endif + + # Make sure both polynomials are in reduced form. + a = polyreduce(a); + b = polyreduce(b); + + b = b/a(1); + a = a/a(1); + + la = length(a); + lb = length(b); + + # Handle special cases here. + if(la == 0 || lb == 0) + k = r = p = e = []; + return; + elseif (la == 1) + k = b/a; + r = p = e = []; + return; + endif + + # Find the poles. + p = roots(a); + lp = length(p); + + # Determine if the poles are (effectively) real. + index = find(abs(imag(p) ./ real(p)) < toler); + if (length(index) != 0) + p(index) = real(p(index)); + endif + + # Find the direct term if there is one. + if(lb>=la) + # Also returns the reduced numerator. + [k, b] = deconv(b,a); + lb = length(b); + else + k = []; + endif + + if(lp == 1) + r = polyval(b,p); + e = 1; + return; + endif + + + # We need to determine the number and multiplicity of the roots. + # D is the number of distinct roots. + # M is a vector of length D containing the multiplicity of each root. + # pr is a vector of length D containing only the distinct roots. + # e is a vector of length lp which indicates the power in the partial + # fraction expansion of each term in p. + + # Set initial values. We'll remove elements from pr as we find + # multiplicities. We'll shorten M afterwards. + e = ones(lp,1); + M = zeros(lp,1); + pr = p; + D = 1; M(1) = 1; + + old_p_index = 1; new_p_index = 2; + M_index = 1; pr_index = 2; + while(new_p_index<=lp) + if(abs(p(new_p_index) - p(old_p_index)) < toler) + # We've found a multiple pole. + M(M_index) = M(M_index) + 1; + e(new_p_index) = e(new_p_index-1) + 1; + # Remove the pole from pr. + pr(pr_index) = []; + else + # It's a different pole. + D++; M_index++; M(M_index) = 1; + old_p_index = new_p_index; pr_index++; + endif + new_p_index++; + endwhile + + # Shorten M to it's proper length + M = M(1:D); + + # Now set up the polynomial matrices. + + MM = max(M); + # Left hand side polynomial + lhs = zeros(MM,lb); + rhs = zeros(MM,lp*lp); + lhs(1,:) = b; + rhi = 1; + dpi = 1; + mpi = 1; + while(dpi<=D) + for ind = 1:M(dpi) + if(mpi>1 && (mpi+ind)<=lp) + cp = [p(1:mpi-1); p(mpi+ind:lp)]; + elseif (mpi==1) + cp = p(mpi+ind:lp); + else + cp = p(1:mpi-1); + endif + rhs(1,rhi:rhi+lp-1) = prepad(poly(cp),lp); + rhi = rhi + lp; + endfor + mpi = mpi + M(dpi); + dpi++; + endwhile + if(MM > 1) + for index = 2:MM + lhs(index,:) = prepad(polyderiv(lhs(index-1,:)),lb); + ind = 1; + for rhi = 1:lp + cp = rhs(index-1,ind:ind+lp-1); + rhs(index,ind:ind+lp-1) = prepad(polyderiv(cp),lp); + ind = ind + lp; + endfor + endfor + endif + + # Now lhs contains the numerator polynomial and as many derivatives as are + # required. rhs is a matrix of polynomials, the first row contains the + # corresponding polynomial for each residue and successive rows are + # derivatives. + + # Now we need to evaluate the first row of lhs and rhs at each distinct + # pole value. If there are multiple poles we will also need to evaluate + # the derivatives at the pole value also. + + B = zeros(lp,1); + A = zeros(lp,lp); + + dpi = 1; + row = 1; + while(dpi<=D) + for mi = 1:M(dpi) + B(row) = polyval(lhs(mi,:),pr(dpi)); + ci = 1; + for col = 1:lp + cp = rhs(mi,ci:ci+lp-1); + A(row,col) = polyval(cp,pr(dpi)); + ci = ci + lp; + endfor + row++; + endfor + dpi++; + endwhile + + # Solve for the residues. + r = A\B; + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/roots-amr.m @@ -0,0 +1,55 @@ +function r = roots(c) +#Find the roots of 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). +# +#roots(c) will return a vector r of length n such that +# +# p(x) = c(1) [ (x - r(1)) * (x - r(2)) * ... * (x - r(n)) ] +# +#roots and poly are inverse functions to within a scaling factor. +# +#SEE ALSO: poly, roots, conv, deconv, residue, filter, polyval, polyvalm, +# polyderiv, polyinteg + +# Author: +# Tony Richardson +# amr@mpl.ucsd.edu +# June 1994 + + if(nargin != 1) + error("usage: roots(vector)"); + endif + + if(is_matrix(c)) + error("argument must be a vector."); + endif + + n = length(c); + + if(is_scalar(c) || n == 0) + r = []; + return; + endif + + # Ensure that c is a row vector. + if(rows(c) > 1) + c = reshape(c,1,n); + endif + + # We could replace this with a call to compan, but it's faster to + # just reproduce the code here. + A = diag(ones(n-2,1),-1); + A(1,:) = -c(2:n)/c(1); + + r = eig(A); + + # Sort roots in order by decreasing magnitude. + [mr i] = sort(abs(r)); + r = r(i(length(i):-1:1)); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/set/complement.m @@ -0,0 +1,31 @@ +function y = complement(a,b) + +# usage: complement(a,b) +# +# Returns the elements of set b that are not in set a. +# +# See - create_set, union, intersection + + if(nargin != 2) + error("usage: complement(a,b)"); + endif + + if(isempty(a)) + y = create_set(b); + elseif(isempty(b)) + y = []; + else + a = create_set(a); + b = create_set(b); + yindex = 1; + y = zeros(1,length(b)); + for index = 1:length(b) + if(all(a != b(index))) + y(yindex++) = b(index); + endif + endfor + y = y(1:(yindex-1)); + endif + +endfunction +
new file mode 100644 --- /dev/null +++ b/scripts/set/create_set.m @@ -0,0 +1,36 @@ +function y = create_set(x) + +# usage: create_set(x) +# +# Returns the unique elements of x, sorted in ascending order. +# +# See - union, intersection, complement + + if ( nargin != 1) + error("usage: create_set(x)"); + endif + + if(isempty(x)) + y = []; + else + [nrx, ncx] = size(x); + nelx = nrx*ncx; + x = reshape(x,1,nelx); + y = zeros(1,nelx); + + x = sort(x); + cur_val = y(1) = x(1); + yindex = xindex = 2; + + while (xindex <= nelx) + if(cur_val != x(xindex)) + cur_val = x(xindex); + y(yindex++) = cur_val; + endif + xindex++; + endwhile + y = y(1:(yindex-1)); + endif + +endfunction +
new file mode 100644 --- /dev/null +++ b/scripts/set/intersection.m @@ -0,0 +1,41 @@ +function y = intersection(a,b) + +# usage: intersection(a,b) +# +# Returns the intersection of sets a and b. +# +# See - create_set, union, complement + + if (nargin != 2) + error("usage: intersection(a,b)"); + endif + + if(isempty(a) || isempty(b)) + y = []; + return; + endif + + a = create_set(a); + b = create_set(b); + + if(length(a) < length(b)) + yindex = 1; + y = zeros(1,length(a)); + for index = 1:length(a) + if(any(b == a(index))) + y(yindex++) = a(index); + endif + endfor + else + yindex = 1; + y = zeros(1,length(b)); + for index = 1:length(b) + if(any(a == b(index))) + y(yindex++) = b(index); + endif + endfor + endif + + y = y(1:(yindex-1)); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/set/union.m @@ -0,0 +1,25 @@ +function y = union(a,b) + +# usage: union(a,b) +# +# Returns the union of sets a and b. +# +# See - create_set, intersection, complement + + if (nargin != 2) + error("usage: union(a,b)"); + endif + + if(isempty(a)) + y = create_set(b); + elseif(isempty(b)) + y = create_set(a); + else + [nra, nca] = size(a); + a = reshape(a,1,nra*nca); + [nrb, ncb] = size(b); + b = reshape(b,1,nrb*ncb); + y = create_set([a, b]); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/signal/sinc.m @@ -0,0 +1,31 @@ +function result = sinc(x) + +# usage: sinc(x) +# +# Returns sin(pi*x)/(pi*x). +# + +# We either need to set the do_fortran_indexing variable to "true" +# or use reshape to convert the input matrix to a vector, so that +# we can use find to determine the elements of x that equal zero. +# I prefer reshaping. + + [nr, nc] = size(x); + + nels = nr*nc; + + x = reshape(x,nels,1); + + # Set result to all ones initially. + result = ones(nels,1); + + # Find non-zero elements in the input matrix. + i = find(x); + + if (!isempty(i)) + result(i) = sin(pi*x(i))./(pi*x(i)); + endif + + result = reshape(result,nr,nc); + +endfunction