Home > voicebox > v_windinfo.m

v_windinfo

PURPOSE ^

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

SYNOPSIS ^

function x=windinfo(w,fs)

DESCRIPTION ^

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): 5183, Jan. 1978.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function x=windinfo(w,fs)
0002 %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): 5183, 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 
0160     subplot(212);
0161     nf=min(max(floor(2*max(x.band6,x.band0)*of*nw/fs)+1,of*8),length(p));
0162     ff=(0:nf-1)*fs/(of*nw);
0163     fqi=[x.enbw x.band3 x.band6]/2;
0164     if ff(end)>2000
0165         ff=ff/1000;
0166         fqi=fqi/1000;
0167         xlab='kHz';
0168     else
0169         xlab='Hz';
0170     end
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,max(dd)-dbrange,'facecolor',[1 0.7 0.7]);
0182     hold on
0183     plot(ffs,dbs,':k',ff3,db3,':k',ff6,db6,':k',ffb,dbb,'r',ff,dd,'b');
0184     legend(['Equiv Noise BW = ' sprintsi(x.enbw,-2) 'Hz'],['Max sidelobe = ' sprintf('%.0f',x.sidelobe) ' dB'],['-3 & -6dB BW = ' sprintsi(x.band3,-2) 'Hz & ' sprintsi(x.band6,-2) 'Hz']);
0185     hold off
0186     axis([0 ff(end) max(dd)-dbrange max(dd)+2]);
0187     ylabel('Gain (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);
0195     ylabel('Window');
0196     xlabel('Time (s)');
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=[tcola ','];
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

Generated on Fri 22-Sep-2017 19:37:38 by m2html © 2003