Add Octave scripts to generate filter coefficients

Erik de Castro Lopo committed 7 years ago

Commit 29b0b171 parent d1054fb5 files · +269 −1Browse files

Unified · Split

5 files

Makefile.am +3 −1
Octave/Readme.md +27 −0
Octave/generate_filter.m +35 −0
Octave/make_filter.m +134 −0
Octave/measure_filter.m +70 −0

@@ -4,7 +4,9 @@ DISTCHECK_CONFIGURE_FLAGS = --enable-werror · 3 unmodified lines
4		4
5	AM_CPPFLAGS = -I$(top_srcdir)/src $(OS_SPECIFIC_INCLUDES)	5	AM_CPPFLAGS = -I$(top_srcdir)/src $(OS_SPECIFIC_INCLUDES)
6		6
7	EXTRA_DIST = autogen.sh libsamplerate.spec.in samplerate.pc.in Make.bat Win32	7	EXTRA_DIST = autogen.sh libsamplerate.spec.in samplerate.pc.in Make.bat Win32 \
		8	Octave/generate_filter.m Octave/make_filter.m Octave/measure_filter.m \
		9	Octave/Readme.md
8		10
9	pkgconfig_DATA = samplerate.pc	11	pkgconfig_DATA = samplerate.pc
10		12

Octave/Readme.md+27−0

@@ -0,0 +1,27 @@
	1	# Gneerating Filter Coefficients
	2
	3	The scripts below are used in [GNU Octave][GNU Octave] which should be available
	4	on most platforms including Linux, Mac and Windows.
	5
	6	On Debian and other Debian derived Linux distributions, Octave can be installed
	7	using: `sudo apt install octave octave-signal`.
	8
	9	Filter coefficients can be generated as follows by running the `octave` command
	10	in the same directory as the scripts and this Readme file:
	11
	12	```
	13	octave:1> pkg load signal
	14	octave:2> f = make_filter (8, 128, 100.3) ;
	15	# f = make_fip_filter (8, 128, 100.3) ;
	16	# Coeff. count : 4922
	17	# Fudge factor : 1.2018769
	18	# Pass band width : 0.0039062499 (should be 0.0039062500)
	19	# Stop band atten. : 100.71 dB
	20	# -3dB band Width : 0.490
	21	# half length : 2464
	22	# increment : 128
	23	```
	24
	25	The parameters above generate the coefficient set used in `src/fastest_coeffs.h`.
	26
	27	[GNU Octave]: https://www.gnu.org/software/octave/

Octave/generate_filter.m+35−0

@@ -0,0 +1,35 @@
	1	function f = generate_filter (cycles, fudge_factor, increment, atten)
	2
	3	if nargin != 4,
	4	error ("Need four args.") ;
	5	endif
	6
	7	# Calclate N and make sure it is even.
	8	N = fix (4 * cycles * fudge_factor * increment) ;
	9
	10	if rem (N, 2) != 0,
	11	N = N - 1 ;
	12	endif
	13
	14	# Generate the Sinc function.
	15
	16	m = -((N-1)/2):((N-1)/2) ;
	17	f = sinc (m / fudge_factor / increment) ;
	18
	19	# Genertate the window function and apply it.
	20
	21	w = kaiser (N, (atten + 0.5) / 10) ;
	22	w = w' ;
	23
	24	f = f .* w ;
	25
	26	f = f / sum (f) ;
	27
	28	endfunction
	29
	30	# Do not edit or modify anything in this comment block.
	31	# The arch-tag line is a file identity tag for the GNU Arch
	32	# revision control system.
	33	#
	34	# arch-tag: 7e57a3cb-3f5c-4346-bfcd-4da1e758e2a7
	35

Octave/make_filter.m+134−0

@@ -0,0 +1,134 @@
	1	function f = make_filter (cycles, increment, atten, filename)
	2
	3	# This works :
	4	#
	5	# f = make_filter (67, 128, 100.3) ;
	6	# f = make_filter (13, 128, 100.5) ;
	7	# f = make_filter (185, 4, 157.0) ;
	8
	9
	10	#=======================================================================
	11
	12	if nargin < 3,
	13	error ('Try make_filter (12, 32, 88, "output.txt")') ;
	14	endif
	15
	16	if nargin < 4,
	17	filename = 0 ;
	18	elseif (isstr (filename) == 0),
	19	error ("Fourth parameter must be a file name.") ;
	20	endif
	21
	22	fudge_factor1 = 1.0 ;
	23	f1 = generate_filter (cycles, fudge_factor1, increment, atten) ;
	24	[stop_atten stop_band_start1 minus_3db] = measure_filter (f1, atten) ;
	25	printf (" fudge_factor : %15.13f stop_band_start : %15.13f 1\n", fudge_factor1, stop_band_start1) ;
	26
	27	fudge_factor2 = 1.25 ;
	28	f2 = generate_filter (cycles, fudge_factor2, increment, atten) ;
	29	[stop_atten stop_band_start2 minus_3db] = measure_filter (f2, atten) ;
	30	printf (" fudge_factor : %15.13f stop_band_start : %15.13f 2\n", fudge_factor2, stop_band_start2) ;
	31
	32	f = f1 ;
	33	fudge_factor = fudge_factor1 ;
	34	stop_band_start = stop_band_start1 ;
	35
	36	while ((stop_band_start1 - stop_band_start2) > 0.0000000001)
	37	if (stop_band_start1 < stop_band_start2)
	38	printf ("stop_band_start1 < stop_band_start2\n") ;
	39	break ;
	40	endif
	41
	42	fudge_factor = fudge_factor1 + (fudge_factor2 - fudge_factor1) / 2 ;
	43	f = generate_filter (cycles, fudge_factor, increment, atten) ;
	44	[stop_atten stop_band_start minus_3db] = measure_filter (f, atten) ;
	45
	46	if (stop_band_start > 1.0)
	47	printf ("A %10.8f %10.8f %10.8f\n", fudge_factor1, fudge_factor, fudge_factor2) ;
	48	continue ;
	49	endif
	50
	51	if (stop_band_start < 0.5 / increment)
	52	f2 = f ;
	53	stop_band_start2 = stop_band_start ;
	54	fudge_factor2 = fudge_factor ;
	55	choice = 2 ;
	56	else
	57	f1 = f ;
	58	stop_band_start1 = stop_band_start ;
	59	fudge_factor1 = fudge_factor ;
	60	choice = 1 ;
	61	endif
	62
	63	printf (" fudge_factor : %15.13f stop_band_start : %15.13f %d\n", fudge_factor, stop_band_start, choice) ;
	64	endwhile
	65
	66	printf ("\n") ;
	67
	68	#-------------------------------------------------------------------------------
	69	# Grab only half the coefficients.
	70
	71	N = length (f) ;
	72
	73	f = increment * f' ;
	74
	75	if rem (length (f), 2) == 0,
	76	index = find (f == max (f)) ;
	77	index = min (index) - 1 ;
	78	half_f = f (index:length (f)) ;
	79	else
	80	error ("Length should be even.") ;
	81	endif
	82
	83	trailing_zeros = 4 - rem (length (half_f), 4) ;
	84
	85	#-------------------------------------------------------------------------------
	86	# Print analysis.
	87
	88	printf ("# f = make_filter (%d, %d, %4.1f) ;\n", cycles, increment, atten) ;
	89	printf ("# Coeff. count : %d\n", N) ;
	90	printf ("# Fudge factor : %9.7f\n", fudge_factor) ;
	91	printf ("# Pass band width : %12.10f (should be %12.10f)\n", stop_band_start, 0.5 / increment) ;
	92	printf ("# Stop band atten. : %5.2f dB\n", abs (stop_atten)) ;
	93	printf ("# -3dB band Width : %5.3f\n", 0.5 / increment / minus_3db) ;
	94	printf ("# half length : %d\n", length (half_f) + trailing_zeros) ;
	95	printf ("# increment : %d\n", increment) ;
	96
	97	if filename,
	98	file = fopen (filename, "w") ;
	99	if file == 0,
	100	str = sprintf ("Error, not able to open '%s'", filename)
	101	error (str) ;
	102	endif
	103
	104	fprintf (file, "/*\n") ;
	105	fprintf (file, "** f = make_filter (%d, %d, %4.1f) ;\n", cycles, increment, atten) ;
	106	fprintf (file, "** Pass band width : %9.7f (should be %9.7f)\n", stop_band_start, 0.5 / increment) ;
	107	fprintf (file, "** Stop band atten. : %5.2f dB\n", abs (stop_atten)) ;
	108	fprintf (file, "** -3dB band width : %5.3f\n", 0.5 / increment / minus_3db) ;
	109	fprintf (file, "** half length : %d\n", length (half_f)) ;
	110	fprintf (file, "** increment : %d\n", increment) ;
	111	fprintf (file, "*/\n\n") ;
	112
	113	for val = half_f,
	114	fprintf (file, "% 24.20e,\n", val) ;
	115	endfor
	116
	117	if trailing_zeros > 0,
	118	for val = 2:trailing_zeros,
	119	fprintf (file, " 0,\n") ;
	120	endfor
	121	fprintf (file, " 0\n") ;
	122	endif
	123
	124	fclose (file) ;
	125	endif
	126
	127	endfunction
	128
	129	# Do not edit or modify anything in this comment block.
	130	# The arch-tag line is a file identity tag for the GNU Arch
	131	# revision control system.
	132	#
	133	# arch-tag: 2f1ff4fa-ea6a-4e54-a5f8-dad55def9834
	134

Octave/measure_filter.m+70−0

@@ -0,0 +1,70 @@
	1	function [stop_atten stop_band_start minus_3db] = measure_filter (f, atten)
	2
	3
	4	spec_len = 400000 ;
	5
	6	# Calculate the spectrum.
	7
	8	spec = 20 * log10 (abs (fft ([f zeros(1, spec_len - length (f))]))) ;
	9
	10	spec = spec (1:spec_len/2) ;
	11
	12	#-------------------------------------------------------------------------------
	13	# Find the first null which starts off the stop band.
	14
	15	first_null = 0 ;
	16	for k = 2:length (spec) - 1,
	17	if spec (k) < -0.8 * atten && spec (k-1) > spec (k) && spec (k) < spec (k + 1),
	18	first_null = k ;
	19	break
	20	endif
	21	endfor
	22
	23	#-------------------------------------------------------------------------------
	24	# Find the stop band minimum attenuation.
	25
	26	stop_atten = max (spec (first_null:length (spec))) ;
	27
	28	#-------------------------------------------------------------------------------
	29	# Find the x position on the transition band which has the same attenuation.
	30
	31	atten_start = 0 ;
	32	for k = 1:first_null,
	33	if spec (k) > stop_atten && spec (k + 1) < stop_atten,
	34	atten_start = k ;
	35	break ;
	36	endif
	37	endfor
	38
	39	atten_start = atten_start - 1 ; # Arrays are 1 based so subtract 1.
	40
	41	stop_band_start = atten_start + (stop_atten - spec (atten_start)) / (spec (atten_start+1) - spec (atten_start)) ;
	42
	43
	44	stop_band_start = stop_band_start / spec_len ;
	45
	46	#-------------------------------------------------------------------------------
	47	# Find -3db point.
	48
	49	minus_3db = 0 ;
	50	for k = 1:first_null,
	51	if spec (k) > -3.0 && spec (k + 1) < -3.0,
	52	minus_3db = k ;
	53	break ;
	54	endif
	55	endfor
	56
	57	minus_3db = minus_3db - 1 ; # Arrays are 1 based so subtract 1.
	58
	59	minus_3db = minus_3db + (stop_atten - spec (minus_3db)) / (spec (minus_3db+1) - spec (minus_3db)) ;
	60
	61	minus_3db = minus_3db / spec_len ;
	62
	63	endfunction
	64
	65	# Do not edit or modify anything in this comment block.
	66	# The arch-tag line is a file identity tag for the GNU Arch
	67	# revision control system.
	68	#
	69	# arch-tag: cc2bc9a2-d387-4fed-aa0a-570e91f17c99
	70

Add Octave scripts to generate filter coefficients

Keyboard shortcuts

General

Issues and pull requests

Issues, pull requests and Inbox

Repository code