firhb.m

很多matlab的源代码

function [hwn,hw,qw] = firhb(ty,a,samp,fp,fs,wind,p1)
% FIRHB Halfband FIR filter design.
%
%	[HW,HWN,QW] = FIRHB(TY,A,SF,fp,fs,WIN,P1) Halfband FIR filter
%       TY = 'lp', 'bp', 'bs', or 'hp',   SF = sampling frequency in HERTZ. 
%       A = [Ap,As] = attn dB   fp,fs=passbad and stopband edge(s) in HERTZ
%	NOTE: SF must equal 2*(fp+fs). If not, it is set to this value.
%       WIN = window name ( as a string) [Default: WIN = 'kaiser']
%	P1 = Optional parameter for some windows.Type HELP WINDOW for more
%       HW = coefficients of actual filter.   HWN = coeffs of FIR LPP.
%	QW is a 4x1 vector of actual values corresponding to 
%	QW = [filter length; LPP cutoff freq; actual Ap; actual As; P1;P2]
%       NOTE: The passband and stopband are made to have arithmetic symmetry 
%       by keeping the passband (for BP) or stopband (for BS) edge(s) fixed.
%	Program returns control to user to change filter length n if desired.
%       Note: The filter length n must be chosen as ODD. 
%
%       FIRHB (with no input arguments) invokes the following example:
%
%       % Design a Kaiser halfband HP filter with attenuation = [3 40]dB, 
%	% passband/stopband edge = 30Hz, 20Hz. Then SF = 2*(20+30) = 100Hz
%       >>hhb = firhb('hp',[3 40],100,30,20);


% 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 firhb,disp('Strike a key to see results of the example')
disp('NOTE: If no changes are needed, hit the enter key at each prompt')
pause,hhb=firhb('hp',[3 40],100,30,20),return,end

p1=[];
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='kaiser';end
if nargin==7,win='windo(wind,n,p1)';
else,win='windo(wind,n)';end
wnd=wind(1:4);ap=a(1);as=a(2);as0=as;
if wnd=='dolp' | wnd=='cheb';win='windo(wind,n,p1)';
if nargin==6,p1=as-8;end,end
if wnd=='kais',win='windo(wind,n,p1)';end

if ty=='lp'| ty=='hp',samp=2*(fp+fs);df=abs(fs-fp);
else,df1=abs(fp(1)-fs(1));df2=abs(fs(2)-fp(2));df=min(df1,df2);
if ty=='bp',fs(1)=fp(1)-df;fs(2)=fp(2)+df;end
if ty=='bs',fp(1)=fs(1)-df;fp(2)=fs(2)+df;end,samp=2*sum(fp);end
[wp,ws,w3]=firlpp(ty,2*pi*fp,2*pi*fs,samp);

wc=(wp+ws)/2;fcn=num2str(wc/2/pi);
dels=10^(-.05*as);delp=(10^(0.05*ap)-1)/(10^(0.05*ap)+1);
if delp<dels,as=-20*log10(delp);end
if as>50,b=.1102*(as-8.7)/pi;
elseif as>21,b=(.5842*(as-21)^.4 + .07886*(as-21))/pi;else,b=0;end
if as>21,d=(as-7.95)/14.36;else,d=0.9222;end
n=ceil(samp*d/df+1);if rem(n,2)==0,n=n+1;end

if wnd=='kais',if nargin<7,p1=b;end,end

as=as0;c=0;
if nargin<7 & wnd=='kais'
n1=1:fix(n/2);h=[sin(wc*n1)./n1/pi];
h=[fliplr(h) wc/pi h];hwin=eval(win);hw=h.*hwin;
else
n=n-2;
clc,a2=0;a3=2*ap;
while a2<as | a3>ap
n=n+2;
home,disp(['Selecting n = ' int2str(n)]);
n1=1:fix(n/2);h=[sin(wc*n1)./n1/pi];
h=[fliplr(h) wc/pi h];hwin=eval(win);hw=h.*hwin;
w0=[wp ws 0:ws/100:.9*pi];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)=[];
atr=a0(1:2);   %%%%a2=1.05*a0(2);a3=.9*a0(1);
a2=1.05*a0(2);a3=.95*a0(1);a0(1:2)=[];l0=length(a0);
if n>50 & c~=1,c=input('n>50, Enter 1 to continue!  ');clc
if c~=1,break,end,end
end
end

inp=1;while inp~=0
%hw=hw/sum(hw);
if ty=='lp',hwn=hw;end,n0=(n-1)/2;nn=-n0:n0;
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]);
num=int2str(n);atp=num2str(max(atr));ats=num2str(min(atr));

plot(f0,a0,ftr,atr,'o'),axis([0 .5 y 0]);grid
xlabel(g),ylabel('Magnitude dB')
title([wind ' window: n = ' num ' Fc = ' fcn ' ap = ' atp ' as = ' ats]),
if ml==3,pause,end
disp([wind ' window: n = ' num ' Fc = ' fcn ' ap = ' atp ' as = ' ats])
inp=input('enter new n or press ENTER to stop :');
if isempty(inp),break,else,n=inp;clc
if rem(n,2)==0,n=n+1;end
n1=1:fix(n/2);h=[sin(wc*n1)./n1/pi];
h=[fliplr(h) wc/pi h];hwin=eval(win);hw=h.*hwin;
end
end
hold off

%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,xl='XLim';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),grid,axis(bx);
xlabel(g),ylabel('Phase deg.')
elseif ty=='lp',%ax=[0 w1/2/pi ma(1)/mx 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,'*'),axis(ax);grid,xlabel(g),ylabel('Passband Magnitude')
else
ax=[w1/2/pi min(ma(1:2)/mx) 1];axis(ax);
plot(f0,mag,ftr,ma,'o'),axis(ax);grid,xlabel(g),ylabel('Passband Magnitude')
end
subplot
if ml==3,
title([wind ' window: n = ' num ' Fc = ' fcn ' ap = ' atp ' as = ' ats]),
hold off,pause,end
qw=[n;wc/2/pi;abs(max(atr));abs(min(atr));p1];