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

[Nicholas Jankowski]

On Jun 2, 2017 2:16 AM, "je suis" <address@hidden> wrote:

>>> h = firls(11, [0 200 300 512]/512, [0 1 0 0], 'd');h'
> ans =
>    0.025308869942737
>   -0.018174108886369
>   -0.083194157264830
>    0.063294222003005
>    0.270747839846565
>    0.167064335397056
>   -0.167064335397056
>   -0.270747839846565
>   -0.063294222003005
>    0.083194157264830
>    0.018174108886369
>   -0.025308869942737

ans =
   0.0269498757347156
  -0.0215227992037768
  -0.0818317682801729
   0.0657678813748652
   0.2689004493151205
   0.1649294765945745
  -0.1649294765945745
  -0.2689004493151205
  -0.0657678813748652
   0.0818317682801729
   0.0215227992037768
  -0.0269498757347156

>>> h = firls(10, [0 200 300 512]/512, [0 1 0 0], 'd');h'
> ans =
>    0.035928864351330
>   -0.067308010314644
>   -0.050456898575185
>    0.185894253348846
>    0.280660645384262
>                    0
>   -0.280660645384262
>   -0.185894253348846
>    0.050456898575185
>    0.067308010314644
>   -0.035928864351330

ans =
   0.036134459973126
  -0.068051420737790
  -0.049405504477138
   0.186115545919129
   0.279580844271503
   0.000000000000000
  -0.279580844271503
  -0.186115545919129
   0.049405504477138
   0.068051420737790
  -0.036134459973126

And, using another example of yours, from previous replies:

>>> 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

this is what I get:

ans =
  -0.0001933377513559
  -0.0002863125210171
   0.0004931347624372
   0.0008827598349943
  -0.0000178929005084
  -0.0014001289943040
  -0.0013489576174388
   0.0007299431745083
   0.0026985362524900
   0.0016862930823029
  -0.0021152971762283
  -0.0045450417812548
  -0.0016992508119037
   0.0046183535728281
   0.0072669680159855
   0.0011461121669188
  -0.0094070870425227
  -0.0124382329716109
   0.0001614104478731
   0.0223140883374915
   0.0362886499966314
   0.0278092007964514
   0.0000000000000000
  -0.0278092007964514
  -0.0362886499966314
  -0.0223140883374915
  -0.0001614104478731
   0.0124382329716109
   0.0094070870425227
  -0.0011461121669188
  -0.0072669680159855
  -0.0046183535728281
   0.0016992508119037
   0.0045450417812548
   0.0021152971762283
  -0.0016862930823029
  -0.0026985362524900
  -0.0007299431745083
   0.0013489576174388
   0.0014001289943040
   0.0000178929005084
  -0.0008827598349943
  -0.0004931347624372
   0.0002863125210171
   0.0001933377513559

So the numbers seem to come closer the higher the order. I don't
understand this and one of the reasons probably involves the "ghost
tweak" (pi) on lines 251 (q -= pi*...) and 322 (b = pi*...) -- which I
got to by accident, see previous replies -- and the workaround for
division by zero on line 209: change 1e-6 to something else and see
the results change. For example, if you change the pi tweak to 4
(lines 251 and 322), for the first two examples in your last reply
(orders 11 and 10) you'll get much closer results than for the 44th
order. Change it back to pi and the higher order benefits, not the
lower ones. This makes me suspect rounding up issues, rather than the
method itself.

> looking back through emails I don't see any others you requested. Let me
> know if any other input variations might be helpful.

You could see what I meant above by using the uploaded script. Still,
I tried plotting the derivative of the abs(FFT) for both Matlab's
results and mine, and it looks like my numbers actually perform better
closer to DC for the 44th order (and tweak=pi), but worse for lower
orders (and tweak=4).

At this point, I'm baffled, but my conclusion is that there are
rounding up issues (for Matlab, too), probably due to workarounds for
division by zero, but in the lack of a clear, reference method, it may
be inconclusive which method fares better or closer to the truth. I
tried using wxMaxima with bfloat(), but it only uses one core (out of
2...) and it takes forever while my CPU spits fire.

Vlad

Have you read through Matlabs public documentation, including the algorithm section? Any insights there? 

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

Do you gave access to either their primary algorithm reference:

[1] Parks, Thomas W., and C. Sidney Burrus. Digital Filter Design. New York: John Wiley & Sons, 1987, pp. 54–83.

Or this supplemental they mention in the algorithm section about some specific cases?

[2] Oppenheim, Alan V., Ronald W. Schafer, and John R. Buck. Discrete-Time Signal Processing. Upper Saddle River, NJ: Prentice Hall, 1999.

Lastly, if you set up a github or other repository others might be willing to help hack your code. Keeping track of it via email is a bit tedious.