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*."
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;
While probably not beautifully efficient, it works. I.e. appears to be
correct.
Regards,
James.
------------------------------------------------------------------------------
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
