Hello everyone
And I hope this isn't the wrong list to send this to. I am trying to make a small program in C++ about filter design, and I am currently trying the least-squares algorithm, based on Selesnick's paper (the same one from firls.m). However, I would like to be able to also make type iii and iv filters, differentiators and Hilbert transformers, too, and I was looking through your firls.m file, where I got stuck at the derivation of "cos_ints2". If it's not too much to ask, can someone please let me know how did it get from the formula in the comment to the one in the function calculating "cos_ints2"? I can't figure it out, not even with wxMaxima by my side. In the extreme case that my beginner programming skills allow it, I would also like to be able to make the weights as a linear variable function, too, so I would guess that the q vector will follow about the same principles as the b vector, so I would like to try it out, even if my small program, most probably, will never leave the screen of my monitor. But, if I get "lucky" (one is allowed to dream), I will let you know. Vlad _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
On Apr 13, 2017 11:11 AM, "je suis" <[hidden email]> wrote: Hello everyone The Signal package maintainer is copied on this email, but the copyright notice in firls.m suggests it was written in 2006 by someone else. I cc'd that person here on the chance the email is still valid and they may be able to take a moment to help you out. Nick J. _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
Yes, I’m still here. It has indeed been a very long time. I’ll try to take a look at the code and see if I am able to answer your question.
_______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
> Yes, I’m still here. It has indeed been a very long time. I’ll try to take a
> look at the code and see if I am able to answer your question. I know not everyone has the time for it, but I would like to know if this is still valid. I, for one, am patient, but if my patience interferes with someone else's life, then I'd rather not continue. Vlad _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
I don’t have the source code readily available. Can you show me where it is now?
> On Apr 23, 2017, at 9:22 AM, je suis <[hidden email]> wrote: > >> Yes, I’m still here. It has indeed been a very long time. I’ll try to take a >> look at the code and see if I am able to answer your question. > > I know not everyone has the time for it, but I would like to know if > this is still valid. I, for one, am patient, but if my patience > interferes with someone else's life, then I'd rather not continue. > > Vlad _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
On Mon, Apr 24, 2017 at 11:28 AM, Quentin Spencer <[hidden email]> wrote: I don’t have the source code readily available. Can you show me where it is now? I believe it's under the inst/ folder _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
In reply to this post by Quentin Spencer-4
________________________________ From: Quentin Spencer <[hidden email]> To: je suis <[hidden email]> Cc: help <[hidden email]>; Mike Miller <[hidden email]> Sent: Monday, April 24, 2017 6:29 PM Subject: Re: About the derivation in firls.m I don’t have the source code readily available. Can you show me where it is now? > On Apr 23, 2017, at 9:22 AM, je suis <[hidden email]> wrote: > >> Yes, I’m still here. It has indeed been a very long time. I’ll try to take a >> look at the code and see if I am able to answer your question. > > I know not everyone has the time for it, but I would like to know if > this is still valid. I, for one, am patient, but if my patience > interferes with someone else's life, then I'd rather not continue. > > Vlad _______________________________________________ The source from signal-1.2.2: ## Copyright (C) 2006 Quentin Spencer <[hidden email]> ## ## This program is free software; you can redistribute it and/or modify it under ## the terms of the GNU General Public License as published by the Free Software ## Foundation; either version 3 of the License, or (at your option) any later ## version. ## ## This program is distributed in the hope that it will be useful, but WITHOUT ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more ## details. ## ## You should have received a copy of the GNU General Public License along with ## this program; if not, see <http://www.gnu.org/licenses/>. ## b = firls(N, F, A); ## b = firls(N, F, A, W); ## ## FIR filter design using least squares method. Returns a length N+1 ## linear phase filter such that the integral of the weighted mean ## squared error in the specified bands is minimized. ## ## F specifies the frequencies of the band edges, normalized so that ## half the sample frequency is equal to 1. Each band is specified by ## two frequencies, to the vector must have an even length. ## ## A specifies the amplitude of the desired response at each band edge. ## ## W is an optional weighting function that contains one value for each ## band that weights the mean squared error in that band. A must be the ## same length as F, and W must be half the length of F. N must be ## even. If given an odd value, firls increments it by 1. ## ## The least squares optimization algorithm for computing FIR filter ## coefficients is derived in detail in: ## ## I. Selesnick, "Linear-Phase FIR Filter Design by Least Squares," ## http://cnx.org/content/m10577 function coef = firls(N, frequencies, pass, weight, str); if (nargin < 3 || nargin > 6) print_usage; elseif (nargin == 3) weight = ones (1, length(pass)/2); str = []; elseif (nargin == 4) if ischar (weight) str = weight; weight = ones (size (pass)); else str = []; endif endif if length (frequencies) ~= length (pass) error("F and A must have equal lengths."); elseif 2 * length (weight) ~= length (pass) error("W must contain one weight per band."); elseif ischar (str) error("This feature is implemented yet"); else N += mod (N, 2); M = N/2; w = kron (weight(:), [-1; 1]); omega = frequencies * pi; i1 = 1:2:length (omega); i2 = 2:2:length (omega); ## Generate the matrix Q ## As illustrated in the above-cited reference, the matrix can be ## expressed as the sum of a Hankel and Toeplitz matrix. A factor of ## 1/2 has been dropped and the final filter coefficients multiplied ## by 2 to compensate. cos_ints = [omega; sin((1:N)' * omega)]; q = [1, 1./(1:N)]' .* (cos_ints * w); Q = toeplitz (q(1:M+1)) + hankel (q(1:M+1), q(M+1:end)); ## The vector b is derived from solving the integral: ## ## _ w ## / 2 ## b = / W(w) D(w) cos(kw) dw ## k / w ## - 1 ## ## Since we assume that W(w) is constant over each band (if not, the ## computation of Q above would be considerably more complex), but ## D(w) is allowed to be a linear function, in general the function ## W(w) D(w) is linear. The computations below are derived from the ## fact that: ## _ ## / a ax + b ## / (ax + b) cos(nx) dx = --- cos (nx) + ------ sin(nx) ## / 2 n ## - n ## cos_ints2 = [omega(i1).^2 - omega(i2).^2; ... cos((1:M)' * omega(i2)) - cos((1:M)' * omega(i1))] ./ ... ([2, 1:M]' * (omega(i2) - omega(i1))); d = [-weight .* pass(i1); weight .* pass(i2)] (:); b = [1, 1./(1:M)]' .* ((kron (cos_ints2, [1, 1]) + cos_ints(1:M+1,:)) * d); ## Having computed the components Q and b of the matrix equation, ## solve for the filter coefficients. a = Q \ b; coef = [ a(end:-1:2); 2 * a(1); a(2:end) ]; endif endfunction --Sergei. _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
I've made a version of firls.m that handles even-length FIRs, too. I
am not versed in Octave, so the code will look awful, at best. The way I did it uses ifft for the b vector, but it needs fairly large sizes in order to become close to accurate. It works, though. Then, since I still don't understand how the integral got to have cos(w)/(n^2*w) from cos(w)/n^2, I just adapted the code as I saw fit. This also means I can't go further and make differentiators and Hilbert transformers, but I'll try the old hammer. It works, but it will hurt your eyes. It would be nice for someone to verify it against Matlab, I don't have it. The name of the function is arbitrary. Here it is: function h = lsfir2(N,f,A,K); % prepare a few helpers M = floor(N/2); m = [1:M+1]; oddN = mod(N, 2); bands = length(f); % calculate q... q = zeros(1, 2*M + 1 + oddN); for i=bands:-1:1 q += (-1)^i*K(floor((i - 1)/2)+1).*f(i)*sinc(f(i)*[0:length(q)-1]); end % ...and the Q matrix Q = toeplitz (q(m)) + hankel (q(m+oddN), q(m+M+oddN)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % I can't figure out cos_ints2 in firls.m so I'll just brute-force this. % Resize the frequency/amplitude/weights vectors for proper ifft, % since Selesnick derives the b[] vector from fp, rather than the usual % (fp+fs)/2. So double the frequency values except the margins, and % make K[] the same size as f[] with 0 inbetween (as the graph shows % in Selesnick's paper). Then make A*K to be used with f in ifft. % weights K2 = zeros(1, 2*length(K)-1); for i=1:2:bands-1 K2(i) = K(floor(i/2)+1); end K2 = [K2; K2](:)'; % amplitude and frequency A2 = zeros(1, 2*bands-2); f2 = A2; for i=1:2*bands-2 f2(i) = f(floor(i/2)+1); A2(i) = A(floor(i/2)+1)*K2(i); end % prepare the magnitude, then ifft, and make the b vector W = interp1(f2, A2, linspace(0,1,1024+1)); if oddN w = real(ifft([W zeros(1, 2*1024) W(end:-1:2)])); b = 2*w(2:2:2*M+2)'; else w = real(ifft([W W(end:-1:2)])); b = w(1:M+1)'; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%% % This is adapted from firls.m, to include the even length filters. % I don't understand how the integral formula from the comment became % (cos(w2)-cos(w1))/(n*(w2-w1))+sinc, so I can't go further to make % differentiators and Hilbert transformers, and general type III and IV FIRs. w = f * pi; i1 = 1:2:length (f); i2 = 2:2:length (f); if oddN cos_ints = sin((0.5:1:M+1)'*w); else cos_ints = [w; sin((1:M)' * w)]; end if ~oddN cos_ints2 = (w(i1).^2 - w(i2).^2)./(2*(w(i2) - w(i1))) else cos_ints2 = []; end cos_ints2 = [cos_ints2; cos(((1:M+oddN)-.5*oddN)' * w(i2)) - ... cos(((1:M+oddN)-.5*oddN)' * w(i1)) ./ ... (([1:M+oddN]-.5*oddN)' * (w(i2) - w(i1)))]; d = [-K .* A(i1); K .* A(i2)] (:); b = []; if ~oddN b = [1, 1./(1:M)]'; else b = [1./((1:M+1)-.5*oddN)]'; end b = b.* ((kron (cos_ints2, [1, 1]) + cos_ints(1:M+1,:)) * d)/pi; %%%%%%%%%%%%%% % the rest of the algorithm a = Q \ b; if(oddN) h = [flipud(a); a]; else h = [ a(end:-1:2); 2*a(1); a(2:end) ]; end _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
I just realized I made some changes while testing and I forgot to
include them. This is how it looks like in the cos_ints2 part: if oddN cos_ints = sin((0.5:1:M+1)'*w); cos_ints2 = [cos((0.5:1:M+1)' * w(i2)) - cos((0.5:1:M+1)' * w(i1))] ./ ... ([0.5:1:M+1]' * (w(i2) - w(i1))); else cos_ints = [w; sin((1:M)' * w)]; cos_ints2 = [w(i1).^2 - w(i2).^2; cos((1:M)' * w(i2)) - cos((1:M)' * w(i1))] ./ ... ([2,1:M]' * (w(i2) - w(i1))); end _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
Free forum by Nabble | Edit this page |