Quantcast

About the derivation in firls.m

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
9 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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: About the derivation in firls.m

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

Re: About the derivation in firls.m

je suis
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
Loading...