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

[je suis]

> Yes, any public facing documentation, help, etc is fine. Output statements
> produced for whatever input statements you want and however many you want
> is also fine to figure out basic functionality, compatibility issues, and
> to tease out odd and corner cases. Just the code they use to get it is out
> if bounds.

Then, since now I have some time for it, please let me abuse this right a bit...

I suspect there may be precision issues related to numerical solvers.
As I said, I use wxMaxima because I am more familiar with it (than
Octave), and then I convert it to Octave, and for the same input, I
see different results in both. For example, the Q matrix has
Q[1,1]=-2.57024 in wxMaxima, and -1.12857 in Octave. The end results,
for lsfir2(11, [0 200 300 512]/512, [0 1 0 0], 'd'), are:

0.03042732113317503" "
-0.03037155658757951" "
-0.07756410762952007" "
0.0726790980511034" "
0.2633671843212612" "
0.1586739460743956" "
-0.1586739460743956" "
-0.2633671843212612" "
-0.0726790980511034" "
0.07756410762952007" "
0.03037155658757951" "
-0.03042732113317503

in wxMaxima, and in Octave:

   0.0360697335542359
  -0.0399248319344572
  -0.0755659462408800
   0.0787714733962686
   0.2597069270567174
   0.1538852220932733
  -0.1538852220932733
  -0.2597069270567174
  -0.0787714733962686
   0.0755659462408800
   0.0399248319344572
  -0.0360697335542359

But the formulas are the same. I know Maxima tries to keep ratios
whenever possible, so that may be one thing, but I explicitly wrote
float(...) in the for loops, to force a numerical evaluation. Also,
instead of a=Q\b, I have a:invert(Q).b;, along with whatever
implementation is under the hood. There is no fpprec set, no bfloats.

However it is, I remade the script (to also count that there is 1/f^2
weighting in the passband, *only*, as per Matlab's explanation). There
are missing ';' at the end of the lines in several places, to reveal
the inner workings. So, could it be possible to run this script in
Matlab, too? I am curious about the results.

There are also some quirks: there is expint(0) which gives Inf+0i (as
it should), and cos(0)/0, which also gives Inf, but to avoid
singularities, some checks need to be made. Mine, after first using a
check at each step, I simply forced the frequency vector to have the
first value 1e-6, but that means cos(1e-6)/1e-6. A more sensible
choice would be 0 (zero), to avoid numerical instabilities, but then
zero is not 1e6... I hope you see my problem.

Well, I attached the new script, please let me know of the results. If
they are the same in both Octave and Matlab then, well, they must use
a different way to do it.

Vlad

lsfir2.m
Description: application/vnd.wolfram.mathematica.package