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