Dear all,
Just in case somebody finds this useful, I'm attaching a proposal of a
modified version of the function window() which includes three new
window functions: Blackman, Blackman-Harris and one out of many
different flat-top windows. It also allows a new argument, opt, which
can be either "per" for periodic option or "sym" or any other value (or
no value) for triggering the default symmetric option.
If there were interest, several other windows could be easily added.
Regards,
Federico Miyara
On 11/02/2021 04:12, Federico Miyara wrote:
Dear All,
I wonder why windowing functions such as Hann, Hamming, etc., provided
by window(), are only symmetric.
When used for spectral analysis by subsequent use of fft(), the
periodic weighting is better than the symmetric one. The symmetric
window is mainly used in the design of FIR filters, which I guess is a
less frequent application than spectral analysis.
While it is true that an easy workaround to get a periodic window of
length n is, for instance
w = window("hn", n+1)(1:$-1);
a syntax such as this
w = window("hn", n, "per");
would be easier.Setting "sym" as the default option, no backward
compatibility issues would possibly arise.
Regards,
Federico Miyara
<https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient>
Libre de virus. www.avast.com
<https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient>
<#DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2>
_______________________________________________
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users
--
El software de antivirus Avast ha analizado este correo electrónico en busca de
virus.
https://www.avast.com/antivirus
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA - 1988 - C. Bunks
//
// Copyright (C) 2012 - 2021 - Scilab Enterprises
//
// This file is hereby licensed under the terms of the GNU GPL v2.0,
// pursuant to article 5.3.4 of the CeCILL v.2.1.
// This file was originally licensed under the terms of the CeCILL v2.1,
// and continues to be available under such terms.
// For more information, see the COPYING file which you should have received
// along with this program.
function [win_l,cwp] = window(wtype,n,par,opt)
//[win_l,cwp] = window(wtype,n,par)
//macro that calculates symmetric window
// wtype :window type (re,tr,hn,hm,kr,ch)
// n :window length
// par :parameter 2-vector (kaiser window: par(1)=beta>0)
// : (chebyshev window:par=[dp,df])
// : dp=main lobe width (0<dp<.5)
// : df=side lobe height (df>0)
// opt :window option: "sym" (symmetric)
// "per" (periodic)
// win :window
// cwp :unspecified Chebyshev window parameter
//!
WT = ["re","tr","hm","hn","kr","ch","bl","bh","ft"]
if isdef("opt") then
if opt=="per" then
n = n+1;
end
end
if and(wtype<>WT) then
error(msprintf(_("%s: Wrong value for input argument #%d: Must be in
the set {%s}.\n"),"window",1,strcat(WT,",")))
end
if type(n)<>1|size(n,"*")<>1|~isreal(n)|size(n,"*")<>1|int(n)<>n|n<1 then
error(msprintf(_("%s: Wrong type for input argument #%d: A positive
integer expected.\n"),"window",2))
end
if or(wtype==["kr","ch"]) then
if type(par)<>1|~isreal(par) then
error(msprintf(_("%s: Wrong type for input argument #%d: A %d-by-%d
real vector expected.\n"),"window",3,1,2))
end
if wtype=="kr" then
if size(par,"*")==0| size(par,"*")>2 then
error(msprintf(_("%s: Wrong size for input argument #%d: A
%d-by-%d real vector expected.\n"),"window",3,1,2))
end
Beta = par(1)
if Beta<0 then
error(msprintf(_("%s: Input argument #%d must be strictly
positive.\n"),"window",3))
end
else //chebychev
if size(par,"*")<>2 then
error(msprintf(_("%s: Wrong size for input argument #%d: A
%d-by-%d real vector expected.\n"),"window",3,1,2))
end
dp = par(1);df = par(2)
if dp>0 then
if df>0 then
error(msprintf(_("%s: Wrong value for input argument #%d:
incorrect element #%d\n"),"window",3,2))
end
if dp>=0.5 then
error(msprintf(_("%s: Wrong value for input argument #%d:
incorrect element #%d\n"),"window",3,1))
end
unknown = "df";
else
if df<=0 then
error(msprintf(_("%s: Wrong value for input argument #%d:
incorrect element #%d\n"),"window",3,2))
end
unknown = "dp";
end
end
end
cwp = -1;
//Pre-calculations
no2 = (n-1)/2;
xt = (-no2:no2);
un = ones(1,n);
//Design window according to window type
select wtype
case "re" then //Rectangular window.
win_l = un;
case "tr" then //Triangular window.
win_l = un-2*abs(xt)/(n+1);
case "hm" then //Hamming window.
win_l = 0.54*un+0.46*cos(2*%pi*xt/(n-1));
case "hn" then //Hann window.
win_l = 0.5*un+0.5*cos(2*%pi*xt/(n-1));
case "bl" then //Blackman window
win_l = 0.42*un + 0.5*cos(2*%pi*xt/(n-1)) + 0.08*cos(4*%pi*xt/(n-1));
case "bh" then //Blackman-Harris window
win_l = 0.35875*un - 0.48829*cos(2*%pi*xt/(n-1)) + ...
0.14128*cos(4*%pi*xt/(n-1)) - ...
0.01168*cos(6*%pi*xt/(n-1));
case "ft" then //Flat-top HFT70 window
//https://pure.mpg.de/pubman/faces/ViewItemFullPage.jsp?itemId=item_152164_2
win_l = un - 1.90796*cos(2*%pi*xt/(n-1)) + ...
1.07349*cos(4*%pi*xt/(n-1)) - ...
0.18199*cos(6*%pi*xt/(n-1));
case "kr" then //Kaiser window with parameter beta (n,beta)
//http://en.wikipedia.org/wiki/Kaiser_window
win_l = besseli(0,Beta*sqrt(1-(2*(0:n-1)/(n-1)-1).^2))/besseli(0,Beta);
case "ch" then //Chebyshev window
m = (n-1)/2;
select unknown
case "dp" then
dp = 1/cosh((n-1)*acosh(1/cos(%pi*df)));
cwp = dp;
case "df" then
df = real(acos(1/cosh(acosh((1+dp)/dp)/(n-1)))/%pi)
cwp = df;
end
//Pre-calculation
np1 = int((n+1)/2);
odd = modulo(n,2)==1
x0 = (3-cos(2*%pi*df))/(1+cos(2*%pi*df));
alpha = (x0+1)/2;
Beta = (x0-1)/2;
c2 = (n-1)/2;
//Obtain the frequency values of the Chebyshev window
f = (0:n-1)/n;
xarg = alpha*cos(2*%pi*f)+Beta;
//Evaluate Chebyshev polynomials
// / cos(n acos(x), |x| <= 1
// p(x) = |
// \ cosh(n acosh(x), |x| > 1
pr = zeros(1,n);//real part
pi = zeros(1,n); //imaginary part
ind = find(xarg<=1); pr(ind) = dp*cos (c2*acos (xarg(ind)));
ind = find(xarg>1); pr(ind) = dp*cosh(c2*acosh(xarg(ind)));
pr = real(pr);
if ~odd then
pr = pr.*cos(%pi*f);
pi = -pr.*sin(%pi*f);
antisym = [1*ones(1:int(n/2)+1),-1*ones(int(n/2)+2:n)];
pr = pr.*antisym;
pi = pi.*antisym;
end
//Calculate the window coefficients using the inverse DFT
twn = 2*%pi/n;
xj = (0:n-1);
for i = 1:np1;
arg = twn*(i-1)*xj
w(i) = sum(pr.*cos(arg)+pi.*sin(arg));
end,
c1 = w(1);
w = w/c1;
if odd then
win_l(np1:n) = w(1:np1);
win_l(1:np1-1) = w(np1-1:-1:1);
else,
win_l(np1+1:n) = w(1:np1);
win_l(1:np1) = w(np1:-1:1);
end
win_l = real(win_l');
end
if isdef("opt") then
if opt=="per" then
win_l = win_l(1:$-1);
end
end
endfunction
_______________________________________________
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users
