V_WINDINFO window information and figures of merit X=(W,FS)
Inputs: W is a vector containing the window
FS is the sampling frequency (default=1)
Outputs: X.len length of the window (samples)
X.nw length of the window (samples)
X.ewgdelay energy centroid delay from first sample (samples)
X.dcgain DC gain (dB)
X.sidelobe maximum sdelobe level in dB relative to DC gain
X.falloff rate at which sidelobes decay (dB/octave)
X.enbw equivalent noise bandwidth (*fs/len Hz)
X.scallop scalloping loss (dB)
X.ploss processing loss (dB)
X.wcploss worst case processing loss (dB)
X.band3 3dB bandwidth (Hz)
X.band6 6 dB bandwidth (Hz)
X.band0 essential bandwidth (to first minimum) (Hz)
X.gain0 gain at first minimum (Hz)
X.olc50 50% overlap correction
X.olc75 75% overlap correction
X.cola overlap factors giving constant overlap add (exlcuding multiples)
X.cola2 as X.cola but for squared window
If no output argument is given, the window and frequency response
will be plotted e.g. windinfo(windows('hamming',256,'ds'),256);
To obtain the figures of merit listed in Table 1 of [1] set
fs = length(W), multiply X.olc50 and X.olc75 by 100%. The "coherent gain
listed in the table is 10^(x.dcgain/20)/(max(w)*length(w)).
[1] F. J. Harris. On the use of windows for harmonic analysis with the
discrete fourier transform. Proc IEEE, 66 (1): 51�83, Jan. 1978.0001 function x=windinfo(w,fs) 0002 %V_WINDINFO window information and figures of merit X=(W,FS) 0003 % 0004 % Inputs: W is a vector containing the window 0005 % FS is the sampling frequency (default=1) 0006 % 0007 % Outputs: X.len length of the window (samples) 0008 % X.nw length of the window (samples) 0009 % X.ewgdelay energy centroid delay from first sample (samples) 0010 % X.dcgain DC gain (dB) 0011 % X.sidelobe maximum sdelobe level in dB relative to DC gain 0012 % X.falloff rate at which sidelobes decay (dB/octave) 0013 % X.enbw equivalent noise bandwidth (*fs/len Hz) 0014 % X.scallop scalloping loss (dB) 0015 % X.ploss processing loss (dB) 0016 % X.wcploss worst case processing loss (dB) 0017 % X.band3 3dB bandwidth (Hz) 0018 % X.band6 6 dB bandwidth (Hz) 0019 % X.band0 essential bandwidth (to first minimum) (Hz) 0020 % X.gain0 gain at first minimum (Hz) 0021 % X.olc50 50% overlap correction 0022 % X.olc75 75% overlap correction 0023 % X.cola overlap factors giving constant overlap add (exlcuding multiples) 0024 % X.cola2 as X.cola but for squared window 0025 % 0026 % If no output argument is given, the window and frequency response 0027 % will be plotted e.g. windinfo(windows('hamming',256,'ds'),256); 0028 % 0029 % To obtain the figures of merit listed in Table 1 of [1] set 0030 % fs = length(W), multiply X.olc50 and X.olc75 by 100%. The "coherent gain 0031 % listed in the table is 10^(x.dcgain/20)/(max(w)*length(w)). 0032 % 0033 % [1] F. J. Harris. On the use of windows for harmonic analysis with the 0034 % discrete fourier transform. Proc IEEE, 66 (1): 51�83, Jan. 1978. 0035 0036 % Copyright (C) Mike Brookes 2009-2014 0037 % Version: $Id: v_windinfo.m 6801 2015-09-12 09:30:42Z dmb $ 0038 % 0039 % VOICEBOX is a MATLAB toolbox for speech processing. 0040 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0041 % 0042 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0043 % This program is free software; you can redistribute it and/or modify 0044 % it under the terms of the GNU General Public License as published by 0045 % the Free Software Foundation; either version 2 of the License, or 0046 % (at your option) any later version. 0047 % 0048 % This program is distributed in the hope that it will be useful, 0049 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0050 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0051 % GNU General Public License for more details. 0052 % 0053 % You can obtain a copy of the GNU General Public License from 0054 % http://www.gnu.org/copyleft/gpl.html or by writing to 0055 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0056 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0057 % 0058 0059 if nargin<2 0060 fs=1; 0061 end 0062 w=w(:); 0063 nw=length(w); 0064 x.len=nw/fs; 0065 x.nw=nw; 0066 % energy weighted group delay = centre of energy 0067 x.ewgdelay=((1:nw)*w.^2/sum(w.^2)-1)/fs; 0068 % now calculate spectrum 0069 of=16; % spectrum oversample factor must be even 0070 nwo=of*nw; 0071 f=rfft(w,nwo); 0072 p=f.*conj(f); 0073 % sidelobe attenuation is maximum peak (note DC peak at p(1) is not found) 0074 [kp,vp]=v_findpeaks(p,'q'); 0075 [kt,vt]=v_findpeaks(-p,'q'); 0076 if ~numel(kp) 0077 x.sidelobe=10*log10(min(p)/p(1)); 0078 else 0079 x.sidelobe=10*log10(max(vp)/p(1)); 0080 end 0081 np=length(kp); 0082 ipa=floor(np/4); 0083 if ~ipa 0084 x.falloff=0; 0085 else 0086 ipb=floor(np/2); 0087 x.falloff=10*log10(vp(ipb)/vp(ipa))/log2((ipb-1)/(ipa-1)); 0088 end 0089 sumw2=sum(w.^2); 0090 sumw=sum(w); 0091 enbwbin=nw*sumw2/sumw^2; 0092 x.enbw=enbwbin*fs/nw; 0093 x.dcgain=20*log10(sumw); 0094 % do linear interpolation in p() to find 3dB and 6dB points 0095 p3=0.5*p(1); 0096 i3=find(p<p3,1); 0097 if ~numel(i3) 0098 x.band3=Inf; 0099 x.band6=Inf; 0100 else 0101 x.band3=2*(i3-(p3-p(i3))/(p(i3-1)-p(i3))-1)/of*fs/nw; 0102 p6=0.25*p(1); 0103 i6=find(p<p6,1); 0104 x.band6=2*(i6-(p6-p(i6))/(p(i6-1)-p(i6))-1)/of*fs/nw; 0105 end 0106 % do linear interpolation in f() to find closest approach to the origin 0107 if~numel(kt) 0108 x.band0=Inf; 0109 x.gain0=0; 0110 else 0111 i0=floor(kt(1)); 0112 df=f(i0+1)-f(i0); 0113 j0=-real(f(i0)*conj(df))/abs(df)^2; 0114 x.band0=2*(i0+j0-1)/of*fs/nw; 0115 p0=abs(f(i0)+j0*df)^2; 0116 if p0>0 0117 x.gain0=10*log10(p0/p(1)); 0118 else 0119 x.gain0=-Inf; 0120 end 0121 end 0122 % overlap factors 0123 i50=round(nw*0.5); 0124 x.olc50=sum(w(1:nw-i50).*w(1+i50:nw))/sumw2; 0125 i75=round(nw*0.25); 0126 x.olc75=sum(w(1:nw-i75).*w(1+i75:nw))/sumw2; 0127 % processing loss and scalloping loss 0128 x.scallop=10*log10(p(1)/p(1+of/2)); 0129 x.ploss=10*log10(enbwbin); 0130 x.wcploss=x.ploss+x.scallop; 0131 co=zeros(1,nw); 0132 co2=co; 0133 w2=w.^2; 0134 for i=1:nw 0135 if i*fix(nw/i)==nw % check i is a factor of the window length 0136 if co(i)==0 0137 co(i)=all(abs(sum(reshape(w',nw/i,i),2)-i*mean(w))<max(abs(w))*1e-4); 0138 if co(i) 0139 co(2*i:i:nw)=-1; % ignore multiples 0140 end 0141 end 0142 if co2(i)==0 0143 co2(i)=all(abs(sum(reshape(w2',nw/i,i),2)-i*mean(w2))<max(w2)*1e-4); 0144 if co2(i) 0145 co2(2*i:i:nw)=-1; % ignore multiples 0146 end 0147 end 0148 end 0149 end 0150 co(co<0)=0; 0151 co2(co2<0)=0; 0152 x.cola=find(co); 0153 x.cola2=find(co2); 0154 % 0155 % now plot it if no output arguments given 0156 % 0157 if ~nargout 0158 clf; 0159 subplot(212); 0160 nf=min(max(floor(2*max(x.band6,x.band0)*of*nw/fs)+1,of*8),length(p)); 0161 ff=(0:nf-1)*fs/(of*nw); 0162 fqi=[x.enbw x.band3 x.band6]/2; 0163 if ff(end)>2000 0164 ff=ff/1000; 0165 fqi=fqi/1000; 0166 xlab='kcyc/L'; 0167 else 0168 xlab='cyc/L'; 0169 end 0170 dbn=20*log10(x.nw); % window width in dB 0171 dbrange=min(100,-1.5*x.sidelobe); 0172 dd=10*log10(max(p(1:nf),p(1)*0.1^(dbrange/10))); 0173 ffs=[0 ff(end)]; 0174 dbs=repmat(x.dcgain+x.sidelobe,1,2); 0175 ffb=[0 fqi(1) fqi(1)]; 0176 dbb=[dd(1) dd(1) dd(1)-dbrange]; 0177 ff3=[0 fqi(2) fqi(2)]; 0178 db3=[dd(1)+db(0.5)/2 dd(1)+db(0.5)/2 dd(1)-dbrange]; 0179 ff6=[0 fqi(3) fqi(3)]; 0180 db6=[dd(1)+db(0.5) dd(1)+db(0.5) dd(1)-dbrange]; 0181 area(ffb,dbb-dbn,max(dd)-dbrange-dbn,'facecolor',[1 0.7 0.7]); 0182 hold on 0183 plot(ffs,dbs-dbn,':k',ff3,db3-dbn,':k',ff6,db6-dbn,':k',ffb,dbb-dbn,'r',ff,dd-dbn,'b'); 0184 legend(['Equiv Noise BW = ' sprintsi(x.enbw,-2) 'cyc/L'],['Max sidelobe = ' sprintf('%.0f',x.sidelobe) ' dB'],['-3 & -6dB BW = ' sprintf('%.2g',(x.band3)) ' & ' sprintf('%.2g',(x.band6)) ' cyc/L']); 0185 hold off 0186 axis([0 ff(end) max(dd)-dbrange-dbn max(dd)+2-dbn]); 0187 ylabel('Gain/N (dB)'); 0188 xlabel(sprintf('Freq (%s)',xlab)); 0189 % 0190 % Now plot the window itself 0191 % 0192 subplot(211); 0193 tax=(0:nw-1)/fs-x.ewgdelay; 0194 area(tax,w,'FaceColor',[0.7 0.7 1]); 0195 ylabel('Window'); 0196 xlabel('Time/L'); 0197 dtax=(tax(end)-tax(1))*0.02; 0198 axv=[tax(1)-dtax tax(end)+dtax min(0,min(w)) max(w)*1.05]; 0199 texthvc(tax(end),max(w),sprintf('N=%d',nw),'rtk'); 0200 if length(x.cola)>3 0201 tcola=sprintf(',%d',x.cola(1:3)); 0202 tcola=[tcola ',...']; 0203 else 0204 tcola=sprintf(',%d',x.cola); 0205 end 0206 if length(x.cola2)>3 0207 tcola2=sprintf(',%d',x.cola2(1:3)); 0208 tcola2=[tcola2 ',...']; 0209 else 0210 tcola2=sprintf(',%d',x.cola2); 0211 end 0212 texthvc(tax(1),max(w),sprintf('COLA=%s\nCOLA^2=%s',tcola(2:end),tcola2(2:end)),'ltk'); 0213 axis(axv); 0214 end 0215 0216 0217