firwind.m

很多matlab的源代码

function [hwn,hw,qw] = firwind(ty,a,samp,fp,fs,wind,p1)
% FIRWIND FIR filter design with windows.
%
%       [HW,HWN,QW] = firwind(TY,A,SF,FP,FS,WIND,P1)
%	FIR filter design using windows.
%	TY = 'lp', 'bp', 'bs', 'hp', A = attenuation array [Ap,As],
%       SF=sampling freq., FP=passband, FS=stopband frequencies (Hz),
%	WIND = window name in quotes (DEFAULT: 'hamming')
%	P1 = optional parameter for some windows.Type HELP WINDOW for more
%       HW returns filter coefficients, 
%	HWN returns coefficients of the LP prototype filter 
%       QW returns actual values of n, cutoff freq and attn as QW=[n;fc,Ap,As]
%        NOTE: The length n is chosen as ODD for 'bs' filters OR Dolph window
%        This routine tries to find the smallest filter order by tweaking fc.
%        TO SEE INTERMEDIATE RESULTS USE SF=[SF tr] where tr>0 [Default: tr=0]
%        Returns control to user to change FC and filter length N, if desired.
%
%        FIRWIND (with no input arguments) invokes the following example:
%
%	 A BP filter with Ap=2dB, As=40dB, sampling freq = 100Hz
%	 passband=[20 30]Hz, stopband =[10 40]Hz, window = 'hamming' 
%         >>hw = firwind('bp',[2 40],100,[20 30],[10 40])


% 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


% USAGE: Highpass:  hw = firwind('hp',[3 40],100,30,15);
%        LP (Kaiser window) hw=firwind('lp',[3 50],100,15,22.5,'kaiser');
%        Bandpass: firwind('bp',[2 40],100,[20 30],[18 33]); minmum n=106!
%        But firwind('bp',[2 40],100,[20 30],[18 33],'kaiser'); minimum n=82!

if nargin==0,help firwind,disp('Strike a key to see results of last example')
disp('NOTE: To make no changes, hit the enter key at each prompt')
pause,hbp=firwind('bp',[2 40],[100],[20 30],[10 40]),return,end

if exist('version')==5,ml=4;else,ml=3;end
ty=ty(1:2);if ty=='bp'|ty=='bs',if length(fp)~=2|length(fs)~=2
error('BP or BS filters require fp and fs as two elements'),return,end,end
if ty=='lp'|ty=='hp',if length(fp)>1|length(fs)>1
error('Too many elements in fp or fs'),return,end,end
e=0;if ty=='lp',if fp>fs,e=1;end
elseif ty=='hp',if fs>fp,e=1;end
elseif ty=='bp',pp=abs(diff(fp));ss=abs(diff(fs));if pp>ss,e=1;end
elseif ty=='bs',pp=abs(diff(fp));ss=abs(diff(fs));if ss>pp,e=1;end
else,error('Unknown type. Use lp hp bp or bs in single quotes'),return,end
if e==1,error('Passband and stopband switched'),return,end

if nargin<6,wind='hamming';end
if nargin==7,win='windo(wind,n,p1)';
else,win='windo(wind,n)';end
wnd=wind(1:4);odd=0;ap=a(1);as=a(2);ass=as;
%if ty=='hp',odd=1;end
if ty=='bs',odd=1;end
if wnd=='dolp' | wnd=='cheb';odd=1;
win='windo(wind,n,p1)';if nargin==6,p1=as-8;end,end
if wnd=='kais',win='windo(wind,n,p1)';end

%b FOR KAISER WINDOW
as0=as;if ty=='bs',as0=as+2;end
dels=10^(-.05*as);delp=(10^(0.05*ap)-1)/(10^(0.05*ap)+1);
if delp<dels,as0=-20*log10(delp);end,if as0>=50,b=.1102*(as0-8.7)/pi;
elseif as0>21 & as0<50,b=(.5842*(as0-21)^.4 + .07886*(as0-21))/pi;else,b=0;end

flg=0;if wnd=='kais',flg=1;if nargin==6,p1=b;end,end

if length(samp)>1,tr=samp(2);samp=samp(1);else,tr=0;end %%%TRACE
[wp,ws,w3]=firlpp(ty,2*pi*fp,2*pi*fs,samp);

%Finding N
n=2*fix(pi*2/(ws-wp)/1.3)-3;n=max(3,n); %%OLD

%LENGTH N (new)
%New Section
if exist('remlpord')>0, %For compatibility with v3.5 
n=fix(eval('remlpord(wp/2/pi,ws/2/pi,delp,dels)'));
else
n=fix(1-(10*log10(delp*dels)+13)/(2.324*(ws-wp)));
end
if rem(n,2)==0,n=n+1;end

wx=(wp+ws)/2;sc0=1.25;wx0=wx/sc0;wc=2;

as=ass;clc,a2=0;a3=2*ap;while a2<as | a3>ap
if odd==0,n=n+1;else,n=n+2;end
n0=fix(n/2);n1=1:n0;home,disp(['Selecting n = ' int2str(n)]);
if rem(n,2)==1,h=sin(wx*n1)./n1/pi;h=[fliplr(h) wx/pi h];
else,h=sin(wx*(n1-.5))./(n1-.5)/pi;h=[fliplr(h) h];end
hwin=eval(win);hw=h.*hwin;w0=[wp ws 0:ws/100:ws];
fr=polyval(hw,exp(sqrt(-1)*w0));
a0=-20*log10(abs(fr));a0=a0-min(a0);ftr=w0(1:2)/2/pi;w0(1:2)=[];f0=w0/2/pi;
atr=a0(1:2);%a2=1.05*a0(2);a3=.95*a0(1);%%OLD
%a2=1.1143*a0(2);a3=0.9344*a0(1);
a2=a0(2);a3=0.922*a0(1);
if n>50,if flg<1,a3=.98*a0(1);end,end %%%NEW
a0(1:2)=[];l0=length(a0);

if ws>2.9*wp,wx=wx0;else,i=find(a0>ap);if isempty(i),k=1;else,k=i(1);end
if k==1,sc=1;else,sc=wp/w0(k);end,wx=wx*sc;end

%if wx>ws,if flag==1,wx=0.5*(wp+ws);else,wx=wx0;end,end
if wx>ws,wx=wx0;end

if rem(n,50)==0 & n>60,c=input(['n>' int2str(n) '. Enter 1 to continue!']);clc
if c~=1,return,end,end
end
vx=matverch;
if vx < 4,eval('clg');else,eval('clf');end
g='Digital Freq F';db='Magnitude dB';


if ml==4,er='erasemode';xo='xor';mk='m*';xd='xdata';yd='ydata';xl='XLim';
eval('hold on,grid on,axis([0 0.5 -2*as 0])')
plt1=eval('plot(ftr,atr,mk,er,xo)');plt2=eval('plot(f0,a0,er,xo)');
if tr>0,xlabel(g);ylabel(db);end,end


wxx=.5*(wp+ws);if ty=='bs',add0=.02;else,add0=-.02;wxx=1.2*wxx;end
wx=wxx;add=add0;
%%%%%%%%%%%%% NEWNEWNEW %%%%%%%%%%%%%%%%%%%
% while ~isempty(n),
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n0=fix(n/2);n1=1:n0;hwin=eval(win);
num=int2str(n);acp=2*ap;ac0=as;ap0=acp;

%FOR ML3: COMMENT THE FOLLOWING LINE IF USING NEWPLOT
if tr==0 | ml==4 %%TRACE REPLACE LINE BELOW
%if ml==4,
home,disp([' Selected filter length n = ' int2str(n)]);end

%%%NEW FOR PLOTTING
if tr>0 %%TRACE
if ml==3
y=-20*round(abs(as)/10);hold off,axis([0 .5 y 0]);
ff0=[0 0];aa0=[0 0];fftr=ff0;aatr=aa0;plot(ff0,aa0,'.w',fftr,aatr,'.w')
grid,xlabel(g),ylabel('Magnitude dB'),
title(['Filter Length n = ' int2str(n) ' Tweaking LPP Cutoff Frequency']);
hold on
end
end %%TRACE
%%%END OF NEW FOR PLOTTING

acomp=0;flag=0;while acomp<as
if flag==0,if as-acomp<10,add=add/10;flag=1;end,end
wx=wx+add;

%COMMENT THE FOLLOWING LINE OUT IF USING NEWPLOT
if tr==0 | ml==4 %%TRACE REPLACE LINE BELOW
%if ml==4,
home,disp(' '),disp([' Tweaking LPP cutoff freq Fc = ' num2str(wx/2/pi)]),end

if acomp-as>10,disp('Cannot meet requirements'),wx=wxx;break,end
if wx<0.8*wp;disp('Cannot meet requirements'),wx=wxx;break,end
ac0=acomp;ap0=acp;if rem(n,2)==1,h=sin(wx*n1)./n1/pi;h=[fliplr(h) wx/pi h];
else,h=sin(wx*(n1-.5))./(n1-.5)/pi;h=[fliplr(h) h];end,hw=h.*hwin;
n0=fix(n/2);if ty=='lp',hwn=hw;end
if rem(n,2)==1,nn=-n0:n0;else,nn=[-n0+.5:-.5 .5:n0-.5];end
if ty=='hp',hwn=hw.*((-1).^nn);end,if ty=='bp',hwn=2*cos(nn*w3).*hw;end
if ty=='bs',hwn=-2*cos(nn*w3).*hw;hwn(n0+1)=1+hwn(n0+1);end
hwn=hwn.*(abs(hwn)>eps);w1=2*pi*fp/samp;w2=2*pi*fs/samp;l=length([w1 w2]);
w0=[w1 w2 0:pi/100:pi];
fr=polyval(hwn,exp(sqrt(-1)*w0));mag=abs(fr);a0=20*log10(mag);
fr(1:l)=[];ftr=w0(1:l)/2/pi;w0(1:l)=[];atr=a0(1:l);a0(1:l)=[];ma=mag(1:l);
mag(1:l)=[];am=max(a0);f0=w0/2/pi;a0=a0-am;atr=atr-am;l0=length(a0);

%%%NEW FOR PLOTTING
cut=num2str(wx/2/pi);num=int2str(n);
if length(atr)<3,ats=num2str(-min(atr));acp=abs(max(atr));as2='';
else,
ats=num2str(-atr(3));as2=num2str(-atr(4));acp=max(abs(atr(1:2)));
end
atp=num2str(acp);

if tr>0 %%TRACE
if ml==3
plot(ff0,aa0,'-i',fftr,aatr,'*i'),
plot(f0,a0,ftr,atr,'*'),grid
end
end %%TRACE
ff0=f0;aa0=a0;fftr=ftr;aatr=atr;
%%%END OF NEW FOR PLOTTING

if length(atr)<3,acomp=abs(min(atr));acp=abs(max(atr));
else,
acomp=min(abs(atr(3:4)));acp=max(abs(atr(1:2)));
end
if tr>0 %%TRACE
if ml==4,eval('set(plt1,xd,ftr,yd,atr),set(plt2,xd,f0,yd,a0),drawnow'),
%title([wind ' n=' num ' fc=' cut ' ap=' atp ' as=(' ats ' ' as2 ')']),
end
end %%TRACE
end
if tr>0 %%TRACE
if ml==3,hold off,end
end %%TRACE

if tr==0,
if ml==3
plot(ff0,aa0,'-i',fftr,aatr,'*i'),
plot(f0,a0,ftr,atr,'*'),grid
end
if ml==4,eval('set(plt1,xd,ftr,yd,atr),set(plt2,xd,f0,yd,a0),drawnow'),end
end
title([wind ' n=' num ' fc=' cut ' ap=' atp ' as=(' ats ' ' as2 ')']),


disp(' '),disp('AUTO DESIGN NOW COMPLETE')
if acp>ap,disp('WARNING: Filter length n may need to be increased'),else
disp('You MAY be able to DECREASE the filter length n'),end
disp('You may now modify n, Fc (try steps of 0.001) etc')


par = [n,wx/2/pi];dis='N,Fc';if odd==1,dis='N(odd),Fc';end
if wnd=='kais',par=[par,p1];dis=[dis,',beta'];end
%par1=par;
par2=par;
dis=['[' dis ']'];disp(' ');
%disp(['Ap=' atp '  As=(' ats ' ' as2 ')'])

disp(['AUTO DESIGN VALUES for ' dis]);

pars=['n=' int2str(n) '   Fc=' sprintf('%.5g',wx/2/pi)];
if wnd=='kais',pars=[pars '   beta=' sprintf('%.5g',p1)];end

flag1=0;par1=-1;
while ~isempty(par1)
% disp(par)
if par1~=-1
par1=input([ dis ' or 0 for AUTO or press ENTER to stop :']);
if isempty(par1),par1=par2;flag1=1;end
else
par1=0;
end

%if par1==0,format long;disp(par);format;par1=par;end
if par1==0,par1=par;end
n=par1(1);if odd==1,if rem(n,2)==0,
n=n+1;par1(1)=n;disp('WARNING: INCREASING n by 1 for odd length'),end,end

if length(par1)>1,wx=par1(2)*2*pi;end
if length(par1)>2,p1=par1(3);end


%%%%%%%%%NEWNEWNEWNEW %%%%%%%%%%%%%%%%%%%%%%%%%%
%while wx~=0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%% ADDED NEW %%%%%%%%%%%%%%%%%%%%%%%%%%
n0=fix(n/2);n1=1:n0;hwin=eval(win);
num=int2str(n);acp=2*ap;ac0=as;ap0=acp;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if rem(n,2)==1,h=sin(wx*n1)./n1/pi;h=[fliplr(h) wx/pi h];
else,h=sin(wx*(n1-.5))./(n1-.5)/pi;h=[fliplr(h) h];end,hw=h.*hwin;
n0=fix(n/2);if ty=='lp',hwn=hw;end
if rem(n,2)==1,nn=-n0:n0;else,nn=[-n0+.5:-.5 .5:n0-.5];end
if ty=='hp',hwn=hw.*((-1).^nn);end,if ty=='bp',hwn=2*cos(nn*w3).*hw;end
if ty=='bs',hwn=-2*cos(nn*w3).*hw;hwn(n0+1)=1+hwn(n0+1);end
hwn=hwn.*(abs(hwn)>eps);w1=2*pi*fp/samp;w2=2*pi*fs/samp;l=length([w1 w2]);
w0=[w1 w2 0:pi/400:pi];
fr=polyval(hwn,exp(sqrt(-1)*w0));mag=abs(fr);a0=20*log10(mag);
fr(1:l)=[];ftr=w0(1:l)/2/pi;w0(1:l)=[];atr=a0(1:l);a0(1:l)=[];ma=mag(1:l);
mag(1:l)=[];am=max(a0);f0=w0/2/pi;a0=a0-am;atr=atr-am;y=-20*round(as/10);
g='Digital Freq F';axis([0 .5 y 0]);cut=num2str(wx/2/pi);num=int2str(n);

if length(atr)<3,ats=num2str(-min(atr));acp=abs(max(atr));as2='';
else,
ats=num2str(-atr(3));as2=num2str(-atr(4));acp=max(abs(atr(1:2)));
end
atp=num2str(acp);
if ml==4,eval('set(plt1,xd,ftr,yd,atr),set(plt2,xd,f0,yd,a0),drawnow')
else,plot(f0,a0,ftr,atr,'*'),axis([0 .5 y 0]);grid,
xlabel(g),ylabel('Magnitude dB'),end

title([wind ' n=' num ' fc=' cut ' ap=' atp ' as=(' ats ' ' as2 ')']),

%%%%%%%%% NEW  NEW %%%%%%%%%%%%%%%%%%%%%
if ml==3,pause,end
%pause
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pars=['n=' int2str(n) '   Fc=' sprintf('%.5g',wx/2/pi)];
if wnd=='kais',pars=[pars '   beta=' sprintf('%.5g',p1)];end
disp(['Ap=' atp '  As=(' ats ' ' as2 ')'])
disp(pars)


% wx1=input('New fc or press ENTER to stop :');
% if isempty(wx1),break,else,wx=wx1*2*pi;end

%%%%%%%%%%%%%%NEWNEWNEW %%%%%%%%%%%%%%%%
%%%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% n2=input('enter new n or press ENTER to stop :');if isempty(n2),break,end
% if rem(n2,2)==0,if odd==1,n2=n2+1;end,end
% wx=wxx;add=add0;flag=0;n=n2;clc

%%%%%%%%%%%%%% NEWNEWNEW %%%%%%%%%%%%%%%%%%%%%
%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if flag1==1,break;else,par2=par1;end
end


%FINAL PLOTS
ang=angle(fr);phu=unwrap(ang,.75*pi)*180/pi;ph=ang*180/pi;
subplot(221),mx=max(mag);mag=mag/mx;ma=ma/mx;
axis([0 .5 y 0]),plot(f0,a0,ftr,atr,'o'),axis([0 .5 y 0]),grid
xlabel(g),ylabel('Magnitude dB')
subplot(222)
axis([0 .5 0 1.2]),plot(f0,mag,ftr,ma,'o'),axis([0 .5 0 1.2]),grid
xlabel(g),ylabel('Magnitude'),if ml==3,hold off,end
subplot(223)
plot(f0,phu),if ml==4,eval('set(gca,xl,[0 0.5])'),end
grid,xlabel(g),ylabel('Unwrapped Phase')
subplot(224),
bx=[0 0.5 -200 200];
if ty=='hp',axis(bx);plot(f0,ph),axis(bx);grid,xlabel(g),ylabel('Phase deg.')
elseif ty=='bs'
axis(bx);plot(f0,ph),axis(bx);grid,xlabel(g),ylabel('Phase deg.')
elseif ty=='lp',%ax=[0 w1/2/pi ma(1) 1];
i=find(f0>=w1/2/pi);i=i(1);mg=mag(1:i);y1=min(mg);y2=max([max(mg) 1 ma]);
ax=[0 w1/2/pi y1 y2];axis(ax);
plot(f0,mag,ftr,ma,'o'),axis(ax),grid,xlabel(g),ylabel('Passband Magnitude')
else
ax=[w1/2/pi min(ma(1:2)) 1];axis(ax);plot(f0,mag,ftr,ma,'o'),axis(ax),
grid,xlabel(g),ylabel('Passband Magnitude')
end
subplot
if ml==3,
title([wind ' n=' num ' fc=' cut ' ap=' atp ' as=(' ats ' ' as2 ')'])
hold off,pause,end
qw=[n;wx/2/pi;atr'];