I have modified firls.m to include all types of FIRs plus HT and diff

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

I have modified firls.m to include all types of FIRs plus HT and diff

je suis
Hello everyone

I managed to modify firls.m such that it now allows all types of FIRs,
I to IV, plus differentiators and Hilbert transformers. It also allows
forbidden filters, such as type II highpass, or similar. Just because
the result won't be a highpass it doesn't mean one should not allow
it: it is assumed that the user knows what he/she/it wants and how to
get it, the option is there; if not, maybe it's better to do something
else. Also, while I could verify in parallel the results in wxMaxima
(slightly more accustomed to), I don't know if the results are,
indeed, real, so it would be nice if someone who has access to Matlab
can verify this. Matlab's site also says something about a 1/f^2
weighting, I didn't bother with that, mainly because I don't know how
to do it, but also because Selesnick's paper says nothing about it.

Since I was not accustomed to Octave's language, I tried to adapt the
code by looking at other scripts, so I don't know how much on par with
your requirements is the result. I'll attach the file. If it proves to
be usable, it would be nice to update the current firls.m, given the
new choices of design. I have built upon the previous user's code,
both to avoid reinventing the wheel and to get more accustomed to
Octave's language. The name of the file is not the same as with the
function's name to avoid name clashes, but the definition inside is
firls(...).

Vlad

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave

lsfir2.m (8K) Download Attachment
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: I have modified firls.m to include all types of FIRs plus HT and diff

NJank
On Tue, May 16, 2017 at 1:28 PM, je suis <[hidden email]> wrote:

> Hello everyone
>
> I managed to modify firls.m such that it now allows all types of FIRs,
> I to IV, plus differentiators and Hilbert transformers. It also allows
> forbidden filters, such as type II highpass, or similar. Just because
> the result won't be a highpass it doesn't mean one should not allow
> it: it is assumed that the user knows what he/she/it wants and how to
> get it, the option is there; if not, maybe it's better to do something
> else. Also, while I could verify in parallel the results in wxMaxima
> (slightly more accustomed to), I don't know if the results are,
> indeed, real, so it would be nice if someone who has access to Matlab
> can verify this.

I could help with that, but I know nothing about the intended function
of the file.  I notice you don't have any self-tests at the end. Could
you peek at some other functions to get an idea of example tests to
add? Ideally they would test for both things expected to work and
things expected to error out. Tests that it produces expected trivial
outputs, and maybe a few examples testing Matlab output compatibility
(e.g., examples from their public documentation, help, etc.) I seem to
recall Octav's quadgk function having a decent set of tests you could
use an examples to make your own.

If you sent a set of test matlab commands I could give you the output
to compare against your function.

_______________________________________________
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: I have modified firls.m to include all types of FIRs plus HT and diff

NJank
On Tue, May 16, 2017 at 3:02 PM, Nicholas Jankowski <[hidden email]> wrote:
> On Tue, May 16, 2017 at 1:28 PM, je suis <[hidden email]> wrote:

ok, Matlab help for firls has the following as test inputs to firls:

% Example of a length 31 lowpass filter.
    h=firls(30,[0 .1 .2 .5]*2,[1 1 0 0]);

% Example of a length 45 lowpass differentiator.
    h=firls(44,[0 .3 .4 1],[0 .2 0 0],'differentiator');

% Example of a length 26 type 4 highpass filter.
    h=firls(25,[0 .4 .5 1],[0 0 1 1],'h');

Your code errors out on all three, so there's a bit of a Matlab
compatibility issue.  Here are the Octave outputs:

>> h=lsfir2(30,[0 .1 .2 .5]*2,[1 1 0 0])'
warning: function name 'firls' does not agree with function filename
'C:\Users\nicholas.jankowski\Deskt
op\test\lsfir2.m'
error: For three input arguments, the frequency and weight vectors
must not have a size greater than 2.
error: called from
    lsfir2 at line 110 column 4

>> h=lsfir2(44,[0 .3 .4 1],[0 .2 0 0],'differentiator')'
warning: implicit conversion from numeric to char
warning: called from
    lsfir2 at line 144 column 3
error: lsfir2: operator *: nonconformant arguments (op1 is 45x4, op2 is 28x1)
error: called from
    lsfir2 at line 148 column 3

>> h=lsfir2(25,[0 .4 .5 1],[0 0 1 1],'h')'
warning: implicit conversion from numeric to char
warning: called from
    lsfir2 at line 144 column 3
error: lsfir2: operator *: nonconformant arguments (op1 is 26x4, op2 is 2x1)
error: called from
    lsfir2 at line 148 column 3

For reference, here are the Matlab 2017a outputs.  (matlab outputs h
in as a row vector. display is rotated to column vector for
convenience)

>> format long
>> h=firls(30,[0 .1 .2 .5]*2,[1 1 0 0]);h'
ans =
   0.001271997764274
   0.001759732488755
  -0.000474050123699
  -0.005076557781510
  -0.007306227791006
  -0.001398349858936
   0.011818230551013
   0.020860039466545
   0.010850963397743
  -0.020194873809769
  -0.050342139829269
  -0.042607589983231
   0.027300296130929
   0.144919701151268
   0.257048286517310
   0.303233271924368
   0.257048286517310
   0.144919701151268
   0.027300296130929
  -0.042607589983231
  -0.050342139829269
  -0.020194873809769
   0.010850963397743
   0.020860039466545
   0.011818230551013
  -0.001398349858936
  -0.007306227791006
  -0.005076557781510
  -0.000474050123699
   0.001759732488755
   0.001271997764274
>> h=firls(44,[0 .3 .4 1],[0 .2 0 0],'differentiator');h'
ans =
  -0.000196566503466
  -0.000331109449953
   0.000451753691513
   0.000917072903695
   0.000068070186189
  -0.001369836371904
  -0.001430920748526
   0.000617750911244
   0.002694234626665
   0.001814598671639
  -0.001991887076241
  -0.004577095042312
  -0.001878777817360
   0.004499973574131
   0.007355400044440
   0.001352920990926
  -0.009317947687786
  -0.012579557571954
  -0.000058756249661
   0.022265490917473
   0.036470607272392
   0.028014514414634
                   0
  -0.028014514414634
  -0.036470607272392
  -0.022265490917473
   0.000058756249661
   0.012579557571954
   0.009317947687786
  -0.001352920990926
  -0.007355400044440
  -0.004499973574131
   0.001878777817360
   0.004577095042312
   0.001991887076241
  -0.001814598671639
  -0.002694234626665
  -0.000617750911244
   0.001430920748526
   0.001369836371904
  -0.000068070186189
  -0.000917072903695
  -0.000451753691513
   0.000331109449953
   0.000196566503466

>> h=firls(25,[0 .4 .5 1],[0 0 1 1],'h');h'
ans =
  -0.004214004475093
   0.010158504328219
   0.010872743291530
  -0.012347270923537
  -0.021784061945701
   0.011072745424400
   0.038152378814080
  -0.002854545303196
  -0.063316155566517
  -0.020963226060468
   0.113467223601173
   0.110544374559510
  -0.482742376597236
   0.482742376597236
  -0.110544374559510
  -0.113467223601173
   0.020963226060468
   0.063316155566517
   0.002854545303196
  -0.038152378814080
  -0.011072745424400
   0.021784061945701
   0.012347270923537
  -0.010872743291530
  -0.010158504328219
   0.004214004475093

_______________________________________________
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: I have modified firls.m to include all types of FIRs plus HT and diff

je suis
> Your code errors out on all three, so there's a bit of a Matlab
> compatibility issue.

Thank you for the answers. I realize now that, as I went along with
deciphering the code and the mathematics behind the original file, I
changed the way I wanted the function to work since I was deciphering
it with wxMaxima, and then retranslating it to Octave. So if it's a
mess, don't hold out on the blames.

> ok, Matlab help for firls has the following as test inputs to firls:

Considering these lines, the way this new firls.m should be used is this:

> % Example of a length 31 lowpass filter.
>     h=firls(30,[0 .1 .2 .5]*2,[1 1 0 0]);

h=firls(30,[0 .2 .4 1],[1 1 0 0],[1 1])
h =
   1.27199776427391e-03
   1.75973248875545e-03
  -4.74050123698976e-04
  -5.07655778150969e-03
  -7.30622779100661e-03
  -1.39834985893641e-03
   1.18182305510126e-02
   2.08600394665460e-02
   1.08509633977440e-02
  -2.01948738097687e-02
  -5.03421398292692e-02
  -4.26075899832328e-02
   2.73002961309274e-02
   1.44919701151268e-01
   2.57048286517311e-01
   3.03233271924369e-01
   2.57048286517311e-01
   1.44919701151268e-01
   2.73002961309274e-02
  -4.26075899832328e-02
  -5.03421398292692e-02
  -2.01948738097687e-02
   1.08509633977440e-02
   2.08600394665460e-02
   1.18182305510126e-02
  -1.39834985893641e-03
  -7.30622779100661e-03
  -5.07655778150969e-03
  -4.74050123698976e-04
   1.75973248875545e-03
   1.27199776427391e-03

> % Example of a length 45 lowpass differentiator.
>     h=firls(44,[0 .3 .4 1],[0 .2 0 0],'differentiator');
h=firls(44,[0.3 0.4],[0.2 0],[1 10])
h =
   0.0001935680319982
   0.0005837361716208
   0.0006318704675258
  -0.0000649881094611
  -0.0010909168811554
  -0.0012746155567639
   0.0000830256283033
   0.0019781831198726
   0.0022112619647656
  -0.0002511958142674
  -0.0034193660983087
  -0.0035798846980063
   0.0006828426031918
   0.0058134841226126
   0.0057701601029457
  -0.0015617139368161
  -0.0103256404335874
  -0.0104086300033443
   0.0030055573421858
   0.0228147964389532
   0.0338124265801903
   0.0250410954472427
   0.0000000000000000
  -0.0250410954472427
  -0.0338124265801903
  -0.0228147964389532
  -0.0030055573421858
   0.0104086300033443
   0.0103256404335874
   0.0015617139368161
  -0.0057701601029457
  -0.0058134841226126
  -0.0006828426031918
   0.0035798846980063
   0.0034193660983087
   0.0002511958142674
  -0.0022112619647656
  -0.0019781831198726
  -0.0000830256283033
   0.0012746155567639
   0.0010909168811554
   0.0000649881094611
  -0.0006318704675258
  -0.0005837361716208
  -0.0001935680319982

I think this one differs because of the unimplemented 1/f^2 weighting.
I used [1 10].

> % Example of a length 26 type 4 highpass filter.
>     h=firls(25,[0 .4 .5 1],[0 0 1 1],'h');
h=firls(25,[0 0.4 0.5 1],[0 0 1 1],[1 1],1)
h =
   0.00480437391569623
   0.01050613053429931
  -0.00340207115372934
  -0.01833526416773617
  -0.00228782077798131
   0.02741907886546163
   0.01471652010061388
  -0.03658489209828520
  -0.03855957915035270
   0.04447014544831961
   0.08964774885148744
  -0.04980513866503054
  -0.31257721229656388
   0.55169152299314106
  -0.31257721229656388
  -0.04980513866503054
   0.08964774885148744
   0.04447014544831961
  -0.03855957915035270
  -0.03658489209828520
   0.01471652010061388
   0.02741907886546163
  -0.00228782077798131
  -0.01833526416773617
  -0.00340207115372934
   0.01050613053429931
   0.00480437391569623

This one is different, I'm wrong somewhere in there, unless there is
something hidden under the hood?

At any rate, I will modify the script to match Matlab's choice of
input, even if, honestly, and without infatuation, the way I did it
seems more practical. But compatibility is compatibility. I'll work on
it tomorrow, now it's late here. Should I also modify the output to be
row vector?

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: I have modified firls.m to include all types of FIRs plus HT and diff

Nicholas Jankowski
On May 16, 2017 4:43 PM, "je suis" <[hidden email]> wrote:
> Your code errors out on all three, so there's a bit of a Matlab
> compatibility issue.

Thank you for the answers. I realize now that, as I went along with
deciphering the code and the mathematics behind the original file, I
changed the way I wanted the function to work since I was deciphering
it with wxMaxima, and then retranslating it to Octave. So if it's a
mess, don't hold out on the blames.

> ok, Matlab help for firls has the following as test inputs to firls:

Considering these lines, the way this new firls.m should be used is this:

> % Example of a length 31 lowpass filter.
>     h=firls(30,[0 .1 .2 .5]*2,[1 1 0 0]);

h=firls(30,[0 .2 .4 1],[1 1 0 0],[1 1])
h =
   1.27199776427391e-03
   1.75973248875545e-03
  -4.74050123698976e-04
  -5.07655778150969e-03
  -7.30622779100661e-03
  -1.39834985893641e-03
   1.18182305510126e-02
   2.08600394665460e-02
   1.08509633977440e-02
  -2.01948738097687e-02
  -5.03421398292692e-02
  -4.26075899832328e-02
   2.73002961309274e-02
   1.44919701151268e-01
   2.57048286517311e-01
   3.03233271924369e-01
   2.57048286517311e-01
   1.44919701151268e-01
   2.73002961309274e-02
  -4.26075899832328e-02
  -5.03421398292692e-02
  -2.01948738097687e-02
   1.08509633977440e-02
   2.08600394665460e-02
   1.18182305510126e-02
  -1.39834985893641e-03
  -7.30622779100661e-03
  -5.07655778150969e-03
  -4.74050123698976e-04
   1.75973248875545e-03
   1.27199776427391e-03

> % Example of a length 45 lowpass differentiator.
>     h=firls(44,[0 .3 .4 1],[0 .2 0 0],'differentiator');
h=firls(44,[0.3 0.4],[0.2 0],[1 10])
h =
   0.0001935680319982
   0.0005837361716208
   0.0006318704675258
  -0.0000649881094611
  -0.0010909168811554
  -0.0012746155567639
   0.0000830256283033
   0.0019781831198726
   0.0022112619647656
  -0.0002511958142674
  -0.0034193660983087
  -0.0035798846980063
   0.0006828426031918
   0.0058134841226126
   0.0057701601029457
  -0.0015617139368161
  -0.0103256404335874
  -0.0104086300033443
   0.0030055573421858
   0.0228147964389532
   0.0338124265801903
   0.0250410954472427
   0.0000000000000000
  -0.0250410954472427
  -0.0338124265801903
  -0.0228147964389532
  -0.0030055573421858
   0.0104086300033443
   0.0103256404335874
   0.0015617139368161
  -0.0057701601029457
  -0.0058134841226126
  -0.0006828426031918
   0.0035798846980063
   0.0034193660983087
   0.0002511958142674
  -0.0022112619647656
  -0.0019781831198726
  -0.0000830256283033
   0.0012746155567639
   0.0010909168811554
   0.0000649881094611
  -0.0006318704675258
  -0.0005837361716208
  -0.0001935680319982

I think this one differs because of the unimplemented 1/f^2 weighting.
I used [1 10].

> % Example of a length 26 type 4 highpass filter.
>     h=firls(25,[0 .4 .5 1],[0 0 1 1],'h');
h=firls(25,[0 0.4 0.5 1],[0 0 1 1],[1 1],1)
h =
   0.00480437391569623
   0.01050613053429931
  -0.00340207115372934
  -0.01833526416773617
  -0.00228782077798131
   0.02741907886546163
   0.01471652010061388
  -0.03658489209828520
  -0.03855957915035270
   0.04447014544831961
   0.08964774885148744
  -0.04980513866503054
  -0.31257721229656388
   0.55169152299314106
  -0.31257721229656388
  -0.04980513866503054
   0.08964774885148744
   0.04447014544831961
  -0.03855957915035270
  -0.03658489209828520
   0.01471652010061388
   0.02741907886546163
  -0.00228782077798131
  -0.01833526416773617
  -0.00340207115372934
   0.01050613053429931
   0.00480437391569623

This one is different, I'm wrong somewhere in there, unless there is
something hidden under the hood?

At any rate, I will modify the script to match Matlab's choice of
input, even if, honestly, and without infatuation, the way I did it
seems more practical. But compatibility is compatibility. I'll work on
it tomorrow, now it's late here. Should I also modify the output to be
row vector?

I would say yes. Also, you should test whether orientation of inputs changes orientation if outputs, etc. Main thing though is to get accurate numerical output.  as much Matlab syntax compatibility as practical is good too.

You can post requests for MATLAB output to the list and one of us will usually turn it around pretty quickly.

_______________________________________________
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: I have modified firls.m to include all types of FIRs plus HT and diff

je suis
> I would say yes. Also, you should test whether orientation of inputs
> changes orientation if outputs, etc. Main thing though is to get accurate
> numerical output.  as much Matlab syntax compatibility as practical is good
> too.

Now that I have to modify it according to Matlab's way, I realize the
differentiator part will never match as long as weights are not
specified, since they are using that 1/f^2 by default. Does anyone
know any details about this? What does it mean? Linear weighting means
using integral(K(w)*cos(w*n),w) (or sin(w*n)), so does 1/f^2 mean
integrate(K(w)^2*cos(w*n),w) ? Or sqrt(K(w))? Or the whole integral
squared? Or integrate(integrate(K(w),w)*cos(w*n),w)? One of my initial
goals was to allow linear weighting per band, so that would mean the
first, but 1/f^2...

> You can post requests for MATLAB output to the list and one of us will
> usually turn it around pretty quickly.

Thank you, I will. Also, if I want to use this code in my program
(C++), what do I have to do?

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: I have modified firls.m to include all types of FIRs plus HT and diff

je suis
There's something else that's unclear to me: the function can be either:

firls(N, f, K), or
firls(N, f, A, 'fType'), or
firls(N, f, A, K, 'fType')

Does this mean that, for the 1st case, A takes the value of the
argument K? Which would also meant that, for the 2nd case, K takes the
value or 'fType', which is a string? Similar for the rest? It's a bit
confusing.

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: I have modified firls.m to include all types of FIRs plus HT and diff

Sergei Steshenko
In reply to this post by je suis





________________________________
From: je suis <[hidden email]>
To: [hidden email]
Cc: help <[hidden email]>
Sent: Wednesday, May 17, 2017 11:05 AM
Subject: Re: I have modified firls.m to include all types of FIRs plus HT and    diff



[snip]
Also, if I want to use this code in my program

(C++), what do I have to do?


Vlad


_______________________________________________



Write the code in C++ directly.

FFTW is in "C", so you can use it directly, and there are C++ bindings for linear algebra stuff. There is also templated C++ code for linear algebra, e.g.http://arma.sourceforge.net/ ("Provides high-level syntax (API) deliberately similar to Matlab").


--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: I have modified firls.m to include all types of FIRs plus HT and diff

NJank
In reply to this post by je suis


On May 17, 2017 5:52 AM, "je suis" <[hidden email]> wrote:
There's something else that's unclear to me: the function can be either:

firls(N, f, K), or
firls(N, f, A, 'fType'), or
firls(N, f, A, K, 'fType')

From the matlab docs I see:

Syntax

b= firls(n,f,a)
b = firls(n,f,a,w)
b = firls(n,f,a,'ftype')
b = firls(n,f,a,w,'ftype')


So there's always an A. w is optional.

https://www.mathworks.com/help/signal/ref/firls.html


Also the octave function reference shows:
Function File: b = firls (n, f, a)

Function File: b = firls (n, f, a, w)
https://octave.sourceforge.io/signal/function/firls.html
So  I'm not sure where you got the A is optional syntax.

_______________________________________________
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: I have modified firls.m to include all types of FIRs plus HT and diff

je suis
In reply to this post by Sergei Steshenko
> Write the code in C++ directly.

I used the wrong words. What I meant was: are there any legal issues
if I use the code from firls.m in my program? The program is almost
done, I got here because I needed to verify my implementation and I
saw that Octave's least-squares didn't handle even lengths, or
diffs/HTs, so, one thing led to another...

> So  I'm not sure where you got the A is optional syntax.

I ogled it wrong, but that doesn't change my concern: it looks like W
(or K) can be either W or fType. At first it looked like the
arguments, themselves, don't matter as long as the number is fine, but
then firls(n,f,a,w) and firls(n,f,a,'ftype') both have 4 arguments, so
w can be either a vector or a string. I just find it confusing, and
it's probably a simple check.

What about the 1/f^2 weighting, does anyone know a starting point for
it? I'd rather not just go fishing and hope to catch a boot.

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: I have modified firls.m to include all types of FIRs plus HT and diff

je suis
I modified the script and I could use some testing. The examples
enumerated in there (without the oddities) should cover it, I think.
If you see something out of place, please let me know. Do note,
however, that using default weights for differentiator will not be a
match (I still don't know how to make the 1/f^2 weighting, I'll ask
around).

Also, I tried to input vertical vectors, but somewhere along the way I
need to multiply with a range. I don't know how to make that vertical
or horizontal to match the input. Any help in this matter is
appreciated.

For even length Hilbert transformers, the frequency bands must be
specified as [x 1] (x being any value), otherwise the response around
Nyquist goes wrong.

Vlad

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave

lsfir2.m (9K) Download Attachment
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: I have modified firls.m to include all types of FIRs plus HT and diff

je suis
I'm sorry for the spam, I may have found the solution to the 1/f^2
weighting: https://dsp.stackexchange.com/questions/41044/1-f2-weighting-for-least-squares-fir
(see the last comment).

I can't verify now (tired, late), but it also looks like it can be
extended to regular filters, too, as it would become cos(nw)/w^2,
which results in the last part of the question in the link, or,
according to Wolfram, integrates to
https://www.wolframalpha.com/input/?i=integrate(cos(n*w)%2Fw%5E2,w) .

The question is: is it worth doing it for filters? If yes, it would be
nice to know Octave does extra stuff than Matlab. :-) Still,
sine/cosine integrals do not come cheap in terms of computation, so
will all this be worth it? Does anyone use this sort of weighting, or
is it just a whim?

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: I have modified firls.m to include all types of FIRs plus HT and diff

Sergei Steshenko





________________________________
[snip]

The question is: is it worth doing it for filters? If yes, it would be
nice to know Octave does extra stuff than Matlab. :-) Still,
sine/cosine integrals do not come cheap in terms of computation, so
will all this be worth it? Does anyone use this sort of weighting, or
is it just a whim?


Vlad

_______________________________________________
Help-octave mailing list
[hidden email]

https://lists.gnu.org/mailman/listinfo/help-octave

When I synthesize filters/fit frequency responses, I often use 1/f weighting.

The point is that in acoustics logarithmic frequency display is used, so number of FFT bins per octave is proportional to frequency, so higher frequencies without the 1/f scaling more affect overall fitting than the lower ones.

The 1/f vs 1/f^2 is, I guess, amplitude vs energy - energy is proportional the the square of amplitude.

--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: I have modified firls.m to include all types of FIRs plus HT and diff

je suis
> When I synthesize filters/fit frequency responses, I often use 1/f
> weighting.
>
> The point is that in acoustics logarithmic frequency display is used, so
> number of FFT bins per octave is proportional to frequency, so higher
> frequencies without the 1/f scaling more affect overall fitting than the
> lower ones.
>
> The 1/f vs 1/f^2 is, I guess, amplitude vs energy - energy is proportional
> the the square of amplitude.

I should've seen this coming, even more so when I used audio-related
stuff a few times. Well, I'll work on it and try to add to both
differentiators and filters, but I can already smell the burn on the
CPU.

In the meantime, except for the "default" weighting for
differentiators, did the new script behave well enough?

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: I have modified firls.m to include all types of FIRs plus HT and diff

NJank


On May 18, 2017 2:50 AM, "je suis" <[hidden email]> wrote:
> When I synthesize filters/fit frequency responses, I often use 1/f
> weighting.
>
> The point is that in acoustics logarithmic frequency display is used, so
> number of FFT bins per octave is proportional to frequency, so higher
> frequencies without the 1/f scaling more affect overall fitting than the
> lower ones.
>
> The 1/f vs 1/f^2 is, I guess, amplitude vs energy - energy is proportional
> the the square of amplitude.

I should've seen this coming, even more so when I used audio-related
stuff a few times. Well, I'll work on it and try to add to both
differentiators and filters, but I can already smell the burn on the
CPU.

In the meantime, except for the "default" weighting for
differentiators, did the new script behave well enough?

Vlad

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave


_______________________________________________
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: I have modified firls.m to include all types of FIRs plus HT and diff

NJank
In reply to this post by je suis


On May 17, 2017 12:18 PM, "je suis" <[hidden email]> wrote:
I modified the script and I could use some testing. The examples
enumerated in there (without the oddities) should cover it, I think.
If you see something out of place, please let me know. Do note,
however, that using default weights for differentiator will not be a
match (I still don't know how to make the 1/f^2 weighting, I'll ask
around).

Also, I tried to input vertical vectors, but somewhere along the way I
need to multiply with a range. I don't know how to make that vertical
or horizontal to match the input. Any help in this matter is
appreciated.

If I understood your previous email correctly, and the issue you mention here it sounds like you need to do some input checking and preconditioning before processing. I'll try to take a look at your code later today, but for the case of the third input either being a vector or a string you would want to have octave check the class of the input for that variable prior to doing any computation. Since there is a variable number of possible inputs, you will want to use 'varargin' to collect them and then do type checking to decide what to do with them and how to run the code. Again, the quadgk function is a good example of this as it has to deal with both a variable number of possible inputs and mixed numerical and string inputs. There is also now a supported inputParser class that supposedly is a more elegant way to do this, but I'm not familiar with it so I just do it the old-fashioned way.

Then for any vectors, if they don't come in the orientation you want, you can have octave change them to either a row or column inside the function as you desire. I find the fastest way to do this is either

A = A(:)

for column output or

A = A(:)'

for row output


_______________________________________________
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: I have modified firls.m to include all types of FIRs plus HT and diff

NJank
On May 18, 2017 7:23 AM, "Nicholas Jankowski" <[hidden email]> wrote:


On May 17, 2017 12:18 PM, "je suis" <[hidden email]> wrote:



<snip>

looking at the code, you're first line

function h = firls(N, F, A, K, ftype)

Should be:

function h = firls(N, F, A, varargin)

That will let it accept any number of arguments after A, and store them in a cell array named varargin.  Then, a code block based in the number and type of varagin elements can determine how to run your code

 Something with the form below might work:

switch numel(nargin)

case 3
<set defaults>

case 4
<Check if #4 is a weight vector or a string, set values, like K = varargin{1}, and defaults for the rest>

case 5
<verify the order is right, assign K= varargin{1}, ftype= varargin{2}>

otherwise
  print_usage()

Then, check that K is valid, ftype matches, set vectors to columns or rows, etc.



_______________________________________________
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: I have modified firls.m to include all types of FIRs plus HT and diff

je suis
Sorry for the delay, I used some more scripts and managed to make it
with "varargin", and also seems to work. I'm surprised the way I did
it before even worked. The input vectors should not matter anymore
now, the output is a row vector.

But now I got stuck a bit at the part with the derivation. I got the
generic results of the integral shown in the block comments, in the
script, but the results are awful. Could any of you please look at it?

In short, what I think: a*Si(x) means a is the slope
(A[k+1]-A[k])/(F[k+1]-F[k]) and Si(x) is the sine integral, defined as
imag(expint(1i*x))+pi/2. This is used as a difference. Then
b*[Ci(x)-sin(x)/x] means b is the intercept, thus A[k]. Where am I
wrong?

The test I used is lsfir2(11, [0 0.3 0.4 1], [0 1 0 0], 'd'). This,
according to Matlab's page, uses the default 1/f^2 weighting. I don't
know what output is supposed to be, but certainly not the one that
comes out now.

Vlad

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave

lsfir2.m (15K) Download Attachment
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: I have modified firls.m to include all types of FIRs plus HT and diff

Ozzy Lash


On Sat, May 20, 2017 at 5:01 AM, je suis <[hidden email]> wrote:
Sorry for the delay, I used some more scripts and managed to make it
with "varargin", and also seems to work. I'm surprised the way I did
it before even worked. The input vectors should not matter anymore
now, the output is a row vector.

But now I got stuck a bit at the part with the derivation. I got the
generic results of the integral shown in the block comments, in the
script, but the results are awful. Could any of you please look at it?

In short, what I think: a*Si(x) means a is the slope
(A[k+1]-A[k])/(F[k+1]-F[k]) and Si(x) is the sine integral, defined as
imag(expint(1i*x))+pi/2. This is used as a difference. Then
b*[Ci(x)-sin(x)/x] means b is the intercept, thus A[k]. Where am I
wrong?

The test I used is lsfir2(11, [0 0.3 0.4 1], [0 1 0 0], 'd'). This,
according to Matlab's page, uses the default 1/f^2 weighting. I don't
know what output is supposed to be, but certainly not the one that
comes out now.

Vlad

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave


I can't really help much, but looking at your function, shouldn't the lines:

   q = (-pi*n'.*imag(expint(1i*n'*w(2:end)) + pi2) - ...

    q = (-pi*n'.*imag(expint(1i*n'*w) + pi2) - cos(n'*w)./F)*K'

have parens around the imag...pi2 section?  That is:

   q = (-pi*n'.*(imag(expint(1i*n'*w(2:end)) + pi2)) - ...

    q = (-pi*n'.*(imag(expint(1i*n'*w) + pi2)) - cos(n'*w)./F)*K'

since the pi/2 is part of the Si(x), I think it should be that way.

_______________________________________________
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: I have modified firls.m to include all types of FIRs plus HT and diff

je suis
> I can't really help much, but looking at your function, shouldn't the
> lines:
>    q = (-pi*n'.*imag(expint(1i*n'*w(2:end)) + pi2) - ...
>     q = (-pi*n'.*imag(expint(1i*n'*w) + pi2) - cos(n'*w)./F)*K'
> have parens around the imag...pi2 section?  That is:
>    q = (-pi*n'.*(imag(expint(1i*n'*w(2:end)) + pi2)) - ...
>     q = (-pi*n'.*(imag(expint(1i*n'*w) + pi2)) - cos(n'*w)./F)*K'
> since the pi/2 is part of the Si(x), I think it should be that way.

Thank you for taking a look, that was a mistake, corrected, but also
my magnificent mathematical skills are to blame. When doing linear
interpolation, I know that the slope is (A2 - A1)/(f2 - f1) and that
the whole comes as slope*(x - f1) + A1, but what I did was consider
the (x - f1) as a whole term, thus, for the continuous version a*x+b,
I made b = A1, when, in fact, it's intercept = -slope*f1 + A1, and
then it comes as slope*x+intercept.

Now I am getting some results but, rather than keeping on uploading
versions, maybe it would be less time consuming if someone would
simply verify the inputs below.

Vlad

---

firls(10, [0, 200, 300, 512]/512, [0, 1, 0, 0], 'd'), plot(abs(fft(h,
1024))(1:512))
0.04235469148728355
-0.08201413012272996
-0.0618119893535477
0.2018661528543597
0.3083554734960205
0
-0.3083554734960205
-0.2018661528543597
0.0618119893535477
0.08201413012272996
-0.04235469148728355

firls(11, [0, 200, 300, 512]/512, [0, 1, 0, 0], 'd'), plot(abs(fft(h,
1024))(1:512))
0.03732594529679734
-0.03433612793625151
-0.09696356417897994
0.07065147268158478
0.2931623606010482
0.1793606524572318
-0.1793606524572318
-0.2931623606010482
-0.07065147268158478
0.09696356417897994
0.03433612793625151
-0.03732594529679734

firls(30, [0, 200, 300, 512]/512, [0, 1, 0, 0], 'd'), plot(abs(fft(h,
1024))(1:512))
-0.001269277091273668
0.003198392108892278
0.001285061684276556
-0.008451971441951756
-0.00162882496078165
0.01772022558321762
0.002475017539728341
-0.03311523898587065
-0.004926273092083377
0.05871977740512246
0.01247729136542475
-0.1054747090506596
-0.04213126736540573
0.2142790732028212
0.2893569317389009
0
-0.2893569317389009
-0.2142790732028212
0.04213126736540573
0.1054747090506596
-0.01247729136542475
-0.05871977740512246
0.004926273092083377
0.03311523898587065
-0.002475017539728341
-0.01772022558321762
0.00162882496078165
0.008451971441951756
-0.001285061684276556
-0.003198392108892278
0.001269277091273668

firls(31, [0, 200, 300, 512]/512, [0, 1, 0, 0], 'd'), plot(abs(fft(h,
1024))(1:512))
-0.00148394238352223
0.002001300556950358
0.003704236399439953
-0.005504263029899814
-0.007962074902567284
0.01147150277313669
0.01551725107452073
-0.02078713878644933
-0.02870851767322025
0.03474477700849787
0.05295249630879661
-0.05567298153479006
-0.1063785931609718
0.08391965177763228
0.2962965843028105
0.1745324683632976
-0.1745324683632976
-0.2962965843028105
-0.08391965177763228
0.1063785931609718
0.05567298153479006
-0.05295249630879661
-0.03474477700849787
0.02870851767322025
0.02078713878644933
-0.01551725107452073
-0.01147150277313669
0.007962074902567284
0.005504263029899814
-0.003704236399439953
-0.002001300556950358
0.00148394238352223

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
123
Loading...