1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
|
function b = fir1(n, w, varargin)
funcprot(0);
if argn(2) < 2 | argn(2) > 5
error("Wrong Number of input arguments");
end
// Assign default window, filter type and scale.
// If single band edge, the first band defaults to a pass band to
// create a lowpass filter. If multiple band edges, the first band
// defaults to a stop band so that the two band case defaults to a
// band pass filter. Ick.
window_in = [];
scale = 1;
ftype = bool2s(length(w)==1);
for i=1:length(varargin)
arg = varargin(i);
if (type(arg)==10)
arg=convstr(arg,"l");
end
if isempty(arg)
continue;
end
select arg
case 'low' then ftype = 1; case 'stop' then ftype = 1; case 'dc-1' then ftype = 1;
case 'high' then ftype = 0; case 'pass' then ftype = 0; case 'bandpass' then ftype = 0; case 'dc-0' then ftype = 0;
case 'scale' then scale=1;
case 'noscale' then scale=0;
else window_in=arg;
end
end
// build response function according to fir2 requirements
bands = length(w)+1;
f = zeros(1,2*bands);
f(2*bands)=1;
f(2:2:2*bands-1) = w;
f(3:2:2*bands-1) = w;
m = zeros(1,2*bands);
m(1:2:2*bands) = modulo([1:bands]-(1-ftype),2);
m(2:2:2*bands) = m(1:2:2*bands);
//Increment the order if the final band is a pass band. Something
// about having a nyquist frequency of zero causing problems.
//
if modulo(n,2)==1 & m(2*bands)==1,
warning("n must be even for highpass and bandstop filters. Incrementing.");
n = n+1;
if isvector(window_in) & isreal(window_in) & ~(type(window_in)==10)
// Extend the window using interpolation
M = length(window_in);
if M == 1,
window_in = [window_in; window_in];
elseif M < 4
window_in = interp1(linspace(0,1,M),window_in,linspace(0,1,M+1),'linear');
else
window_in = interp1(linspace(0,1,M),window_in,linspace(0,1,M+1),'spline');
end
end
end
// compute the filter
b = fir2(n, f, m, [], 2, window_in);
// normalize filter magnitude
if scale == 1
// find the middle of the first band edge
// find the frequency of the normalizing gain
if m(1) == 1
// if the first band is a passband, use DC gain
w_o = 0;
elseif f(4) == 1
// for a highpass filter,
// use the gain at half the sample frequency
w_o = 1;
else
// otherwise, use the gain at the center
// frequency of the first passband
w_o = f(3) + (f(4)-f(3))/2;
end
// compute |h(w_o)|^-1
if ~(isvector(b) | isempty(b)) // Check input is a vector
error('Invalid');
end
x=exp(-1*%i*%pi*w_o)
// z=[1 -exp(-1*%i*%pi*w_o)];
disp(x)
nc = length(b);
if(isscalar(x) & nc>0 & (x~=%inf) & or(b(:)~=%inf))
// Make it scream for scalar x. Polynomial evaluation can be
// implemented as a recursive digital filter.
q=b;
k = filter(1,[1 -real(x)],q);
k=k(nc);
end
k=abs(k);
renorm = 1/k
// normalize the filter
b = renorm*b;
end
endfunction
|
