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

 
 

patch #8535: [octave forge] (signal) firpm: new function, replacement for remez

Submitter:  Rob Sykes <robsykes>
Submitted:  Sun 14 Sep 2014 06:21:27 AM UTC
   
 
Category:  Forge : new function Priority:  5 - Normal
Status:  Done Privacy:  Public
Assigned to:  lostbard Open/Closed:  Closed
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 17 May 2023 05:53:44 PM UTC, comment #21: 

v1.4.4 released - closing ticket

John Donoghue <lostbard>
Group Member
Wed 10 May 2023 12:13:27 PM UTC, comment #20: 

done

John Donoghue <lostbard>
Group Member
Wed 10 May 2023 04:24:28 AM UTC, comment #19: 

Hello John, thanks for having a look at this patch.  Would it be possible to use file #51196 for the signal repo?  It contains several improvements over the earlier submission.

Rob Sykes <robsykes>
Mon 08 May 2023 07:26:57 PM UTC, comment #18: 

Added to the signal repo

John Donoghue <lostbard>
Group Member
Mon 31 May 2021 08:57:09 AM UTC, comment #17: 

With the addition of respfn and res, this is now feature-complete w.r.t. Matlab AFAIK.  Code-style status is: .m & .cc files use Octave style; mmfir.[ch] use K&R style but reside is a sub-directory (to be considered perhaps, as a 3rd-party library).  If this sufficient, or would it help to run mmfir.[ch] through gnu indent?

Rob Sykes <robsykes>
Thu 18 Jan 2018 09:32:50 PM UTC, comment #16: 

Updated per Octave API changes

(file #42968)

Rob Sykes <robsykes>
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>
Group Member
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>
Group Member
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,

remez(59, [0 0.2 0.3 1], [0 1 0 0], 'differentiator')

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>
Group Member
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 Files:
   
   
Comment:
   

Attached Files
file #51196:  firpm.tar.xz added by robsykes (23KiB - application/x-xz - Now supports Matlab's 'function function' syntax and res return)
file #42968:  firpm.tar.xz added by robsykes (20KiB - application/x-xz - Updated per Octave API changes)
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
  • -email is unavailable- added by lostbard (Posted a comment)
  • -email is unavailable- added by roland65 (Posted a comment)
  • -email is unavailable- added by laurentg (Posted a comment)
  • -email is unavailable- added by robsykes (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only logged-in users can vote.

     

    Follow 19 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2023-05-17 lostbard StatusReady For Test Done
        Open/ClosedOpen Closed
    2023-05-08 lostbard StatusNone Ready For Test
        Assigned tomtmiller lostbard
    2021-04-05 robsykes Attached File- Added firpm.tar.xz, #51196
    2019-03-07 mtmiller Carbon-CopyRemoved 80942 -
    2018-03-30 mtmiller CategoryNone Forge : new function
        Summary[signal] firpm (new remez) [octave forge] (signal) firpm: new function, replacement for remez
    2018-01-18 robsykes Attached File- Added firpm.tar.xz, #42968
    2015-10-11 robsykes Attached File- Added firpm.tar.xz, #35137
    2015-05-21 robsykes Attached File- Added firpm.tar.xz, #34072
    2015-05-04 robsykes Attached File- Added firpm.tar.xz, #33895
    2015-04-26 robsykes Attached File- Added firpm.tar.xz, #33792
    2015-02-11 robsykes Attached File- Added firpm.tar.xz, #33055
    2015-01-08 robsykes Attached File- Added firpm.tar.xz, #32784
    2014-09-16 robsykes Attached File- Added firpm.tar.xz, #32118
    2014-09-16 robsykes Attached File- Added firpm.tar.xz, #32117
    2014-09-15 mtmiller Assigned toNone mtmiller
    2014-09-14 robsykes Attached File- Added firpm.cc.xz, #32091

    Back to the top

    Powered by Savane 3.13-3230.
    Corresponding source code