# About the derivation in firls.m

9 messages
Open this post in threaded view
|
Report Content as Inappropriate

## About the derivation in firls.m

 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
Open this post in threaded view
|
Report Content as Inappropriate

## Re: About the derivation in firls.m

 On Apr 13, 2017 11:11 AM, "je suis" <[hidden email]> wrote: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 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
Open this post in threaded view
|
Report Content as Inappropriate

## Re: About the derivation in firls.m

Open this post in threaded view
|
Report Content as Inappropriate

## Re: About the derivation in firls.m

 > 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
Open this post in threaded view
|
Report Content as Inappropriate

## 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 _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|
Report Content as Inappropriate

## Re: About the derivation in firls.m

 On Mon, Apr 24, 2017 at 11:28 AM, Quentin Spencer wrote:I don’t have the source code readily available. Can you show me where it is now? Signal repository:https://sourceforge.net/p/octave/signal/ci/default/tree/I believe it's under the inst/ folder _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|
Report Content as Inappropriate

## Re: About the derivation in firls.m

Open this post in threaded view
|
Report Content as Inappropriate

## Re: About the derivation in firls.m

 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