patchGNU Octave - Patches: patch #8535, [signal] firpm (new remez)

 
 

patch #8535: [signal] firpm (new remez)

Submitted by:  Rob Sykes <robsykes>
Submitted on:  Sun 14 Sep 2014 06:21:27 AM UTC  
 
Category: NonePriority: 5 - Normal
Status: NonePrivacy: Public
Assigned to: Mike Miller <mtmiller>Open/Closed: Open

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)

Thu 28 Jul 2016 07:47:06 AM UTC, comment #15:

It worked, thanks!

I tried your firpm on a 8000 coefficient filter and it works very well (and fast).

I hope this stuff will be added soon to the signal package.

Roland Baudin <roland65>
Wed 27 Jul 2016 04:29:07 PM UTC, comment #14:

Please try the following:

if (arg < nargin)
{
std::ostringstream os;
- args(arg).print (os, true);
+ octave_value tmp(args(arg)); /* Need tmp since print is no longer const. */
+ tmp.print (os, true);
std::string s(os.str ());
size_t p = s.find_last_not_of (" \t\r\n");
if (p != std::string::npos)

Rob Sykes <robsykes>
Wed 27 Jul 2016 12:03:11 PM UTC, comment #13:

I tried to compile firpm with g++ 5.4.0 (in Ubuntu 16.04) and I got the following error message:

firpm.cc: In function ‘octave_value_list Ffirpm(const octave_value_list&, int)’:
firpm.cc:392:32: error: passing ‘const octave_value’ as ‘this’ argument discards qualifiers [-fpermissive]
args(arg).print (os, true);
^
In file included from /usr/include/octave-4.0.3/octave/../octave/oct-obj.h:33:0,
from /usr/include/octave-4.0.3/octave/../octave/ov-fcn.h:31,
from /usr/include/octave-4.0.3/octave/../octave/ov-builtin.h:28,
from /usr/include/octave-4.0.3/octave/../octave/defun-int.h:28,
from /usr/include/octave-4.0.3/octave/../octave/defun-dld.h:30,
from /usr/include/octave-4.0.3/octave/../octave/oct.h:36,
from firpm.cc:21:
/usr/include/octave-4.0.3/octave/../octave/ov.h:1030:8: note: in call to ‘void octave_value::print(std::ostream&, bool)’
void print (std::ostream& os, bool pr_as_read_syntax = false)

Any idea about fixing thi issue?

Roland Baudin <roland65>
Sun 11 Oct 2015 07:38:02 AM UTC, comment #12:
  • added an example of precise minimum order determination

(file #35137)

Rob Sykes <robsykes>
Thu 21 May 2015 06:31:13 AM UTC, comment #11:

Changes:

  • Can now use openmp (scales well)
  • Added test for the 'persistence' property

(file #34072)

Rob Sykes <robsykes>
Mon 04 May 2015 01:31:25 PM UTC, comment #10:

Changes between file #33055 and file #33895:

mmfir.c

  • Fix negative invalid DC/nyquist amplitude not being detected
  • Fix incorrect sign of non-differentiator type III/IV filter
  • Move allocation for filter coefs client-side
  • Unembed from firpm.cc & move to cruft

firpm.cc:

  • Use Gnu-style white-space
  • Cleaner/extensible syntax using properties
  • Support Matlab's lgrid syntax and emulate its effect
  • More tests
  • Improved demos, & documentation

New companion functions:

  • firpmord.m
  • remezord.m (still in Matlab)

firpm.cc, *.m tested with Octave 3.8.1 on x86_64 GNU/Linux
Also, mmfir.c tested with MSVS 2008 on Win32

Rob Sykes <robsykes>
Sun 22 Feb 2015 05:34:20 PM UTC, comment #9:

The ideal outcome would be if the files were rewritten to conform to Octave's coding style and use STL and Octave library functions where appropriate :)

The same applies to ultrwin really, I would like to see that use Octave's coding style and library functions as well if possible.

If that is a non-starter, I'll think about other ways we can work with this. Maybe add the code to an internal library "cruft" directory and link the oct-files against that, like the "cruft" directory that Octave maintains with various numerical libraries stuffed inside.

Mike Miller <mtmiller>
Project AdministratorIn charge of this item.
Sun 22 Feb 2015 09:23:00 AM UTC, comment #8:

Mike, yes, the submitted .cc file comprises three components: mmfir.c (the filter designing code), mmfir.h (it's interface), and the Octave 'firpm' wrapper. Note that, for maximum portability, mmfir.c/h are written in the C language, so use of std::max, for example, is not appropriate there. (The situation is very similar to that of Jake Janovetz's C implementation and its Octave remez wrapper.) The C part is, in effect, library code, but I know from previous experience that Octave prefers the simplicity of a single source file. If you think it's better to do so, the C code could be given a separate home as a library on SourceForge or Github, etc.

Rob Sykes <robsykes>
Thu 12 Feb 2015 06:02:21 PM UTC, comment #7:

Thanks Laurent for evaluating this patch and providing feedback. I still haven't had the time to look at this in much depth yet.

At a glance, though, I would like it if C++ submissions to the signal package conform to the Octave and GNU coding standards. It looks to me like you want to use this code elsewhere, and have a script to combine it into one file for compilation as a .oct file. That's fine, but code distributed with the signal package should be consistent with the coding style of Octave. If it isn't submitted that way, it will be reformatted later eventually, and won't be comparable to your original source.

It would also be preferable to use standard library and Octave library functions where possible, e.g. use std::max() and xround() instead of defining your own :)

Mike Miller <mtmiller>
Project AdministratorIn charge of this item.
Wed 11 Feb 2015 09:28:29 PM UTC, comment #6:

It seems that, though Matlab documentation says 1/f, it actually applies 2/f weighting, so I've modified it here to do the same. Also, the problematic max macro has now been eliminated.

(file #33055)

Rob Sykes <robsykes>
Sun 08 Feb 2015 01:48:17 PM UTC, comment #5:

The update looks pretty good. I was a bit surprised at the relative weight of the passband vs. stopband.
Not much of an issue since the user can provide weights to obtained whatever passband/stopband tradeoff they desire, but to get the same results as remez (and some literature results) without playing with the weights, you'd actually have to tweak W(f) to divide by (f/2) instead of f.
I guess the question is what does Matlab generate in that case (if we want to match Matlab's output for the same arguments). I unfortunately don't have Matlab to check.

Strangely enough the issue with VS (which I reproduced with VS2010, 2012 and 2013 express) messes up the second argument to the max macro, and it's still an issue with the renamed max_ variable
(see [this microsoft connect issue](https://connect.microsoft.com/VisualStudio/feedback/details/1101733/incorrect-nested-macro-expansion)).

Laurent Goulet <laurentg>
Thu 08 Jan 2015 01:01:19 PM UTC, comment #4:

Thanks for this feedback. There's an updated version addressing differentiator stop-band. I couldn't repeat the VS issue with the version I have here but have renamed a variable in the hope that it might help.

(file #32784)

Rob Sykes <robsykes>
Wed 07 Jan 2015 01:46:33 AM UTC, comment #3:

This looks like a good alternative to remez.
I did however notice a few minor things:

Differentiator design from the original PM program used a 1/f weight only in the passband. This typically results in a lower relative error in the passband.
For example,
with the patch for Bug #38134
yields a filter with relative error of 0.074% in the passband and an absolute equiripple error of 0.0074 in the stopband
(which corresponds to some textbook examples of PM program execution).
firpm on the other hand yields a filter with a relative error of 0.2%
in the passband and a maximum absolute error (not equiripple) of 0.01 in the stopband.

Also, the code seems to trigger a preprocessing bug in Visual Studio with the max macro arguments at lines 620 & 643 of file #32118
(a simple workaround being to use a temporary, eg. tmp=fabs(E(i)),max=max(max,tmp)).
I realize that this package is typically not built with Visual Studio so that should generally not be an issue.

Laurent Goulet <laurentg>
Tue 16 Sep 2014 08:02:18 PM UTC, comment #2:

Yes, updated patch includes the wrapper (also added a few tests for firpm argument errors).

(file #32117)

Rob Sykes <robsykes>
Mon 15 Sep 2014 03:10:07 PM UTC, comment #1:

Thanks for your contribution! I will definitely take a look at incorporating this into the signal package, the remez function definitely has some problems as it currently stands.

Will it also be possible to rewrite remez as a thin wrapper m-file function calling firpm?

Mike Miller <mtmiller>
Project AdministratorIn charge of this item.
Sun 14 Sep 2014 06:21:27 AM UTC, original submission:

This is a new implementation of McClellan-Parks-Rabiner method.
Compared to the current Octave 'remez', it has:

New features:

  • Weight per band-edge.
  • Simple syntax for amplitude per band.
  • Can increase number of iterations.
  • Returns minimax error.

Performance improvements:

  • Allows longer filters.
  • Several times faster.
  • Produces more accurate filters.

Other improvements:

  • Improved help.
  • Demos.

Fewer bugs AFAICT:
The following filters fail to converge with remez, presumably because of bugs:

b=remez(40, [0 .25 .3 .6 .65 1],[0 0 1 1 0 0]);
b=remez(40, [0 2.5/pi], [0 2.5], 'differentiator');
b=remez(31, [0 .5 .7 1], [0 0 1 1], 'hilbert');
b=remez(120, [0 .5 .55 .6 .601 .75 .751 1], [1 1 0 0 0 0 0 0], [1 8 32 64]);

For a speed & accuracy improvements example, try:

tic; b=remez(1500, [0 .0267 .033 1], [1 1 0 0]); toc;[h,f]=freqz(b,1,2^17); plot(f/pi, 20*log10(abs(h))); axis([0.032 .045 -81.5 -80.5]);

then the same with firpm.

Rob Sykes <robsykes>

 

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

Attach File(s):
   
   
Comment:
   

Attached Files
file #35137:  firpm.tar.xz added by robsykes (20KiB - application/x-xz - adds an example of precise minimum order determination)
file #34072:  firpm.tar.xz added by robsykes (20KiB - application/x-xz - adds openmp, persistence test)
file #33895:  firpm.tar.xz added by robsykes (20KiB - application/x-xz - add firpmord; fix negative invalid DC/nyquist amplitude not being detected; restore differentiator sign; improve tests, demos, doc.; separate firpm/mmfir)
file #33792:  firpm.tar.xz added by robsykes (15KiB - application/x-xz - fix sign of type III/IV filter; move allocation for filter coefs client-side; Octave wrapper: cleaner/extensible syntax using properties; Gnu-style white-space; doc. & example improvements)
file #33055:  firpm.tar.xz added by robsykes (14KiB - application/x-xz - Updates in light of comment # 5; minor clean-ups)
file #32784:  firpm.tar.xz added by robsykes (14KiB - application/x-xz - Updates in light of comment # 3)
file #32118:  firpm.tar.xz added by robsykes (14KiB - application/x-xz - fix bug in test in file #32117)
file #32117:  firpm.tar.xz added by robsykes (14KiB - application/x-xz - with remez wrapper and tests for firpm argument errors)
file #32091:  firpm.cc.xz added by robsykes (13KiB - application/x-xz)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -unavailable- added by roland65 (Posted a comment)
  • -unavailable- added by laurentg (Posted a comment)
  • -unavailable- added by mtmiller (Posted a comment)
  • -unavailable- added by robsykes (Submitted the item)
  •  

    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 logged-in users can vote.

     

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

     

     

    Follow 10 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Sun 11 Oct 2015 07:38:02 AM UTCrobsykesAttached File-=>Added firpm.tar.xz, #35137
    Thu 21 May 2015 06:31:13 AM UTCrobsykesAttached File-=>Added firpm.tar.xz, #34072
    Mon 04 May 2015 01:29:30 PM UTCrobsykesAttached File-=>Added firpm.tar.xz, #33895
    Sun 26 Apr 2015 08:18:46 PM UTCrobsykesAttached File-=>Added firpm.tar.xz, #33792
    Wed 11 Feb 2015 09:28:29 PM UTCrobsykesAttached File-=>Added firpm.tar.xz, #33055
    Thu 08 Jan 2015 01:01:19 PM UTCrobsykesAttached File-=>Added firpm.tar.xz, #32784
    Tue 16 Sep 2014 08:42:33 PM UTCrobsykesAttached File-=>Added firpm.tar.xz, #32118
    Tue 16 Sep 2014 08:02:18 PM UTCrobsykesAttached File-=>Added firpm.tar.xz, #32117
    Mon 15 Sep 2014 03:10:07 PM UTCmtmillerAssigned toNone=>mtmiller
    Sun 14 Sep 2014 06:21:27 AM UTCrobsykesAttached File-=>Added firpm.cc.xz, #32091

    Back to the top


    Powered by Savane 3.1-cleanup1