bugGNU Octave - Bugs: bug #51310, [octave forge] (signal) firls.m...

 
 

bug #51310: [octave forge] (signal) firls.m modification to include all 4 FIR types, Hilbert transformer, and differentiator

Submitted by:  Vlad <archbugaboo>
Submitted on:  Sun 25 Jun 2017 07:00:13 PM UTC  
 
Category: Octave Forge PackageSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Feature Request
Status: NoneAssigned to: None
Originator Name: Open/Closed: Open
Release: 4.2.1Operating System: Any

Add a New Comment(Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

(Jump to the original submission Jump to the original submission)

Sat 01 Jul 2017 11:11:51 AM UTC, comment #8:

> Ok. Was just thinking that if the output doesn't match MATLAB, someone in the future will probably report incompatibility. Some third party to compare against to justify 'goodness' of answer would be worthwhile


Well, for what it's worth, I'll upload two of the results to match these two (for now):

firls (42, [0, 0.3, 0.4, 1], [0, 0.3, 0, 0]*pi, 'd')
firls (43, [0, 0.3, 0.4, 1], [0, 0.3, 0, 0]*pi, [30, 1], 'd')

The rest of the examples are the same.

I also changed the script to allow a 'MatlabCompat' argument that takes 'no' or 'n', meaning it will not be compatible with Matlab's firls.m. If unspecified, it defaults to compatibility.

(file #41095, file #41096)

Vlad <archbugaboo>
Fri 30 Jun 2017 04:30:51 PM UTC, comment #7:

Ok. Was just thinking that if the output doesn't match MATLAB, someone in the future will probably report incompatibility. Some third party to compare against to justify 'goodness' of answer would be worthwhile

Also, I'd recommend a more informative name than 'x'. 'prevent_order_increment' or... something? And I haven't looked but help and %!test should be updated to cover the new input if they haven't already.

Nicholas Jankowski <nrjank>
Fri 30 Jun 2017 03:50:28 PM UTC, comment #6:

> I think a posting a couple maxima comparisons would be instructive.


It might not be so, since the code in wxMaxima follows the one in Octave (actually the other way around, but it doesn't matter). It would only matter if wxMaxima had Matlab's implementation, which it doesn't. The only thing it will prove is the accuracy of expint().

I modified the script to allow an extra argument, 'x'. If left alone (no extra argument), the script will increment N for type II and IV FIRs with non-zero amplitude@Nyquist. If added, it will not increment it. The behaviour now matches Matlab's, while being able to be changed as desired. The changed version can also be found at the github address (in the OP), but I'll also upload it.

Vlad <archbugaboo>
Thu 29 Jun 2017 01:43:02 PM UTC, comment #5:

I think a posting a couple maxima comparisons would be instructive.

Nicholas Jankowski <nrjank>
Thu 29 Jun 2017 05:30:19 AM UTC, comment #4:

The more I think of it, the better suited this seems for the increment N part: adding an extra parameter, say 'x', that, when added, makes the behaviour as originally intended, otherwise (no argument, default) mimicks Matlab. Later this day I'll make an example with what I mean.

About (2), I see that the results seem to get worse as N increases, at least based on the existent examples in the discussion list. For example, according to the results in https://lists.gnu.org/archive/html/help-octave/2017-05/msg00179.html and the followup https://lists.gnu.org/archive/html/help-octave/2017-05/msg00180.html, if one considers "h" to be Octave's output and "g" to be Matlab's, then plotting this should give an idea:

plot(diff(abs(fft(h,1024)(1:512))),"",diff(abs(fft(g,1024)(1:512))))

Zooming in on the passband region it can be seen that, as the order increases, the 1/f^2 weighting is clearly consistent in for "h", but not for "g". Going further with https://lists.gnu.org/archive/html/help-octave/2017-05/msg00138.html (the 2nd example, order 44) and the followup https://lists.gnu.org/archive/html/help-octave/2017-05/msg00139.html, the differences are more visible, and going even further with https://lists.gnu.org/archive/html/help-octave/2017-06/msg00259.html (the "firls_expected_output.txt" attachement) it can be seen that the passband for "g" gets more and more distorted, but the general magnitude is close to the requirements.

This is why I think Matlab must use some very small number for F(1)=0, probably 1e-6 or similar, which, when used for cos(1e-6)/1e-6 results in a six orders of magnitude number, which is then used to build the Q matrix which, in turn, is used for a=Q\b. This must have a negative effect on precision.

The results I had with wxMaxima were numerical confirmations, since their expintegral_e1(x) is accurate and can be used for bigfloats, as well.

Vlad <archbugaboo>
Wed 28 Jun 2017 02:38:05 PM UTC, comment #3:

Mike? thoughts on either (2) or (3)?

Unless we really want to force compatibility here, my suggestion on (3) is to have the function produce a warning when the non-compatibile event occurs, document it in the help text, and move on. One alternative would be would be to provide an optional compatibility flag so that it forces compatible behavior. (i.e., if (matlab_compatible) n++ )

See the help text for the currently released version of cov [1] for an example of such documentation.

Then, for (2) you said you were benchmarking against Maxima before. Can you provide a couple examples here showing that what you're referring to on your code produces a more accurate filter? I'm not disputing this, but having it documented here would be enough justification in my opinion to leave it as is and move on. If someone later has a better idea for the algorithm they could always propose an update. What's more important is to be confident your version doesn't produce any wildly incorrect answers (for odd inputs, corner cases, etc.). I think another note in the help would be useful stating that because of algorithmic differences the output may not exactly match Matlab's but still meets "the definition of 'good' according to ??_".

Then if format/function are good, seems it's about ready to propose rolling it up into a patch for the signal package.

[1] https://octave.sourceforge.io/octave/function/cov.html

Nicholas Jankowski <nrjank>
Mon 26 Jun 2017 05:49:19 AM UTC, comment #2:

> So, apart from the accuracy issues in (2) and the compatibility decision in (3), is there anything else that needs to be done/fixed for it to be completely functional?


No, unless some bugs crept up the code, the script is functional.

> Is there another set of data we can evaluate (2) against? Maybe a symbolic routine for so.e simpler cases to which we can compare yours and matlabs outputs?


I'm afraid the results will always come visibly close, but not quite. I think, at this point, the question is: would it be preferred to mimick a wrong answer, or improve it? (because, based on my tests, I dare say Matlab is wrong, though not by much)

> Would maybe be worth get the package mmaintainers opinion on (3)


That's why I said it's a personal view. However, the results are correct for types I and III, which would be the resulting FIRs after the increment, so the results are the same, it's just that the script doesn't forcefully choose for you. Consider the educational point of view: you're a teacher and want to show your students what a type III or IV lowpass FIR looks like; could you do it with Matlab's script?

At any rate, if this contravenes so much with the intended way the script is meant to be working, then I'll change it. Maybe add an extra parameter?

Vlad <archbugaboo>
Sun 25 Jun 2017 11:23:11 PM UTC, comment #1:

So, apart from the accuracy issues in (2) and the compatibility decision in (3), is there anything else that needs to be done/fixed for it to be completely functional?

Is there another set of data we can evaluate (2) against? Maybe a symbolic routine for so.e simpler cases to which we can compare yours and matlabs outputs?

Would maybe be worth get the package mmaintainers opinion on (3)

Nicholas Jankowski <nrjank>
Sun 25 Jun 2017 07:00:13 PM UTC, original submission:

Following the discussions in the Octave mailing lists, namely https://lists.gnu.org/archive/html/help-octave/2017-05/msg00136.html and https://lists.gnu.org/archive/html/help-octave/2017-06/msg00058.html, I managed to get this version for firls.m that allows creating all four types of FIR, plus Hilbert transformers and differentiators, the latter allowing either weighted bands or a special 1/f^2 weighting, all according to Matlab's description. The script has some differences between Matlab's, as follows:

1. It uses the more recent advances in Octave's builtin expint() (see http://hg.savannah.gnu.org/hgweb/octave/file/tip/scripts/specfun/expint.m), which corrects the previous function for completely imaginary arguments.

2. The derivation for the least-squares FIR is different than those mentioned in Matlab's online help (https://nl.mathworks.com/help/signal/ref/firls.html#bqih5y1-1 and https://nl.mathworks.com/help/signal/ref/firls.html#bqih5y9-1), and also the one in the original firls.m script (Ivan Selesnick's paper), and it involves sine and cosine integrals, and cos(x)/x terms. Ci(0) and cos(0)/0 result in Inf, and Matlab must use some sort of avoiding (most probably a small number for the frequency at zero) in orcer to circumvent the Inf problem. In my derivation, I considered their values as zero, in the same way type III FIRs require the impulse response be zero at h[M]. As a consequence, after verifying the resposnes side by side, and the derivatives of their FFTs, I conclude that my numbers are closer to reality, thus the results between Octave and Matlab are different, but not necessarily wrong.

3. This is more of a personal view. Types II and IV FIRs have a builtin zero at Nyquist, while types III an additional zero at DC, so having a type II highpass FIR will come with inherent downsides. But this doesn't mean that one cannot build a type II highpass FIR, since the books clearly warn you that making one will not yield what you would expect, but you can create one if you so wish. Because of this, my script does not increment the order behind your back in order to give what the script might think I want to see, instead, it lets you create your type II highpass FIR and enjoy the crippled magnitude, if that is what you wish. In short, it does not think it knows your thoughts better, but considers you know what you're doing (since you're using the scrpt), and you're a grown up that knows how to handle the responsabilities. If not, whatever would you be doing here in the first place?

The original script did its job, but only allowed the creation of type I FIRs.

As for the uploading of the script, I'm in doubt. I uploaded the script, but it can also be found at https://github.com./archbugaboo/Octave_signal.

Vlad <archbugaboo>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach File(s):
   
   
Comment:
   

Attached Files
file #41095:  firls.m added by archbugaboo (39KiB - text/x-objcsrc - changed argument 'x' to something hopefully clearer; two comparative results in wxMaxima.)
file #41096:  wxMaxima_test.txt added by archbugaboo (2KiB - text/plain - changed argument 'x' to something hopefully clearer; two comparative results in wxMaxima.)
file #41083:  firls.m added by archbugaboo (37KiB - text/x-objcsrc - modified for extra argument, 'x')
file #41044:  firls.m added by archbugaboo (34KiB - text/x-objcsrc)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -unavailable- added by jwe (Updated the item)
  • -unavailable- added by nrjank (Posted a comment)
  • -unavailable- added by nrjank (Signal package maintainer)
  • -unavailable- added by archbugaboo (Submitted the item)
  • -unavailable- added by archbugaboo
  •  

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 7 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Sat 12 Aug 2017 12:47:13 AM UTCjweSummaryfirls.m modification to include all 4 FIR types, Hilbert transformer, and differentiator=>[octave forge] (signal) firls.m modification to include all 4 FIR types, Hilbert transformer, and differentiator
    Sat 01 Jul 2017 11:11:51 AM UTCarchbugabooAttached File-=>Added firls.m, #41095
      Attached File-=>Added wxMaxima_test.txt, #41096
    Fri 30 Jun 2017 03:51:09 PM UTCarchbugabooAttached File-=>Added firls.m, #41083
    Sun 25 Jun 2017 11:23:11 PM UTCnrjankCarbon-Copy-=>Added mtmiller
    Sun 25 Jun 2017 07:00:13 PM UTCarchbugabooAttached File-=>Added firls.m, #41044
      Carbon-Copy-=>Added nrjank

    Back to the top


    Powered by Savane 3.1-cleanup1