History

Move octave scripts to Octave dir.

Erik de Castro Lopo committed 19 years ago

Commit 497f1c61 parent ad6e0fd3 files · +0 −231Browse files

Unified · Split

3 files
  1. generate_filter.m +0 −35
  2. make_src_filter.m +0 −126
  3. measure_filter.m +0 −70
generate_filter.m+0−35
@@ -1,35 +0,0 @@
1function f = generate_filter (cycles, fudge_factor, increment, atten)
2
3if nargin != 4,
4 error ("Need four args.") ;
5 endif
6
7# Calclate N and make sure it is odd.
8N = fix (4 * cycles * fudge_factor * increment) ;
9
10if rem (N, 2) != 1,
11 N = N - 1 ;
12 endif
13
14# Generate the Sinc function.
15
16m = -((N-1)/2):((N-1)/2) ;
17f = sinc (m / fudge_factor / increment) ;
18
19# Genertate the window function and apply it.
20
21w = kaiser (N, (atten + 0.5) / 10) ;
22w = w' ;
23
24f = f .* w ;
25
26f = f / sum (f) ;
27
28endfunction
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
make_src_filter.m+0−126
@@ -1,126 +0,0 @@
1function f = make_src_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
12if nargin < 3,
13 error ('Try make_filter (12, 32, 88, "output.txt")') ;
14 endif
15
16if nargin < 4,
17 filename = 0 ;
18elseif (isstr (filename) == 0),
19 error ("Fourth parameter must be a file name.") ;
20 endif
21
22fudge_factor1 = 1.0 ;
23f1 = generate_filter (cycles, fudge_factor1, increment, atten) ;
24[stop_atten stop_band_start1 minus_3db] = measure_filter (f1, atten) ;
25printf (" fudge_factor : %15.13f stop_band_start : %15.13f 1\n", fudge_factor1, stop_band_start1) ;
26
27fudge_factor2 = 1.25 ;
28f2 = generate_filter (cycles, fudge_factor2, increment, atten) ;
29[stop_atten stop_band_start2 minus_3db] = measure_filter (f2, atten) ;
30printf (" fudge_factor : %15.13f stop_band_start : %15.13f 2\n", fudge_factor2, stop_band_start2) ;
31
32f = f1 ;
33fudge_factor = fudge_factor1 ;
34stop_band_start = stop_band_start1 ;
35
36while ((stop_band_start1 - stop_band_start2) > 0.00000001)
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
66printf ("\n") ;
67
68#-------------------------------------------------------------------------------
69# Grab only half the coefficients.
70
71N = length (f) ;
72
73f = increment * f' ;
74
75if rem (length (f), 2) == 1,
76 index = find (f == max (f)) ;
77 half_f = f (index:length (f)) ;
78else
79 error ("Length should be odd.") ;
80 half_f = f ((length(f)/2):length(f)) ;
81 endif
82
83#-------------------------------------------------------------------------------
84# Print analysis.
85
86printf ("# f = make_filter (%d, %d, %4.1f) ;\n", cycles, increment, atten) ;
87printf ("# Coeff. count : %d\n", N) ;
88printf ("# Fudge factor : %9.7f\n", fudge_factor) ;
89printf ("# Pass band width : %12.10f (should be %12.10f)\n", stop_band_start, 0.5 / increment) ;
90printf ("# Stop band atten. : %5.2f dB\n", abs (stop_atten)) ;
91printf ("# -3dB band Width : %5.3f\n", 0.5 / increment / minus_3db) ;
92printf ("# half length : %d\n", length (half_f)) ;
93printf ("# increment : %d\n", increment) ;
94
95if filename,
96 file = fopen (filename, "w") ;
97 if file == 0,
98 str = sprintf ("Error, not able to open '%s'", filename)
99 error (str) ;
100 endif
101
102 fprintf (file, "/*\n") ;
103 fprintf (file, "** f = make_filter (%d, %d, %4.1f) ;\n", cycles, increment, atten) ;
104 fprintf (file, "** Pass band width : %9.7f (should be %9.7f)\n", stop_band_start, 0.5 / increment) ;
105 fprintf (file, "** Stop band atten. : %5.2f dB\n", abs (stop_atten)) ;
106 fprintf (file, "** -3dB band width : %5.3f\n", 0.5 / increment / minus_3db) ;
107 fprintf (file, "** half length : %d\n", length (half_f)) ;
108 fprintf (file, "** increment : %d\n", increment) ;
109 fprintf (file, "*/\n\n") ;
110
111 for val = half_f,
112 fprintf (file, "% 24.20e,\n", val) ;
113 endfor
114
115 fprintf (file, " 0\t\t\t/* Need a final zero coefficient */\n\n") ;
116 fclose (file) ;
117 endif
118
119endfunction
120
121# Do not edit or modify anything in this comment block.
122# The arch-tag line is a file identity tag for the GNU Arch
123# revision control system.
124#
125# arch-tag: 2f1ff4fa-ea6a-4e54-a5f8-dad55def9834
126
measure_filter.m+0−70
@@ -1,70 +0,0 @@
1function [stop_atten stop_band_start minus_3db] = measure_filter (f, atten)
2
3
4spec_len = 400000 ;
5
6# Calculate the spectrum.
7
8spec = 20 * log10 (abs (fft ([f zeros(1, spec_len - length (f))]))) ;
9
10spec = spec (1:spec_len/2) ;
11
12#-------------------------------------------------------------------------------
13# Find the first null which starts off the stop band.
14
15first_null = 0 ;
16for 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
26stop_atten = max (spec (first_null:length (spec))) ;
27
28#-------------------------------------------------------------------------------
29# Find the x position on the transition band which has the same attenuation.
30
31atten_start = 0 ;
32for 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
39atten_start = atten_start - 1 ; # Arrays are 1 based so subtract 1.
40
41stop_band_start = atten_start + (stop_atten - spec (atten_start)) / (spec (atten_start+1) - spec (atten_start)) ;
42
43
44stop_band_start = stop_band_start / spec_len ;
45
46#-------------------------------------------------------------------------------
47# Find -3db point.
48
49minus_3db = 0 ;
50for 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
57minus_3db = minus_3db - 1 ; # Arrays are 1 based so subtract 1.
58
59minus_3db = minus_3db + (stop_atten - spec (minus_3db)) / (spec (minus_3db+1) - spec (minus_3db)) ;
60
61minus_3db = minus_3db / spec_len ;
62
63endfunction
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

Keyboard shortcuts

General

Show keyboard shortcuts
?
Focus the current search
/

Issues and pull requests

Open filters
f

Issues, pull requests and Inbox

Select next item
j↓
Select previous item
k↑
Open selected item
Enter
Dismiss filters or clear selection
Esc

Repository code

Go to file
t