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