ref: 29b0b173de4e3a9841db98526b8621d9bb88f0d8
parent: d1054fbde6bf7bbaf5d24af9a80ad7a8cbe261db
author: Erik de Castro Lopo <[email protected]>
date: Sun Apr 21 21:09:07 EDT 2019
Add Octave scripts to generate filter coefficients Closes: https://github.com/erikd/libsamplerate/issues/52
--- a/Makefile.am
+++ b/Makefile.am
@@ -4,7 +4,9 @@
AM_CPPFLAGS = -I$(top_srcdir)/src $(OS_SPECIFIC_INCLUDES)
-EXTRA_DIST = autogen.sh libsamplerate.spec.in samplerate.pc.in Make.bat Win32
+EXTRA_DIST = autogen.sh libsamplerate.spec.in samplerate.pc.in Make.bat Win32 \
+ Octave/generate_filter.m Octave/make_filter.m Octave/measure_filter.m \
+ Octave/Readme.md
pkgconfig_DATA = samplerate.pc
--- /dev/null
+++ b/Octave/Readme.md
@@ -1,0 +1,27 @@
+# Gneerating Filter Coefficients
+
+The scripts below are used in [GNU Octave][GNU Octave] which should be available
+on most platforms including Linux, Mac and Windows.
+
+On Debian and other Debian derived Linux distributions, Octave can be installed
+using: `sudo apt install octave octave-signal`.
+
+Filter coefficients can be generated as follows by running the `octave` command
+in the same directory as the scripts and this Readme file:
+
+```
+octave:1> pkg load signal
+octave:2> f = make_filter (8, 128, 100.3) ;
+# f = make_fip_filter (8, 128, 100.3) ;
+# Coeff. count : 4922
+# Fudge factor : 1.2018769
+# Pass band width : 0.0039062499 (should be 0.0039062500)
+# Stop band atten. : 100.71 dB
+# -3dB band Width : 0.490
+# half length : 2464
+# increment : 128
+```
+
+The parameters above generate the coefficient set used in `src/fastest_coeffs.h`.
+
+[GNU Octave]: https://www.gnu.org/software/octave/
--- /dev/null
+++ b/Octave/generate_filter.m
@@ -1,0 +1,35 @@
+function f = generate_filter (cycles, fudge_factor, increment, atten)
+
+if nargin != 4,
+ error ("Need four args.") ;+ endif
+
+# Calclate N and make sure it is even.
+N = fix (4 * cycles * fudge_factor * increment) ;
+
+if rem (N, 2) != 0,
+ N = N - 1 ;
+ endif
+
+# Generate the Sinc function.
+
+m = -((N-1)/2):((N-1)/2) ;
+f = sinc (m / fudge_factor / increment) ;
+
+# Genertate the window function and apply it.
+
+w = kaiser (N, (atten + 0.5) / 10) ;
+w = w' ;
+
+f = f .* w ;
+
+f = f / sum (f) ;
+
+endfunction
+
+# Do not edit or modify anything in this comment block.
+# The arch-tag line is a file identity tag for the GNU Arch
+# revision control system.
+#
+# arch-tag: 7e57a3cb-3f5c-4346-bfcd-4da1e758e2a7
+
--- /dev/null
+++ b/Octave/make_filter.m
@@ -1,0 +1,134 @@
+function f = make_filter (cycles, increment, atten, filename)
+
+# This works :
+#
+# f = make_filter (67, 128, 100.3) ;
+# f = make_filter (13, 128, 100.5) ;
+# f = make_filter (185, 4, 157.0) ;
+
+
+#=======================================================================
+
+if nargin < 3,
+ error ('Try make_filter (12, 32, 88, "output.txt")') ;+ endif
+
+if nargin < 4,
+ filename = 0 ;
+elseif (isstr (filename) == 0),
+ error ("Fourth parameter must be a file name.") ;+ endif
+
+fudge_factor1 = 1.0 ;
+f1 = generate_filter (cycles, fudge_factor1, increment, atten) ;
+[stop_atten stop_band_start1 minus_3db] = measure_filter (f1, atten) ;
+printf (" fudge_factor : %15.13f stop_band_start : %15.13f 1\n", fudge_factor1, stop_band_start1) ;+
+fudge_factor2 = 1.25 ;
+f2 = generate_filter (cycles, fudge_factor2, increment, atten) ;
+[stop_atten stop_band_start2 minus_3db] = measure_filter (f2, atten) ;
+printf (" fudge_factor : %15.13f stop_band_start : %15.13f 2\n", fudge_factor2, stop_band_start2) ;+
+f = f1 ;
+fudge_factor = fudge_factor1 ;
+stop_band_start = stop_band_start1 ;
+
+while ((stop_band_start1 - stop_band_start2) > 0.0000000001)
+ if (stop_band_start1 < stop_band_start2)
+ printf ("stop_band_start1 < stop_band_start2\n") ;+ break ;
+ endif
+
+ fudge_factor = fudge_factor1 + (fudge_factor2 - fudge_factor1) / 2 ;
+ f = generate_filter (cycles, fudge_factor, increment, atten) ;
+ [stop_atten stop_band_start minus_3db] = measure_filter (f, atten) ;
+
+ if (stop_band_start > 1.0)
+ printf ("A %10.8f %10.8f %10.8f\n", fudge_factor1, fudge_factor, fudge_factor2) ;+ continue ;
+ endif
+
+ if (stop_band_start < 0.5 / increment)
+ f2 = f ;
+ stop_band_start2 = stop_band_start ;
+ fudge_factor2 = fudge_factor ;
+ choice = 2 ;
+ else
+ f1 = f ;
+ stop_band_start1 = stop_band_start ;
+ fudge_factor1 = fudge_factor ;
+ choice = 1 ;
+ endif
+
+ printf (" fudge_factor : %15.13f stop_band_start : %15.13f %d\n", fudge_factor, stop_band_start, choice) ;+ endwhile
+
+printf ("\n") ;+
+#-------------------------------------------------------------------------------
+# Grab only half the coefficients.
+
+N = length (f) ;
+
+f = increment * f' ;
+
+if rem (length (f), 2) == 0,
+ index = find (f == max (f)) ;
+ index = min (index) - 1 ;
+ half_f = f (index:length (f)) ;
+else
+ error ("Length should be even.") ;+ endif
+
+trailing_zeros = 4 - rem (length (half_f), 4) ;
+
+#-------------------------------------------------------------------------------
+# Print analysis.
+
+printf ("# f = make_filter (%d, %d, %4.1f) ;\n", cycles, increment, atten) ;+printf ("# Coeff. count : %d\n", N) ;+printf ("# Fudge factor : %9.7f\n", fudge_factor) ;+printf ("# Pass band width : %12.10f (should be %12.10f)\n", stop_band_start, 0.5 / increment) ;+printf ("# Stop band atten. : %5.2f dB\n", abs (stop_atten)) ;+printf ("# -3dB band Width : %5.3f\n", 0.5 / increment / minus_3db) ;+printf ("# half length : %d\n", length (half_f) + trailing_zeros) ;+printf ("# increment : %d\n", increment) ;+
+if filename,
+ file = fopen (filename, "w") ;
+ if file == 0,
+ str = sprintf ("Error, not able to open '%s'", filename)+ error (str) ;
+ endif
+
+ fprintf (file, "/*\n") ;
+ fprintf (file, "** f = make_filter (%d, %d, %4.1f) ;\n", cycles, increment, atten) ;
+ fprintf (file, "** Pass band width : %9.7f (should be %9.7f)\n", stop_band_start, 0.5 / increment) ;
+ fprintf (file, "** Stop band atten. : %5.2f dB\n", abs (stop_atten)) ;
+ fprintf (file, "** -3dB band width : %5.3f\n", 0.5 / increment / minus_3db) ;
+ fprintf (file, "** half length : %d\n", length (half_f)) ;
+ fprintf (file, "** increment : %d\n", increment) ;
+ fprintf (file, "*/\n\n") ;
+
+ for val = half_f,
+ fprintf (file, "% 24.20e,\n", val) ;
+ endfor
+
+ if trailing_zeros > 0,
+ for val = 2:trailing_zeros,
+ fprintf (file, " 0,\n") ;
+ endfor
+ fprintf (file, " 0\n") ;
+ endif
+
+ fclose (file) ;
+ endif
+
+endfunction
+
+# Do not edit or modify anything in this comment block.
+# The arch-tag line is a file identity tag for the GNU Arch
+# revision control system.
+#
+# arch-tag: 2f1ff4fa-ea6a-4e54-a5f8-dad55def9834
+
--- /dev/null
+++ b/Octave/measure_filter.m
@@ -1,0 +1,70 @@
+function [stop_atten stop_band_start minus_3db] = measure_filter (f, atten)
+
+
+spec_len = 400000 ;
+
+# Calculate the spectrum.
+
+spec = 20 * log10 (abs (fft ([f zeros(1, spec_len - length (f))]))) ;
+
+spec = spec (1:spec_len/2) ;
+
+#-------------------------------------------------------------------------------
+# Find the first null which starts off the stop band.
+
+first_null = 0 ;
+for k = 2:length (spec) - 1,
+ if spec (k) < -0.8 * atten && spec (k-1) > spec (k) && spec (k) < spec (k + 1),
+ first_null = k ;
+ break
+ endif
+ endfor
+
+#-------------------------------------------------------------------------------
+# Find the stop band minimum attenuation.
+
+stop_atten = max (spec (first_null:length (spec))) ;
+
+#-------------------------------------------------------------------------------
+# Find the x position on the transition band which has the same attenuation.
+
+atten_start = 0 ;
+for k = 1:first_null,
+ if spec (k) > stop_atten && spec (k + 1) < stop_atten,
+ atten_start = k ;
+ break ;
+ endif
+ endfor
+
+atten_start = atten_start - 1 ; # Arrays are 1 based so subtract 1.
+
+stop_band_start = atten_start + (stop_atten - spec (atten_start)) / (spec (atten_start+1) - spec (atten_start)) ;
+
+
+stop_band_start = stop_band_start / spec_len ;
+
+#-------------------------------------------------------------------------------
+# Find -3db point.
+
+minus_3db = 0 ;
+for k = 1:first_null,
+ if spec (k) > -3.0 && spec (k + 1) < -3.0,
+ minus_3db = k ;
+ break ;
+ endif
+ endfor
+
+minus_3db = minus_3db - 1 ; # Arrays are 1 based so subtract 1.
+
+minus_3db = minus_3db + (stop_atten - spec (minus_3db)) / (spec (minus_3db+1) - spec (minus_3db)) ;
+
+minus_3db = minus_3db / spec_len ;
+
+endfunction
+
+# Do not edit or modify anything in this comment block.
+# The arch-tag line is a file identity tag for the GNU Arch
+# revision control system.
+#
+# arch-tag: cc2bc9a2-d387-4fed-aa0a-570e91f17c99
+