v_windinfo

PURPOSE ^

V_WINDINFO window information and figures of merit X=(W,FS)

SYNOPSIS ^

function x=windinfo(w,fs)

DESCRIPTION ^

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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated by m2html © 2003