-
Notifications
You must be signed in to change notification settings - Fork 1
/
firls_my.m
292 lines (271 loc) · 9.62 KB
/
firls_my.m
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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
function [h,a]=firls_my(N,F,M,W,ftype)
% FIRLS Linear-phase FIR filter design using least-squares error minimization.
% B=FIRLS(N,F,A) returns a length N+1 linear phase (real, symmetric
% coefficients) FIR filter which has the best approximation to the
% desired frequency response described by F and A in the least squares
% sense. F is a vector of frequency band edges in pairs, in ascending
% order between 0 and 1. 1 corresponds to the Nyquist frequency or half
% the sampling frequency. A is a real vector the same size as F
% which specifies the desired amplitude of the frequency response of the
% resultant filter B. The desired response is the line connecting the
% points (F(k),A(k)) and (F(k+1),A(k+1)) for odd k; FIRLS treats the
% bands between F(k+1) and F(k+2) for odd k as "transition bands" or
% "don't care" regions. Thus the desired amplitude is piecewise linear
% with transition bands. The integrated squared error is minimized.
%
% For filters with a gain other than zero at Fs/2, e.g., highpass
% and bandstop filters, N must be even. Otherwise, N will be
% incremented by one. Alternatively, you can use a trailing 'h' flag to
% design a type 4 linear phase filter and avoid incrementing N.
%
% B=FIRLS(N,F,A,W) uses the weights in W to weight the error. W has one
% entry per band (so it is half the length of F and A) which tells
% FIRLS how much emphasis to put on minimizing the integral squared error
% in each band relative to the other bands.
%
% B=FIRLS(N,F,A,'Hilbert') and B=FIRLS(N,F,A,W,'Hilbert') design filters
% that have odd symmetry, that is, B(k) = -B(N+2-k) for k = 1, ..., N+1.
% A special case is a Hilbert transformer which has an approx. amplitude
% of 1 across the entire band, e.g. B=FIRLS(30,[.1 .9],[1 1],'Hilbert').
%
% B=FIRLS(N,F,A,'differentiator') and B=FIRLS(N,F,A,W,'differentiator')
% also design filters with odd symmetry, but with a special weighting
% scheme for non-zero amplitude bands. The weight is assumed to be equal
% to the inverse of frequency, squared, times the weight W. Thus the
% filter has a much better fit at low frequency than at high frequency.
% This designs FIR differentiators.
%
% % Example of a length 31 lowpass filter.
% h=firls(30,[0 .1 .2 .5]*2,[1 1 0 0]);
% fvtool(h);
%
% % Example of a length 45 lowpass differentiator.
% h=firls(44,[0 .3 .4 1],[0 .2 0 0],'differentiator');
% fvtool(h);
%
% % Example of a length 26 type 4 highpass filter.
% h=firls(25,[0 .4 .5 1],[0 0 1 1],'h');
% fvtool(h);
%
% See also FIRPM, FIR1, FIR2, FREQZ and FILTER.
% Author(s): T. Krauss
% History: 10-18-91, original version
% 3-30-93, updated
% 9-1-95, optimize adjacent band case
% Copyright 1988-2011 The MathWorks, Inc.
%
% check number of arguments, set up defaults.
error(nargchk(3,5,nargin,'struct'));
if (max(F)>1) || (min(F)<0)
error(message('signal:firls:InvalidRange'))
end
if (rem(length(F),2)~=0)
error(message('signal:firls:MustHaveEvenLength', 'F'));
end
if (length(F) ~= length(M))
error(message('signal:firls:UnequalLengths', 'F', 'A'));
end
if (nargin==3),
W = ones(length(F)/2,1);
ftype = '';
end
if (nargin==4),
if ischar(W),
ftype = W; W = ones(length(F)/2,1);
else
ftype = '';
end
end
if (nargin==5),
if isempty(W),
W = ones(length(F)/2,1);
end
end
if isempty(ftype)
ftype = 0; differ = 0;
else
ftype = lower(ftype);
if strcmpi(ftype,'h') || strcmpi(ftype,'hilbert')
ftype = 1; differ = 0;
elseif strcmpi(ftype,'d') || strcmpi(ftype,'differentiator')
ftype = 1; differ = 1;
else
error(message('signal:firls:InvalidEnum'))
end
end
% Check for valid filter length
%[N,msg1,msg2,msgobj] = firchk(N,F(end),M,ftype);
%if ~isempty(msg1), error(msgobj); end;
%if ~isempty(msg2), warning(msgobj); end;
N = N+1; % filter length
F=F(:)/2; M=M(:); W=sqrt(W(:)); % make these guys columns
dF = diff(F);
if (length(F) ~= length(W)*2)
error(message('signal:firls:InvalidDimensions'));
end
if any(dF<0),
error(message('signal:firls:InvalidFreqVec'))
end
% Fix for 67187
if all(dF(2:2:length(dF)-1)==0) && length(dF) > 1,
fullband = 1;
else
fullband = 0;
end
if all((W-W(1))==0)
constant_weights = 1;
else
constant_weights = 0;
end
L=(N-1)/2;
Nodd = rem(N,2);
if (ftype == 0), % Type I and Type II linear phase FIR
% basis vectors are cos(2*pi*m*f) (see m below)
if ~Nodd
m=(0:L)+.5; % type II
else
m=(0:L); % type I
end
k=m';
need_matrix = (~fullband) || (~constant_weights);
if need_matrix
I1=k(:,ones(size(m)))+m(ones(size(k)),:); % entries are m + k
I2=k(:,ones(size(m)))-m(ones(size(k)),:); % entries are m - k
G=zeros(size(I1));
end
if Nodd
k=k(2:length(k));
b0=0; % first entry must be handled separately (where k(1)=0)
end;
b=zeros(size(k));
for s=1:2:length(F),
m=(M(s+1)-M(s))/(F(s+1)-F(s)); % slope
b1=M(s)-m*F(s); % y-intercept
if Nodd
b0 = b0 + (b1*(F(s+1)-F(s)) + m/2*(F(s+1)*F(s+1)-F(s)*F(s)))...
* abs(W((s+1)/2)^2) ;
end
b = b+(m/(4*pi*pi)*(cos(2*pi*k*F(s+1))-cos(2*pi*k*F(s)))./(k.*k))...
* abs(W((s+1)/2)^2);
b = b + (F(s+1)*(m*F(s+1)+b1)*sinc(2*k*F(s+1)) ...
- F(s)*(m*F(s)+b1)*sinc(2*k*F(s))) ...
* abs(W((s+1)/2)^2);
if need_matrix
G = G + (.5*F(s+1)*(sinc(2*I1*F(s+1))+sinc(2*I2*F(s+1))) ...
- .5*F(s)*(sinc(2*I1*F(s))+sinc(2*I2*F(s))) ) ...
* abs(W((s+1)/2)^2);
end
end;
if Nodd
b=[b0; b];
end;
if need_matrix
a=G\b;
else
a=(W(1)^2)*4*b;
if Nodd
a(1) = a(1)/2;
end
end
if Nodd
h=[a(L+1:-1:2)/2; a(1); a(2:L+1)/2].';
else
h=.5*[flipud(a); a].';
end;
elseif (ftype == 1), % Type III and Type IV linear phase FIR
% basis vectors are sin(2*pi*m*f) (see m below)
if (differ), % weight non-zero bands with 1/f^2
do_weight = ( abs(M(1:2:length(M))) + abs(M(2:2:length(M))) ) > 0;
else
do_weight = zeros(size(F));
end
if Nodd
m=(1:L); % type III
else
m=(0:L)+.5; % type IV
end;
k=m';
b=zeros(size(k));
need_matrix = (~fullband) || (any(do_weight)) || (~constant_weights);
if need_matrix
I1=k(:,ones(size(m)))+m(ones(size(k)),:); % entries are m + k
I2=k(:,ones(size(m)))-m(ones(size(k)),:); % entries are m - k
G=zeros(size(I1));
end
for s=1:2:length(F),
if (do_weight((s+1)/2)), % weight bands with 1/f^2
if F(s) == 0, F(s) = 1e-5; end % avoid singularities
m=(M(s+1)-M(s))/(F(s+1)-F(s));
b1=M(s)-m*F(s);
snint1 = sineint(2*pi*k*F(s+1)) - sineint(2*pi*k*F(s));
%snint1 = (-1/2/i)*(expint(1i*2*pi*k*F(s+1)) ...
% -expint(-i*2*pi*k*F(s+1)) -expint(1i*2*pi*k*F(s)) ...
% +expint(-i*2*pi*k*F(s)) );
% csint1 = cosint(2*pi*k*F(s+1)) - cosint(2*pi*k*F(s)) ;
csint1 = (-1/2)*(expint(1i*2*pi*k*F(s+1))+expint(-1i*2*pi*k*F(s+1))...
-expint(1i*2*pi*k*F(s)) -expint(-1i*2*pi*k*F(s)) );
b=b + ( m*snint1 ...
+ b1*2*pi*k.*( -sinc(2*k*F(s+1)) + sinc(2*k*F(s)) + csint1 ))...
* abs(W((s+1)/2)^2);
snint1 = sineint(2*pi*F(s+1)*(-I2));
snint2 = sineint(2*pi*F(s+1)*I1);
snint3 = sineint(2*pi*F(s)*(-I2));
snint4 = sineint(2*pi*F(s)*I1);
G = G - ( ( -1/2*( cos(2*pi*F(s+1)*(-I2))/F(s+1) ...
- 2*snint1*pi.*I2 ...
- cos(2*pi*F(s+1)*I1)/F(s+1) ...
- 2*snint2*pi.*I1 )) ...
- ( -1/2*( cos(2*pi*F(s)*(-I2))/F(s) ...
- 2*snint3*pi.*I2 ...
- cos(2*pi*F(s)*I1)/F(s) ...
- 2*snint4*pi.*I1) ) ) ...
* abs(W((s+1)/2)^2);
else % use usual weights
m=(M(s+1)-M(s))/(F(s+1)-F(s));
b1=M(s)-m*F(s);
b=b+(m/(4*pi*pi)*(sin(2*pi*k*F(s+1))-sin(2*pi*k*F(s)))./(k.*k))...
* abs(W((s+1)/2)^2) ;
b = b + (((m*F(s)+b1)*cos(2*pi*k*F(s)) - ...
(m*F(s+1)+b1)*cos(2*pi*k*F(s+1)))./(2*pi*k)) ...
* abs(W((s+1)/2)^2) ;
if need_matrix
G = G + (.5*F(s+1)*(sinc(2*I1*F(s+1))-sinc(2*I2*F(s+1))) ...
- .5*F(s)*(sinc(2*I1*F(s))-sinc(2*I2*F(s)))) * ...
abs(W((s+1)/2)^2);
end
end;
end
if need_matrix
a=G\b;
else
a=-4*b*(W(1)^2);
end
if Nodd
h=.5*[flipud(a); 0; -a].';
else
h=.5*[flipud(a); -a].';
end
if differ, h=-h; end
end
if nargout > 1
a = 1;
end
%----------------------------------------------------------------------------
function y = sineint(x)
% SINEINT (a.k.a. SININT) Numerical Sine Integral
% Used by FIRLS in the Signal Processing Toolbox.
% Untested for complex or imaginary inputs.
%
% See also SININT in the Symbolic Toolbox.
% Was Revision: 1.5, Date: 1996/03/15 20:55:51
i1 = find(real(x)<0); % this equation is not valid if x is in the
% left-hand plane of the complex plane.
% use relation Si(-z) = -Si(z) in this case (Eq 5.2.19, Abramowitz
% & Stegun).
x(i1) = -x(i1);
y = zeros(size(x));
ind = find(x);
% equation 5.2.21 Abramowitz & Stegun
% y(ind) = (1/(2*i))*(expint(1i*x(ind)) - expint(-i*x(ind))) + pi/2;
y(ind) = imag(expint(1i*x(ind))) + pi/2;
y(i1) = -y(i1);