On 05/07/2012 09:29 PM, James Steward wrote:
> On 08/05/12 11:54, Daniel J Sebald wrote:
>> On 05/07/2012 07:39 PM, James Steward wrote:
>>> On 08/05/12 09:39, James Steward wrote:
>>>> Hi,
>>>>
>>>> Regarding the fir1 function and the "blackmanharris" window;
>>>>
>>>> If I specify an even number of coefficients (n) it produces a symmetric
>>>> filter with an odd number of coefficients (n+1).
>>>>
>>>> If I specify an odd number of coefficients it produces a non symmetric
>>>> filter with an even number of coefficients.
>>
>> N is the order, not the number of coefficients. E.g., second order
>> equation would have three coefficients.
>
> That is evident. It is the symmetry of the filter that is in question,
> which seems to be affected by the number of coefficients being odd or even.
>
>>>> octave:39> version
>>>> ans = 3.0.5
>>>> octave:40> b=fir1(8, 0.5, "blackmanharris")
>>>> b =
>>>>
>>>> Columns 1 through 8:
>>>>
>>>> -7.1851e-08 -2.8280e-03 2.6042e-04 2.7158e-01 6.1193e-01
>>>> 2.7158e-01 2.6042e-04 -2.8280e-03
>>>>
>>>> Column 9:
>>>>
>>>> -7.1851e-08
>>>>
>>>> octave:41> b=fir1(7, 0.5, "blackmanharris")
>>>> b =
>>>>
>>>> -3.8206e-04 -5.2516e-04 2.6785e-02 3.8140e-01 6.1393e-01
>>>> 1.2791e-01 -1.5875e-02 -3.8206e-04
>>>>
>>>> octave:42>
>>>>
>>>> Is there some trick to producing a symmetric "blackmanharris" windowed
>>>> filter with an even number of coefficients?
>>>>
>>>
>>> In reply to myself, I looked at the wikipedia page [1] that describes
>>> the Blackman-Harris filter, and wrote my own blackmanharris.m function
>>> from that. It seems to work just fine, though with no parameter
>>> sensibility checking.
>>>
>>> Then out of curiosity I checked my installed signal library function [2]
>>> and found that it is certainly different.
>>>
>>> [1]
>>> http://en.wikipedia.org/wiki/Window_function#Blackman.E2.80.93Harris_window
>>>
>>> [2] /usr/share/octave/packages/signal-1.0.10/blackmanharris.m
>>
>> Good catch James. The first thing I recognize is that the formula
>> doesn't appear to match what is listed in the Wikipedia reference. In
>> the file [2] is the code
>>
>> a0 = 0.35875;
>> a1 = 0.48829;
>> a2 = 0.14128;
>> a3 = 0.01168;
>> n = -ceil(N/2):N/2;
>> w = a0 + a1.*cos(2.*pi.*n./N) + a2.*cos(4.*pi.*n./N) +
>> a3.*cos(6.*pi.*n./N);
>>
>> But according to [1], it should be "- a1" and "- a3". HOWEVER, I made
>> that change and ended up with a window where the tail is weighted more
>> than the center, so I think that Wikipedia might have the wrong
>> formula. Agreed?
>
> It worked for me, but the definition of n I used was also different.
>
>> Also, this formula inside of blackmanharris.m:
>>
>> n = -ceil(N/2):N/2;
>>
>> does not produce the desired result for N odd.
>
> Correct.
>
>>
>> octave:33> N=7
>> N = 7
>> octave:34> n = -ceil(N/2):N/2
>> n =
>>
>> -4 -3 -2 -1 0 1 2 3
>>
>> Instead of defining 'n' the way the code does, it might be better to
>> use the description from the Wiki page that says "w(n) = w_0(n -
>> (N-1)/2)", this N not being the same as the N in the lines above. Note
>> that blackmanharris.m uses an input L which is the length of window
>> and N = L - 1.
>
> The definition of "n" in the wiki page says "n is an integer, with
> values *0 ≤ n ≤ N-1*."

The N in the Wiki is the length, equivalent to L in blackmanharris.m. 
So the N in blackmanharras.m is the same as N-1 in the Wiki.  Confusing, 
yes.


>>
>> I think the definition of n inside blackmanharris.m should be:
>>
>> n = (0:N) - N/2;
>>
>> Does that produce the result you are expecting James?
>
> I haven't tested this inside the signal toolbox blackmanharris.m file,
> but it yields non integer values of "n" for L = even numbers, N being
> then odd, and N/2 resulting in a non integer value.
>
> My version of blackmanharris.m, that I simply overloaded the system
> default with, looks like this;
>
> function w = blackmanharris (N)
>
> w = 1;
>
> a0 = 0.35875;
> a1 = 0.48829;
> a2 = 0.14128;
> a3 = 0.01168;
>
> for n=0:N-1
> w(n+1) = a0 - a1 * cos(2*pi()*n/(N-1)) + a2 * cos(4*pi()*n/(N-1)) - a3 *
> cos(6*pi()*n/(N-1));
> endfor
> endfunction;

OK, so you are using N as the length.  That matches with the Wiki then. 
  However, loops are discouraged.  Did you use a loop to avoid a 
conditional test that N be greater than zero?  Otherwise, I'd suggest 
the following:

function w = blackmanharris (N)

     if (N < 1)
         error('Window length must be 1 or greater');
     elseif (N == 1)
         w = 1;
     else
         a0 = 0.35875;
         a1 = 0.48829;
         a2 = 0.14128;
         a3 = 0.01168;
         n=0:N-1;
         w = a0 - a1 * cos(2*pi*n./(N-1)) + a2 * cos(4*pi*n./(N-1)) - a3 
* cos(6*pi*n./(N-1));
     endif

endfunction;

The result looks similar to what I was getting from the change I 
suggested in the previous post.  Note one peculiar artifact, which is 
that changing the phasing of the inputs (theoretically correct) results 
in a small round off error as reported by Octave:

octave:181>     for n=0:N-1
  * cos(6*pi()*n/(N-1)); a1 * cos(2*pi()*n/(N-1)) + a2 * 
cos(4*pi()*n/(N-1)) - a3
>     endfor
octave:182> w
w =

  Columns 1 through 6:

    6.0000e-05   5.5645e-02   5.2057e-01   1.0000e+00   5.2058e-01 
5.5645e-02

  Column 7:

    6.0000e-05

The above is indicating 5.2057e-01 and 5.2058e-01 for coefficients, but 
when I list them in long format, why the rounding occurs isn't evident. 
  The rounding effect will likely show up somewhere no matter the formula.

As to the question Mike raised, I'll follow up with another email.

Dan


------------------------------------------------------------------------------
Live Security Virtual Conference
Exclusive live event will cover all the ways today's security and 
threat landscape has changed and how IT managers can respond. Discussions 
will include endpoint security, mobile security and the latest in malware 
threats. http://www.accelacomm.com/jaw/sfrnl04242012/114/50122263/
_______________________________________________
Octave-dev mailing list
Octave-dev@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/octave-dev

Reply via email to