winspec.m

很多matlab的源代码

function [st,w,wh] = winspec(ty,beta)
% WINSPEC Spectral measures of window functions
%
%         [ST,W,WH] = WINSPEC(TY,BETA) computes spectral measures of 
%         window functions specified by TY where TY is a string which can be
%         'rect', 'boxcar', 'vonhann', 'hamming', 'cosine','lanczos','harris'
%	  'bartlett', 'par1','par2','papoulis', 'blackman','kaiser','bickmore'
%	  BETA is an optional parameter for Kaiser window [default: 1]
%	  For harris windows, BETA can be a number from 0 to 7
%	  For help on windows, type HELP WINDOW
%
%         ST is a 6 element vector which returns:
%	  [Peak; PeakSL/Peak; pSL/p(dB);lobewidth; width to PeakSL; w6dB; w3dB]
%         W is the magnitude for F = 0:0.01:10 (1001 points).
%	  WH is the high frequency decay (dB) for FH = 50:0.2:100.
%         Use PLOT(F,W) and PLOT(FH,WH) for the spectra.
%
%         WINSPEC (with no input arguments) invokes the following example:
%	  %Spectral measures of the vonHann window
%         >>whann = winspec('vonhann')


% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: [e-mail]
% Copyright (c) 1998


if nargin==0,help winspec,disp('Strike a key to see results of the example')
pause,whann=winspec('vonhann'),return,end

win=ty;ty=ty(1:4);f=0:.01:10;n=1001;fh=50:.2:100;l=length(fh);
if ty=='rect' | ty=='boxc'
w=(sinc(f));w0=w(1);w=w/w0;n1=102;n2=100;
wh=abs(sinc(fh))/w0;
elseif ty=='hann'|ty=='vonh'
w=(.5*sinc(f)+.25*sinc(f-1)+.25*sinc(f+1));
w0=w(1);w=w/w0;n1=202;n2=200;
wh=abs(.5*sinc(fh)+.25*sinc(fh-1)+.25*sinc(fh+1))/w0;
elseif ty=='hamm'
w=(.54*sinc(f)+.23*sinc(f-1)+.23*sinc(f+1));
w0=w(1);w=w/w0;n1=202;n2=200;
wh=abs(.54*sinc(fh)+.23*sinc(fh-1)+.23*sinc(fh+1))/w0;
elseif ty=='cosi'
ty='cosn';w=(.5*sinc(f+.5)+.5*sinc(f-.5)); %cosine tip
w0=w(1);w=w/w0;n1=152;n2=150;
wh=abs(.5*sinc(fh+.5)+.5*sinc(fh-.5))/w0;
elseif ty=='lanc'
w=.5*(si(pi*f+pi)-si(pi*f-pi))/pi; %lanczos
w0=w(1);w=w/w0;n1=166;n2=164;
wh=.5*abs(si(pi*fh+pi)-si(pi*fh-pi))/pi/w0;
elseif ty=='bart'
w=sinc(f/2);w=.5*w.*w;
w0=w(1);w=w/w0;n1=202;n2=200;
wh=sinc(fh/2);wh=.5*wh.*wh/w0;
elseif ty=='par1'
w=sinc(f/4);w=3*w.*w.*w.*w/8;
w0=w(1);w=w/w0;n1=402;n2=400;
wh=sinc(fh/4);wh=3*wh.*wh.*wh.*wh/8/w0;
elseif ty=='papo'
w=sinc(.5*(f-1))+sinc(.5*(f+1));w=.25*w.*w;
w0=w(1);w=w/w0;n1=302;n2=300;
wh=sinc(.5*(fh-1))+sinc(.5*(fh+1));wh=.25*wh.*wh/w0;
elseif ty=='par2'
f=0.01:.01:10;w=2*(sinc(f)-cos(pi*f))./(f.*f+eps)/pi/pi;w=[2/3 w];
f=[0 f];
w0=w(1);w=w/w0;n1=145;n2=143;
wh=2*abs(sinc(fh)-cos(pi*fh))./(fh.*fh+eps)/pi/pi/w0;
elseif ty=='blac'
w=.42*sinc(f)+.25*sinc(f-1)+.25*sinc(f+1);
w=(w+.04*sinc(f-2)+.04*sinc(f+2));
w0=w(1);w=w/w0;n1=302;n2=300;
wh=.42*sinc(fh)+.25*sinc(fh-1)+.25*sinc(fh+1);
wh=abs(wh+.04*sinc(fh-2)+.04*sinc(fh+2))/w0;
elseif ty=='harr' ,if nargin<2,beta=0;end
if beta==0,n1=302;b=[.375 .5 .125 0];end
if beta==1,n1=402;b=[10/32 15/32 6/32 1/32];end
if beta==2,n1=302;b=[.40897 .5 .09103 0];end
if beta==3,n1=302;b=[.42323 .49755 .07922 0];end
if beta==4,n1=302;b=[.4243801 .4973406 .0782793 0];end
if beta==5,n1=402;b=[.338946 .481973 .161054 .018027];end
if beta==6,n1=402;b=[.355768 .487396 .144232 .012604];end
if beta==7,n1=402;b=[.3635819 .4891775 .1365995 .0106411];end
n2=n1-2;w = b(1)*sinc(f)+0.5*b(2)*(sinc(f-1)+sinc(f+1));
w=(w+0.5*b(3)*(sinc(f-2)+sinc(f+2))+0.5*b(4)*(sinc(f-3)+sinc(f+3)));
w0=w(1);w=w/w0;
wh=b(1)*sinc(fh)+0.5*b(2)*(sinc(fh-1)+sinc(fh+1));
wh=(wh+0.5*b(3)*(sinc(fh-2)+sinc(fh+2))+0.5*b(4)*(sinc(fh-3)+sinc(fh+3)));
wh=abs(wh)/w0;
elseif ty=='kais'
if nargin<2,beta=1;end
a=pi*sqrt(f.*f-beta*beta);w=(sinc(a/pi)/besin(0,pi*beta)); %Kaiser, beta=1
w0=w(1);w=w/w0;
k=fix(150*sqrt(1+beta*beta));df=diff(w(1:k));i=find(df>=0);n2=i(1);n1=n2+2;
a=pi*sqrt(fh.*fh-beta*beta);wh=abs(sinc(a/pi)/besin(0,pi*beta))/w0;
elseif ty=='mkai'
if nargin<2,beta=1;end
a=pi*sqrt(f.*f-beta*beta);w=((sinc(a/pi)-sinc(f))/(besin(0,pi*beta)-1));
w0=w(1);w=w/w0;
k=fix(150*sqrt(1+beta*beta));df=diff(w(1:k));i=find(df>=0);n2=i(1);n1=n2+2;
a=pi*sqrt(fh.*fh-beta*beta);
wh=abs((sinc(a/pi)-sinc(fh))/(besin(0,pi*beta)-1))/w0;
else
error('Unrecognized window type'),return
end
if nargin==1,wt=windo(ty,201);end,if nargin==2,wt=windo(ty,201,beta);end
t=-.5:.005:.5;b=max(abs(w(n1:n)));bdb=20*log10(b);b1=w(1:n2);
b1=20*log10(b1);m=find(b1<=bdb);wid=(m(1)-2)*.01;
b0=-20*log10(sqrt(2));m=find(b1<=b0);b3=(m(1)-2)*.01;
m=find(b1<=2*b0);b6=(m(1)-2)*.01;st=[w0;b;bdb;n2*.01;wid;b6;b3];
%bmax=ceil(1200*b)/1000;
bmax=1.1*b;
i=find(w(n1:n)<0);if ~isempty(i),nb=-bmax;else,nb=0;end
%PLOTTING
disp('STRIKE A KEY BETWEEN PLOTS'), if exist('version')~=5,pause(3),end
 vx=matverch;
 if vx < 4, eval('clg');else,eval('clf');end
plot(t,wt),title([win ' time window']),pause
subplot(121),
axis([0 10 nb w(1)*w0]),plot(f,w*w0),axis([0 10 nb w(1)*w0])
grid,title(['Amplitude Spectrum of ' win ' window'])

subplot(122),axis([0 10 nb bmax]);plot(f,w),axis([0 10 nb bmax]);
grid,title('     Blowup: Normalized'),pause,hold off
w1=abs(w);
i=find(w1==0);if ~isempty(i),[mm,nn]=size(i);w1(i)=1e-7*ones(mm,nn);end
 if vx < 4, eval('clg');else,eval('clf');end
subplot(121),axis([0 10 -120 0]);plot(f,20*log10(w1));axis([0 10 -120 0]);
grid,title('dB magnitude')
wh=20*log10(wh);l1=fix(n2/10);wm=5*ceil(.2*max(wh(1:l1)));
wn=5*floor(.2*max(wh(l-l1:l)))-5;
subplot(122),axis([50 100 wn wm]);plot(fh,wh),axis([50 100 wn wm]);
grid,title('hi freq decay (dB)'),hold off,subplot
if exist('version')~=5,pause,end