Quantcast

About the derivation in firls.m

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

About the derivation in firls.m

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

Re: About the derivation in firls.m

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

Re: About the derivation in firls.m

Quentin Spencer-4
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.


On Apr 13, 2017, at 10:19 AM, Nicholas Jankowski <[hidden email]> wrote:

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

Re: About the derivation in firls.m

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

Re: About the derivation in firls.m

Quentin Spencer-4
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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: About the derivation in firls.m

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

Re: About the derivation in firls.m

Sergei Steshenko
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
Loading...