# 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:

- v_findpeaks V_FINDPEAKS finds peaks with optional quadratic interpolation [K,V]=(Y,M,W,X)

This function is called by:

- v_windows V_WINDOWS Generate a standard windowing function (TYPE,N,MODE,P,OV)

## 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 Lesser General Public License as published by
    0046 %   the Free Software Foundation; either version 3 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 Lesser General Public License for more details.
    0053 %
    0054 %   You can obtain a copy of the GNU Lesser General Public License from
    0055 %   https://www.gnu.org/licenses/ .
    0056 %
    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