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

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

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

-- 
mike