From bug-octave-request at bevo dot che dot wisc dot edu Thu Nov 16 17:48:02 2000 Subject: freqz extensions From: Paul Kienzle To: bug-octave at bevo dot che dot wisc dot edu, pkienzle@kienzle.powernet.co.uk Date: Thu, 16 Nov 2000 22:43:46 +0000 To: bug-octave at bevo dot che dot wisc dot edu Cc: pkienzle Subject: freqz extensions Bug report for Octave 2.1.31 configured for %OCTAVE_CANONICAL_HOST_TYPE% Description: ----------- I've made a number of updates to the freqz function which I find useful. The following are compatible with Matlab's freqz: - if nargout=0, produce mag/phase plot and don't return values - add sampling frequency Fs as optional fifth parameter - evaluate at specific frequencies if a vector of frequencies is given instead of n. I don't know what matlab's behaviour is for complex filters, but displaying the whole frequency response by default makes sense to me: - default to 'whole' if complex b or a Other updates: - allow [] for default paramter values - print usage if not enough arguments - minor optimization: don't compute fft if the poly has no roots I haven't implemented the following undocumented features: freqz(...,'magnitude','linear') % linear rather than logarithmic mag freqz(...,'phase','no') % just plot the magnitude graph I don't know what other properties are associated with freqz graphs. Repeat-By: --------- Fix: --- *** ChangeLog 2000/11/16 01:52:25 1.1 --- ChangeLog 2000/11/16 01:55:44 *************** *** 1,3 **** --- 1,12 ---- + 2000-11-15 Paul Kienzle * control/base/dkalman.m: New file. *** signal/freqz.m 2000/11/16 01:03:38 1.1 --- signal/freqz.m 2000/11/16 22:36:25 *************** *** 53,107 **** ## ## For fastest computation, at var{n} should factor into a small number of ## small primes. ## at end deftypefn ## Author: jwe ??? ! function [h, w] = freqz(b,...) ! if (nargin == 1) ## Response of an FIR filter. ! a = 1; ! n = 512; ! region = "half"; elseif (nargin == 2) ## Response of an IIR filter ! a = va_arg(); ! n = 512; ! region = "half"; elseif (nargin == 3) ! a = va_arg(); ! n = va_arg(); ! region = "half"; elseif (nargin == 4) ! a = va_arg(); ! n = va_arg(); ! region = va_arg(); endif la = length(a); a = reshape(a,1,la); lb = length(b); b = reshape(b,1,lb); - k = max([la, lb]); ! if (n >= k) ! if (strcmp(region,"whole")) ! h = fft(postpad(b,n)) ./ fft(postpad(a,n)); ! w = 2*pi*[0:(n-1)]/n; else ! h = fft(postpad(b,2*n)) ./ fft(postpad(a,2*n)); ! h = h(1:n); ! w = pi*[0:(n-1)]/n; endif else ! if (strcmp(region,"whole")) ! w = 2*pi*[0:(n-1)]/n; else ! w = pi*[0:(n-1)]/n; endif ! h = polyval(postpad(b,k),exp(j*w)) ./ polyval(postpad(a,k),exp(j*w)); endif endfunction --- 53,212 ---- ## ## For fastest computation, at var{n} should factor into a small number of ## small primes. + ## + ## at deftypefnx {Function File} {@var{h} =} freqz (@var{b}, @var{a}, @var{w}) + ## Evaluate the response at the specific frequencies in the vector at {w} dot + ## + ## at deftypefnx {Function File} {... =} freqz (..., @var{Fs}) + ## Return frequencies in Hz instead of radians assuming a sampling rate + ## at {Fs} dot If you are evaluating the response at specific frequencies + ## at var{w}, those frequencies should be requested in Hz rather than radians. + ## + ## at deftypefnx {Function File} freqz(...) + ## Plot the pass band, stop band and phase response of at var{h} rather + ## than returning them. + ## ## at end deftypefn ## Author: jwe ??? ! function [h_r, w_r] = freqz(b, a, n, region, Fs) ! if (nargin<1 || nargin>5) ! usage("[h w]=freqz(b, a, n [, 'whole'] [, Fs])"); ! elseif (nargin == 1) ## Response of an FIR filter. ! a=[]; n=[]; region=[]; Fs=[]; elseif (nargin == 2) ## Response of an IIR filter ! n=[]; region=[]; Fs=[]; elseif (nargin == 3) ! region=[]; Fs=[]; elseif (nargin == 4) ! Fs=[]; ! if !isstr(region) && !isempty(region) ! Fs = region; region=[]; ! endif ! endif ! ! if isempty(a) a=1; endif ! if isempty(n) n=512; endif ! if isempty(region) ! if isreal(b) && isreal(a) ! region = "half"; ! else ! region = "whole"; ! endif ! endif ! if isempty(Fs) ! if (nargout==0) Fs = 2; else Fs = 2*pi; endif endif la = length(a); a = reshape(a,1,la); lb = length(b); b = reshape(b,1,lb); k = max([la, lb]); ! if !is_scalar(n) ! if nargin==4 ## Fs was specified ! w = 2*pi*n/Fs; else ! w = n; endif + n = length(n); + extent = 0; + elseif strcmp(region,"whole") + w = 2*pi*[0:(n-1)]/n; + extent = n; else ! w = pi*[0:(n-1)]/n; ! extent = 2*n; ! endif ! ! if length(b) == 1 ! if length(a) == 1 ! hb=b*ones(1,n); else ! hb = b; endif ! elseif extent >= k ! hb = fft(postpad(b,extent)); ! hb = hb(1:n); ! else ! hb = polyval(postpad(b,k),exp(j*w)); endif + if length(a) == 1 + ha = a; + elseif extent >= k + ha = fft(postpad(a,extent)); + ha = ha(1:n); + else + ha = polyval(postpad(a,k),exp(j*w)); + endif + h = hb./ha; + w = Fs*w/(2*pi); + + if nargout != 0, # return values and don't plot + h_r = h; + w_r = w; + else # plot and don't return values + ## ## exclude zero-frequency + ## h = h(2:length(h)); + ## w = w(2:length(w)); + ## n = n-1; + mag = 20*log10(abs(h)); + phase = unwrap(arg(h)); + maxmag = max(mag); + + unwind_protect # protect graph state + if gnuplot_has_multiplot + subplot(311); + gset lmargin 10; + axis("labely"); + xlabel(""); + endif + grid("on"); + axis([w(1), w(n), maxmag-3, maxmag]); + plot(w, mag, ";Pass band (dB);"); + + if gnuplot_has_multiplot + subplot(312); + axis("labely"); + title(""); + xlabel(""); + gset tmargin 0; + else + input("press any key for the next plot: "); + endif + grid("on"); + if (maxmag - min(mag) > 100) + axis([w(1), w(n), maxmag-100, maxmag]); + else + axis("autoy"); + endif + plot(w, mag, ";Stop band (dB);"); + + if gnuplot_has_multiplot + subplot(313); + axis("label"); + title(""); + else + input("press any key for the next plot: "); + endif + grid("on"); + axis("autoy"); + xlabel(["Frequency (Fs=", num2str(Fs), ")"]); + axis([w(1), w(n)]); + plot(w, phase/(2*pi), ";Phase (radians/2pi);"); + + unwind_protect_cleanup # restore graph state + grid("off"); + axis("auto","label"); + gset lmargin; + gset tmargin; + oneplot(); + end_unwind_protect + end endfunction Configuration (please do not edit this section): ----------------------------------------------- uname output: Linux kienzle 2.2.17 #1 Sun Sep 10 17:51:52 BST 2000 i586 unknown configure opts: --prefix=/usr --datadir=/usr/share --libdir=/usr/lib --libexecdir=/usr/lib --infodir=/usr/share/info --mandir=/usr/share/man --with-g77 --with-fastblas --enable-dl --enable-shared --enable-lite-kernel --disable-static --host i386-linux Fortran compiler: g77 FFLAGS: -O2 F2C: F2CFLAGS: FLIBS: -lg2c -lm -L/usr/lib/gcc-lib/i386-linux/2.95.2 -lm CPPFLAGS: INCFLAGS: -I. -I. -I./liboctave -I./src -I./libcruft/misc -I./glob -I./glob C compiler: gcc, version 2.95.2 20000220 (Debian GNU/Linux) CFLAGS: -O2 CPICFLAG: -fPIC C++ compiler: c++, version 2.95.2 20000220 (Debian GNU/Linux) CXXFLAGS: -O2 CXXPICFLAG: -fPIC LDFLAGS: -s LIBFLAGS: -L. RLD_FLAG: -Xlinker -rpath -Xlinker /usr/lib/octave-2.1.31 TERMLIBS: -lncurses LIBS: LEXLIB: LIBPLPLOT: LIBDLFCN: LIBGLOB: ./glob/glob.o ./glob/fnmatch.o DEFS: -DOCTAVE_SOURCE=1 -DSEPCHAR=':' -DSEPCHAR_STR=":" -DUSE_READLINE=1 -D__NO_MATH_INLINES=1 -DCXX_NEW_FRIEND_TEMPLATE_DECL=1 -DHAVE_LIBM=1 -DHAVE_LIBZ=1 -DF77_APPEND_UNDERSCORE=1 -DOCTAVE_LITE=1 -DSIZEOF_SHORT=2 -DSIZEOF_INT=4 -DSIZEOF_LONG=4 -DSIZEOF_LONG_LONG=8 -DHAVE_ALLOCA_H=1 -DHAVE_ALLOCA=1 -DNPOS=std::string::npos -DSTDC_HEADERS=1 -DHAVE_DIRENT_H=1 -DTIME_WITH_SYS_TIME=1 -DHAVE_SYS_WAIT_H=1 -DHAVE_ASSERT_H=1 -DHAVE_CURSES_H=1 -DHAVE_DLFCN_H=1 -DHAVE_FCNTL_H=1 -DHAVE_FLOAT_H=1 -DHAVE_FNMATCH_H=1 -DHAVE_GLOB_H=1 -DHAVE_GRP_H=1 -DHAVE_LIMITS_H=1 -DHAVE_MEMORY_H=1 -DHAVE_NCURSES_H=1 -DHAVE_POLL_H=1 -DHAVE_PWD_H=1 -DHAVE_SGTTY_H=1 -DHAVE_STDLIB_H=1 -DHAVE_STRING_H=1 -DHAVE_SYS_IOCTL_H=1 -DHAVE_SYS_PARAM_H=1 -DHAVE_SYS_POLL_H=1 -DHAVE_SYS_RESOURCE_H=1 -DHAVE_SYS_SELECT_H=1 -DHAVE_SYS_STAT_H=1 -DHAVE_SYS_TIME_H=1 -DHAVE_SYS_TIMES_H=1 -DHAVE_SYS_TYPES_H=1 -DHAVE_SYS_UTSNAME_H=1 -DHAVE_TERMCAP_H=1 -DHAVE_TERMIO_H=1 -DHAVE_UNISTD_H=1 -DHAVE_VARARGS_H=1 -DHAVE_ATEXIT=1 -DHAVE_BCOPY=1 -DHAVE_BZERO=1 -DHAVE_DUP2=1 -DHAVE_ENDGRENT=1 -DHAVE_ENDPWENT=1 -DHAVE_EXECVP=1 -DHAVE_FCNTL=1 -DHAVE_FORK=1 -DHAVE_GETCWD=1 -DHAVE_GETEGID=1 -DHAVE_GETEUID=1 -DHAVE_GETGID=1 -DHAVE_GETGRENT=1 -DHAVE_GETGRGID=1 -DHAVE_GETGRNAM=1 -DHAVE_GETHOSTNAME=1 -DHAVE_GETPGRP=1 -DHAVE_GETPID=1 -DHAVE_GETPPID=1 -DHAVE_GETPWENT=1 -DHAVE_GETPWNAM=1 -DHAVE_GETPWUID=1 -DHAVE_GETTIMEOFDAY=1 -DHAVE_GETUID=1 -DHAVE_GETWD=1 -DHAVE_LOCALTIME_R=1 -DHAVE_LSTAT=1 -DHAVE_MEMMOVE=1 -DHAVE_MKDIR=1 -DHAVE_MKFIFO=1 -DHAVE_ON_EXIT=1 -DHAVE_PIPE=1 -DHAVE_POLL=1 -DHAVE_PUTENV=1 -DHAVE_RENAME=1 -DHAVE_RINDEX=1 -DHAVE_RMDIR=1 -DHAVE_SELECT=1 -DHAVE_SETGRENT=1 -DHAVE_SETPWENT=1 -DHAVE_SETVBUF=1 -DHAVE_SIGACTION=1 -DHAVE_SIGPENDING=1 -DHAVE_SIGPROCMASK=1 -DHAVE_SIGSUSPEND=1 -DHAVE_STAT=1 -DHAVE_STRCASECMP=1 -DHAVE_STRDUP=1 -DHAVE_STRERROR=1 -DHAVE_STRFTIME=1 -DHAVE_STRNCASECMP=1 -DHAVE_STRPTIME=1 -DHAVE_TEMPNAM=1 -DHAVE_UMASK=1 -DHAVE_UNLINK=1 -DHAVE_USLEEP=1 -DHAVE_VFPRINTF=1 -DHAVE_VSPRINTF=1 -DHAVE_VSNPRINTF=1 -DHAVE_WAITPID=1 -DHAVE_LIBDL=1 -DHAVE_DLOPEN=1 -DHAVE_DLSYM=1 -DHAVE_DLERROR=1 -DHAVE_DLCLOSE=1 -DWITH_DL=1 -DWITH_DYNAMIC_LINKING=1 -DHAVE_TIMEVAL=1 -DHAVE_FINITE=1 -DHAVE_ISNAN=1 -DHAVE_ISINF=1 -DHAVE_ACOSH=1 -DHAVE_ASINH=1 -DHAVE_ATANH=1 -DHAVE_ERF=1 -DHAVE_ERFC=1 -DHAVE_ST_BLKSIZE=1 -DHAVE_ST_BLOCKS=1 -DHAVE_ST_RDEV=1 -DHAVE_TM_ZONE=1 -DHAVE_GR_PASSWD=1 -DEXCEPTION_IN_MATH=1 -DRETSIGTYPE=void -DSYS_SIGLIST_DECLARED=1 -DHAVE_SYS_SIGLIST=1 -DHAVE_POSIX_SIGNALS=1 -DHAVE_GETRUSAGE=1 -DHAVE_TIMES=1 -DGNUPLOT_HAS_MULTIPLOT=1 -DGNUPLOT_HAS_FRAMES=1 User-preferences (please do not edit this section): -------------------------------------------------- EDITOR = "vi" EXEC_PATH = ":/scratch/programs/audio/package/pipewave.1.3/bin:/scratch/programs/audio/pitch/POWERpv1.1G/bin:/home/pkienzle/pantome/bin:/home/pkienzle/sfs/bin:/home/pkienzle/bin:/usr/local/bin:/usr/bin:/bin:/usr/X11R6/bin:/usr/games" IMAGEPATH = ".:/usr/share/octave/2.1.31/imagelib//" INFO_FILE = "/usr/share/info/octave.info" INFO_PROGRAM = "info" LOADPATH = "/home/pkienzle/matcompat21//:/home/pkienzle/octave//::" PAGER = "less" PS1 = "\\s:\\#> " PS2 = "> " PS4 = "+ " automatic_replot = 0 beep_on_error = 0 completion_append_char = " " default_eval_print_flag = 1 # default_global_variable_value = default_return_value = [] default_save_format = "ascii" define_all_return_values = 0 do_fortran_indexing = 0 echo_executing_commands = 0 empty_list_elements_ok = "warn" fixed_point_format = 0 gnuplot_binary = "gnuplot" gnuplot_command_end = "\n" gnuplot_command_plot = "pl" gnuplot_command_replot = "rep" gnuplot_command_splot = "sp" gnuplot_command_title = "t" gnuplot_command_using = "u" gnuplot_command_with = "w" gnuplot_has_frames = 1 gnuplot_has_multiplot = 1 history_file = "/home/pkienzle/.octave_hist" history_size = 1024 ignore_function_time_stamp = "system" implicit_num_to_str_ok = 0 implicit_str_to_num_ok = 0 initialize_global_variables = 0 max_recursion_depth = 256 ok_to_lose_imaginary_part = "warn" output_max_field_width = 10 output_precision = 5 page_output_immediately = 0 page_screen_output = 1 prefer_column_vectors = 1 print_answer_id_name = 1 print_empty_dimensions = 1 print_rhs_assign_val = 0 propagate_empty_matrices = 1 resize_on_range_error = 1 return_last_computed_value = 0 save_precision = 15 saving_history = 1 silent_functions = 0 split_long_rows = 1 string_fill_char = " " struct_levels_to_print = 2 suppress_verbose_help_message = 1 treat_neg_dim_as_zero = 0 warn_assign_as_truth_value = 1 warn_divide_by_zero = 1 warn_function_name_clash = 1 warn_future_time_stamp = 1 warn_missing_semicolon = 0 warn_variable_switch_label = 0 whitespace_in_literal_matrix = ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web: http://www.octave.org How to fund new projects: http://www.octave.org/funding.html Subscription information: http://www.octave.org/archive.html -------------------------------------------------------------