Re: [OctDev] signal package: possible bug on fir2 with patch (please review)

[Mike M. <mtm...@ie...>]

On Tue, Aug 14, 2012 at 5:33 PM, Carnë Draug <car...@gm...> wrote:
> On 14 August 2012 22:19, Mike Miller <mtm...@ie...> wrote:
>> On Tue, Aug 14, 2012 at 4:54 PM, Carnë Draug wrote:
>>> On 14 August 2012 21:33, Mike Miller <mtm...@ie...> wrote:
>>>> On Tue, Aug 14, 2012 at 3:59 PM, Carnë Draug wrote:
>>>>> On 14 August 2012 20:02, Mike Miller <mtm...@ie...> wrote:
>>>>>> On Tue, Aug 14, 2012 at 1:53 PM, Carnë Draug wrote:
>>>>>>> On 9 August 2012 03:11, Mike Miller <mtm...@ie...> wrote:
>>>>>>>>
>>>>>>>> Finally got back to this, it's now closer to what you showed here.
>>>>>>>> It's clear that Matlab is forcing the 4th argument to always be a
>>>>>>>> power of 2, even though that's not in the help text. I implemented
>>>>>>>> that, added an error message for the case of the argument being too
>>>>>>>> small (after conversion to a power of 2), and added a test case based
>>>>>>>> on your examples above. Thanks for those!
>>>>>>>
>>>>>>> Hi Mike
>>>>>>>
>>>>>>> I have just noticed that this introduced another bug (which I think
>>>>>>> was already there before, this just made it more apparent. For
>>>>>>> example, if the minimum value for 4th argument is more than the
>>>>>>> default, it also automatically (and silently) increases though I'm not
>>>>>>> sure for how much. For example:
>>>>>>
>>>>>> Thanks for checking, I'll take a look at this again today when I get
>>>>>> back to my dev box.
>>>>>>
>>>>>>>>> f = [0 0.1 0.1 0.2 0.2 1]; m = [0 0 1 1 0 0];
>>>>>>>>> bdef = fir2(1024, f, m);
>>>>>>>
>>>>>>> should make the 4th argument default to 512 which would be too low and
>>>>>>> return an error (current behaviour with Octave). However, this does
>>>>>>> not happen in Matlab. But if the default value is specified, then it
>>>>>>> gives an error. And it doesn't default to the next possible value
>>>>>>> (1024) as I first expected, it defaults to 2048.
>>>>>>>
>>>>>>> [...]
>>>>>>>
>>>>>>> Let me know if you need more tests.
>>>>>>
>>>>>> Yes, can you repeat your tests using 1023 and 2047 instead? The limit
>>>>>> is calculated based on the length of the filter which is order+1, I
>>>>>> think that's why you're seeing it go up to the next power of two.
>>>>>>
>>>>>> bdef = fir2(1023, f, m);
>>>>>> b512 = fir2(1023, f, m, 512); % should not error now
>>>>>
>>>>> Yes, it doesn't give an error.
>>>>>
>>>>>> b1024 = fir2(1023, f, m, 1024);
>>>>>> isequal(bdef, b512) % should be true
>>>>>
>>>>> No, it's false
>>>>>
>>>>>> bdef = fir2(2047, f, m);
>>>>>> b1024 = fir2(2047, f, m, 1024);
>>>>>> isequal(bdef, b1024) % should be true
>>>>>
>>>>> No. Also false.
>>>>
>>>> Hmm, ok, so it allows it to be as low as 1/2 but defaults to a larger
>>>> value, how about these, let's try and find the cutoff:
>>>
>>> isequal (fir2 (511, f, m), fir2 (511, f, m, 512))   # true
>>> isequal (fir2 (512, f, m), fir2 (512, f, m, 512))   # true
>>> isequal (fir2 (513, f, m), fir2 (513, f, m, 512))   # true
>>> isequal (fir2 (512, f, m), fir2 (512, f, m, 1024))  # false
>>> isequal (fir2 (513, f, m), fir2 (513, f, m, 1024))   #false
>>> isequal (fir2 (1022, f, m), fir2 (1022, f, m, 512))  #true
>>> isequal (fir2 (1022, f, m), fir2 (1022, f, m, 1024))  # false
>>> isequal (fir2 (1023, f, m), fir2 (1023, f, m, 1024))  #true
>>> isequal (fir2 (2046, f, m), fir2 (2046, f, m, 2048))  #true
>>> isequal (fir2 (2046, f, m), fir2 (2046, f, m, 2048))  #true
>>> isequal (fir2 (2047, f, m), fir2 (2047, f, m, 2048))  #true
>>
>> Great, I think I have enough to fix that one now.
>>
>>>>>> If you have patience for any more, I can come up with some tests to
>>>>>> verify the 5th argument too.
>>>>>
>>>>> Sure, no problem. Give me whatever tests you need and I'll run them.
>>>>> Same for other functions of the package or if you want to take a
>>>>> closer look at differences between results.
>>>>
>>>> Based on how our fir2 is working now, the 5th arg defaults to the 4th
>>>> arg / 20 (and ML's help says the defaults are 512 and 25), so:
>>>> isequal (fir2 (63, f, m, 512), fir2 (63, f, m, 512, 25))
>>>> isequal (fir2 (63, f, m, 1024), fir2 (63, f, m, 1024, 51))
>>>> isequal (fir2 (63, f, m, 2048), fir2 (63, f, m, 2048, 102))
>>>>
>>>> I don't like my chances that these are right out of the box :)
>>>
>>> You are right to not make any bets. All 3 are false in Matlab 2011b.
>>
>> Ok, the docs are wrong yet again, brute force it:
>>
>> for i = [512, 1024, 2048, 4096]
>>   for j = 1:fix(i/10)
>>     if isequal (fir2 (63, f, m, i), fir2 (63, f, m, i, j))
>>       sprintf('default for %d is %d\n', i, j)
>>     end
>>   end
>> end
>
> There you go:
>
> default for 512 is 20
> default for 1024 is 40
> default for 2048 is 81
> default for 4096 is 163
>
> From this, it seems to me that the default of j is "floor (i/25)".

Agreed.

> Where on documentation did you saw that it was /20?

Only piecing together these bits:
* Matlab help [1] for fir2 says "By default, lap is 25."
* Our fir2 had ramp_n = grid_n/20 since before I looked at it

So I'll fix that and fix the help while I'm at it. Thanks for all the test runs.

[1] http://www.mathworks.com/help/toolbox/signal/ref/fir2.html

-- 
mike