freqz produces incorrect output
Frederick_Umminger at playstation.sony.com
Frederick_Umminger at playstation.sony.com
Wed Dec 24 15:21:08 CST 2008
I believe the following code fixes the problem. I apologize for it not
being a diff.
-Frederick
## Copyright (C) 1994, 1995, 1996, 1997, 1999, 2000, 2002, 2005, 2006,
## 2007 John W. Eaton
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING. If not, see
## <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{h}, @var{w}] =} freqz (@var{b},
@var{a}, @var{n}, "whole")
## Return the complex frequency response @var{h} of the rational IIR
filter
## whose numerator and denominator coefficients are @var{b} and @var{a},
## respectively. The response is evaluated at @var{n} angular frequencies
## between 0 and
## @ifinfo
## 2*pi.
## @end ifinfo
## @iftex
## @tex
## $2\pi$.
## @end tex
## @end iftex
##
## @noindent
## The output value @var{w} is a vector of the frequencies.
##
## If the fourth argument is omitted, the response is evaluated at
## frequencies between 0 and
## @ifinfo
## pi.
## @end ifinfo
## @iftex
## @tex
## $\pi$.
## @end tex
## @end iftex
##
## If @var{n} is omitted, a value of 512 is assumed.
##
## If @var{a} is omitted, the denominator is assumed to be 1 (this
## corresponds to a simple FIR filter).
##
## For fastest computation, @var{n} should factor into a small number of
## small primes.
##
## @deftypefnx {Function File} {@var{h} =} freqz (@var{b}, @var{a},
@var{w})
## Evaluate the response at the specific frequencies in the vector
@var{w}.
## The values for @var{w} are measured in radians.
##
## @deftypefnx {Function File} {[@dots{}] =} freqz (@dots{}, @var{Fs})
## Return frequencies in Hz instead of radians assuming a sampling rate
## @var{Fs}. If you are evaluating the response at specific frequencies
## @var{w}, those frequencies should be requested in Hz rather than
radians.
##
## @deftypefnx {Function File} {} freqz (@dots{})
## Plot the pass band, stop band and phase response of @var{h} rather
## than returning them.
## @end deftypefn
## Author: jwe ???
function [h_r, f_r] = freqz (b, a, n, region, Fs)
if (nargin < 1 || nargin > 5)
print_usage ();
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 (! ischar (region) && ! isempty (region))
Fs = region;
region = [];
endif
endif
if (isempty (b))
b = 1;
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
a = a(:);
b = b(:);
if (! isscalar (n)) ## Explicit frequency vector given
w = f = n;
if (nargin == 4) ## Sampling rate Fs was specified
w = 2*pi*f/Fs;
endif
k = max (length (b), length (a));
hb = polyval (postpad (b, k), exp (j*w));
ha = polyval (postpad (a, k), exp (j*w));
else
## polyval(fliplr(P),exp(jw)) is O(p n) and fft(x) is O(n log(n)),
## where p is the order of the polynomial P. For small p it
## would be faster to use polyval but in practice the overhead for
## polyval is much higher and the little bit of time saved isn't
## worth the extra code.
if (strcmp (region, "whole"))
N= n;
else
N = 2*n;
endif
f = Fs * (0:n-1)' / N;
sz = max(size(b,1), size(a,1));
pad_sz = N*ceil(sz/N);
b = postpad(b,pad_sz);
a = postpad(a,pad_sz);
hb = zeros(n,1);
ha = zeros(n,1);
for i = 1:N:pad_sz
hb = hb + fft (postpad (b(i:i+N-1), N))(1:n);
ha = ha + fft (postpad (a(i:i+N-1), N))(1:n);
endfor
endif
h = hb ./ ha;
if (nargout != 0), # return values and don't plot
h_r = h;
f_r = f;
else # plot and don't return values
freqz_plot (f, h);
endif
endfunction
%!test # correct values and fft-polyval consistency
%! # butterworth filter, order 2, cutoff pi/2 radians
%! b = [0.292893218813452 0.585786437626905 0.292893218813452];
%! a = [1 0 0.171572875253810];
%! [h,w] = freqz(b,a,32);
%! assert(h(1),1,10*eps);
%! assert(abs(h(17)).^2,0.5,10*eps);
%! assert(h,freqz(b,a,w),10*eps); # fft should be consistent with polyval
%!test # whole-half consistency
%! b = [1 1 1]/3; # 3-sample average
%! [h,w] = freqz(b,1,32,'whole');
%! assert(h(2:16),conj(h(32:-1:18)),20*eps);
%! [h2,w2] = freqz(b,1,16,'half');
%! assert(h(1:16),h2,20*eps);
%! assert(w(1:16),w2,20*eps);
%!test # Sampling frequency properly interpreted
%! b = [1 1 1]/3; a = [1 0.2];
%! [h,f] = freqz(b,a,16,320);
%! assert(f,[0:15]'*10,10*eps);
%! [h2,f2] = freqz(b,a,[0:15]*10,320);
%! assert(f2,[0:15]*10,10*eps);
%! assert(h,h2.',20*eps);
%! [h3,f3] = freqz(b,a,32,'whole',320);
%! assert(f3,[0:31]'*10,10*eps);
Frederick Umminger/R&D/SCEA
12/24/2008 12:53 PM
To
bug at octave.org
cc
Frederick Umminger/R&D/SCEA at Playstation
Subject
freqz produces incorrect output
Bug report for Octave 3.0.3 configured for i686-pc-msdosmsvc
Description:
-----------
The command [h,w] = freqz(b) produces an incorrect answer if the vector
b has length greater than 512.
Similarly for [h,w] = freqz(b,a) , [h,w] = freqz(b,a,n) if the length of
b or a is greater than 512, n respectively.
Repeat-By:
---------
% Compute the spectrum of an audio file
[audio, fs] = auload("InputSample.wav");
[h_audio, w_audio] = freqz(audio );
Fix:
---
This is because the branch of the code that uses the fft as a
replacement for polyval has an implicit assumption about relative sizes.
If this assumption is violated, then the call to postpad truncates b and
a.
In order to be correct, this code must be rewritten to sum over multiple
ffts when either b or a is larger than n.
Configuration (please do not edit this section):
-----------------------------------------------
uname output: Windows
configure opts:
Fortran compiler: fc-msvc
FFLAGS: -O2
F2C: @F2C@
F2CFLAGS: @F2CFLAGS@
FLIBS: -lhdf5 -lzlib -lf2c -lkernel32
CPPFLAGS: -I. -Ic:/Software/VCLibs/include
INCFLAGS: -I. -I. -I./liboctave -I./src -I./libcruft/misc
C compiler: cc-msvc, version
CFLAGS: -O2 -MD
CPICFLAG:
C++ compiler: cc-msvc, version
CXXFLAGS: -O2 -EHs -MD
CXXPICFLAG:
LD_CXX: cc-msvc
LDFLAGS:
LIBFLAGS: -L.
RLD_FLAG:
BLAS_LIBS: -llapack -lblas
FFTW_LIBS: -lfftw3
LIBS: -lreadline -lncurses -lhdf5 -lzlib -lws2_32 -lkernel32
LEXLIB:
LIBGLOB: -lglob
SED: /usr/bin/sed
DEFS:
-DPACKAGE_NAME="" -DPACKAGE_TARNAME="" -DPACKAGE_VERSION=""
-DPACKAGE_STRING="" -DPACKAGE_BUGREPORT="" -DOCTAVE_SOURCE=1
-D_GNU_SOURCE=1 -DSTDC_HEADERS=1 -DHAVE_SYS_TYPES_H=1 -DHAVE_SYS_STAT_H=1
-DHAVE_STDLIB_H=1 -DHAVE_STRING_H=1 -DHAVE_MEMORY_H=1 -DHAVE_INTTYPES_H=1
-DHAVE_STDINT_H=1 -DHAVE_UNISTD_H=1 -DSEPCHAR=';' -DSEPCHAR_STR=";"
-D__NO_MATH_INLINES=1 -DCXX_NEW_FRIEND_TEMPLATE_DECL=1
-DCXX_ISO_COMPLIANT_LIBRARY=1 -DCXX_ABI=unknown -DHAVE_QHULL=1
-DHAVE_PCRE=1 -DHAVE_ZLIB_H=1 -DHAVE_ZLIB=1 -DHAVE_HDF5_H=1 -DHAVE_HDF5=1
-DHAVE_H5GGET_NUM_OBJS=1 -D_HDF5USEDLL_=1 -DHAVE_FFTW3=1 -DHAVE_GLPK_H=1
-DHAVE_GLPK=1 -DHAVE_CURL_CURL_H=1 -DHAVE_CURL=1
-DHAVE_IEEE754_DATA_FORMAT=1 -DF77_FUNC(name,NAME)=name ## _
-DF77_FUNC_(name,NAME)=name ## __ -DHAVE_BLAS=1
-DHAVE_SUITESPARSE_UMFPACK_H=1 -DHAVE_UMFPACK=1 -DUMFPACK_SEPARATE_SPLIT=1
-DHAVE_SUITESPARSE_COLAMD_H=1 -DHAVE_COLAMD=1
-DHAVE_SUITESPARSE_CCOLAMD_H=1 -DHAVE_CCOLAMD=1
-DHAVE_SUITESPARSE_CHOLMOD_H=1 -DHAVE_CHOLMOD=1 -DHAVE_SUITESPARSE_CS_H=1
-DHAVE_CXSPARSE=1 -Dmode_t=int -Dpid_t=int -Duid_t=int -Dgid_t=int
-DHAVE_DEV_T=1 -DHAVE_INO_T=1 -DHAVE_LONG_LONG_INT=1
-DHAVE_UNSIGNED_LONG_LONG_INT=1 -DHAVE_SIG_ATOMIC_T=1 -DSIZEOF_SHORT=2
-DSIZEOF_INT=4 -DSIZEOF_LONG=4 -DSIZEOF_LONG_LONG=8 -DHAVE_ALLOCA=1
-DNPOS=std::string::npos -DHAVE_PLACEMENT_DELETE=1 -DSTDC_HEADERS=1
-DHAVE_ASSERT_H=1 -DHAVE_DIRECT_H=1 -DHAVE_FCNTL_H=1 -DHAVE_FLOAT_H=1
-DHAVE_INTTYPES_H=1 -DHAVE_LIMITS_H=1 -DHAVE_LOCALE_H=1 -DHAVE_MEMORY_H=1
-DHAVE_STDINT_H=1 -DHAVE_STDLIB_H=1 -DHAVE_STRING_H=1 -DHAVE_SYS_STAT_H=1
-DHAVE_SYS_TYPES_H=1 -DHAVE_SYS_UTIME_H=1 -DHAVE_UNISTD_H=1
-DHAVE_VARARGS_H=1 -DHAVE_SSTREAM=1 -DHAVE_GLOB_H=1 -DHAVE_FNMATCH_H=1
-DHAVE_CONIO_H=1 -DHAVE_ATEXIT=1 -DHAVE_CHMOD=1 -DHAVE_DUP2=1
-DHAVE_EXECVP=1 -DHAVE_GETCWD=1 -DHAVE_GETPID=1 -DHAVE__KBHIT=1
-DHAVE_MEMMOVE=1 -DHAVE_MKDIR=1 -DHAVE_PUTENV=1 -DHAVE_RAISE=1
-DHAVE_RENAME=1 -DHAVE_RMDIR=1 -DHAVE_SETLOCALE=1 -DHAVE_SETVBUF=1
-DHAVE_STAT=1 -DHAVE_STRDUP=1 -DHAVE_STRERROR=1 -DHAVE_STRICMP=1
-DHAVE_STRNICMP=1 -DHAVE_TEMPNAM=1 -DHAVE_UMASK=1 -DHAVE_UNLINK=1
-DHAVE_UTIME=1 -DHAVE_VFPRINTF=1 -DHAVE_VSPRINTF=1 -DHAVE_VSNPRINTF=1
-DHAVE__CHMOD=1 -DHAVE__SNPRINTF=1 -DHAVE__UTIME32=1
-DOCTAVE_HAVE_BROKEN_STRPTIME=1 -D_USE_MATH_DEFINES=1
-DHAVE_LOADLIBRARY_API=1 -DENABLE_DYNAMIC_LINKING=1 -DHAVE__FINITE=1
-DHAVE__ISNAN=1 -DHAVE__COPYSIGN=1 -DHAVE_DECL_SIGNBIT=0
-DHAVE_STRUCT_STAT_ST_RDEV=1 -DHAVE_DECL_TZNAME=1 -DHAVE_TZNAME=1
-DCLOSEDIR_VOID=1 -DMKDIR_TAKES_ONE_ARG=1 -DUSE_READLINE=1
-DRETSIGTYPE=void -DHAVE_DECL_SYS_SIGLIST=0 -DMUST_REINSTALL_SIGHANDLERS=1
-DRETSIGTYPE_IS_VOID=1 -DGNUPLOT_BINARY="pgnuplot"
User-preferences (please do not edit this section):
EDITOR = C:\Program Files\Octave\tools\wscite\SciTE.exe
EXEC_PATH = C:\Program Files\Octave\msys\bin;C:\Program
Files\Octave\libexec\octave\3.0.3\site\exec\i686-pc-msdosmsvc;C:\Program
Files\Octave\libexec\octave\api-v32\site\exec\i686-pc-msdosmsvc;C:\Program
Files\Octave\libexec\octave\site\exec\i686-pc-msdosmsvc;C:\Program
Files\Octave\libexec\octave\3.0.3\exec\i686-pc-msdosmsvc;C:\Program
Files\Octave\bin;C:\Program
Files\Octave;C:\usr\local\cell\\host-win32\bin;C:\usr\local\cell\\host-win32\Cg\bin;C:\Program
Files\SN Systems\Common\VSI\bin;C:\Program Files\SN
Systems\PS3\bin;C:\Program Files\SN
Systems\Common\bin;C:\usr\local\cell\\host-win32\sn\bin;C:\usr\local\cell\\host-win32\spu\bin;C:\usr\local\cell\\host-win32\ppu\bin;C:\Perl\site\bin;C:\Perl\bin;C:\WINDOWS\system32;C:\WINDOWS;C:\WINDOWS\System32\Wbem;C:\Program
Files\Intel\DMIX;C:\Program Files\Common Files\Roxio
Shared\DLLShared\;c:\Program Files\Microsoft SQL
Server\90\Tools\binn\;C:\Program Files\CollabNet Subversion
Server;C:\Program Files\QuickTime\QTSystem\;C:\Program
Files\TortoiseSVN\bin;C:\Program Files\sox-14.2.0
IMAGE_PATH = .;C:\Program Files\Octave\share\octave\3.0.3\imagelib
PAGER = less
PS1 = \s:\#>
PS2 = >
PS4 = +
beep_on_error = 0
completion_append_char =
crash_dumps_octave_core = 1
echo_executing_commands = 0
fixed_point_format = 0
gnuplot_binary = pgnuplot
gnuplot_command_end =
gnuplot_command_plot = pl
gnuplot_command_replot = rep
gnuplot_command_splot = sp
gnuplot_command_title = t
gnuplot_command_using = u
gnuplot_command_with = w
history_file = H:\.octave_hist
history_size = 1024
ignore_function_time_stamp = system
info_file = C:\Program Files\Octave\share\info\octave.info
info_program = info
makeinfo_program = makeinfo
max_recursion_depth = 256
output_max_field_width = 5
output_precision = 5
page_output_immediately = 0
page_screen_output = 1
print_answer_id_name = 1
print_empty_dimensions = 1
save_precision = 16
saving_history = 1
sighup_dumps_octave_core = 1
sigterm_dumps_octave_core = 1
silent_functions = 0
split_long_rows = 1
string_fill_char =
struct_levels_to_print = 2
suppress_verbose_help_message = 0
-------------- next part --------------
An HTML attachment was scrubbed...
URL: https://www-old.cae.wisc.edu/pipermail/bug-octave/attachments/20081224/e57d6922/attachment-0001.html
More information about the Bug-octave
mailing list