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 |
Free forum by Nabble | Edit this page |