firmf.m

很多matlab的源代码

function hn = firmf(fp,fs,pl)
%FIRMF	Maximally flat FIR filter design
%
%	HN=FIRMF(FP,FS,pl) MAXIMALLY FLAT FIR filter impulse response HN
%	FP and FS = DIGITAL frequencies over which response drops from 
%	0.95 and 0.05 (Corresponds to attenuations Ap=0.4455dB, As=26.0206dB)
%	Pl= plotting argument. Pl=0 omits plots of hn and spectrum (DEF. pl=1)
%
%	FIRMF (with no input arguments) invokes the following example
%
%	% Design a maximally flat FIR filter with band edges FP=0.2, FS=0.4
%	  >>hmf=firmf(0.2,0.4);


% 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 firmf,disp('Strike a key to see results of the example')
pause,hmf=firmf(0.2,0.4),return,end

if nargin<3,pl=1;end 
delta=abs(fs-fp);beta=0.5*(fs+fp);n=ceil(1+.5/delta/delta);
if rem(n,2)==0,n=n+1;end
alpha=cos(pi*beta);alpha=alpha*alpha;
np=((n-1)/2:n-1);
for k=1:n;m(k)=min(abs(alpha-k./np));end
k=find(m==min(m));k=min(k);
np=round(k/alpha);l=np-k;n=2*(k+l)-1;

nn=0:l-1;bi=bicoeff(k-1+nn,nn);jj=0;
for f=(0:n-1)/n;
jj=jj+1;cs=(cos(pi*f))^(2*k);kk=0;
hw(jj)=cs*sum((sin(pi*f).^(2*nn)).*bi);
end
ht=ifft(hw);htr=real(ht);
hn=[fliplr(htr(2:(n+1)/2)) htr(1:(n+1)/2)];

if pl==0,return,end
dtplot(0:n-1,hn,'.')
title(['h[n] vs n,  length n=' int2str(n)]),pause
f1=[fp fs];f2=0:.02:0.5;
h=polyval(hn,exp(sqrt(-1)*2*pi*[f1 f2]));
h=20*log10(abs(h));h1=h(1:2);h(1:2)=[];
plot(f2,h,f1,h1,'o'),
grid,title('Magnitude Spectrum, dB vs f')
if exist('version')~=5,hold off,pause,end