firdiff.m

很多matlab的源代码

function hw = firdiff(n,fc,wind,p1)
% FIRDIFF FIR differentiator design using windows.
%
%	HW = FIRDIFF(N,FC,WIN,P1) Designs an FIR differentiator
%       N = filter length
%	FC = DIGITAL cutoff frequency [Default: FC = 0.5]
%	WIN = window name as a string [Default: WIN = 'hamming']. 
%	P1 = optional parameter for some windows. Type HELP WINDOW for info 
%	HW returns the coefficients of the designed filter.
%
%	The program allows the user to change N and FC if desired.
%
%    	FIRDIFF (with no input arguments) invokes the following example:
%
%	% Design a 13-point FIR differentiator with FC =0.5 and vonHann window
%      	  >>hdif = firdiff(13,0.4,'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 firdiff,disp('Strike a key to see results of the example')
disp('NOTE: This routine returns control to you. To make no changes,')
disp('Hit the enter key at each prompt')
pause,hdif=firdiff(13,0.4,'vonhann'),return,end

if exist('version')==5,ml=4;else,ml=3;end
if nargin<2,fc=0.5;end
if fc>0.5,error('Digital cutoff frequency must be <= 0.5'),return,end

if nargin<3,wind='hamming';end
if nargin==4,win='windo(wind,n,p1)';
else,win='windo(wind,n)';end
wnd=wind(1:4);odd=0;if wnd=='dolp' | wnd=='cheb';odd=1;end

disp('You may now manually tweak the cutoff and/or change n')
w=[0:pi/200:pi];f=w/2/pi;g='Digital Freq F';
 vx=matverch;
 if vx < 4, eval('clg');else,eval('clf');end

if ml==4,plot([0 0.5],[0 pi],'--'),axis([0 0.5 0 3.5]);hold on
eval('grid on')
hh=0*f;e1='erasemode';x1='xor';plt=eval('plot(f,abs(hh),e1,x1)');
end
while ~isempty(n),
num=int2str(n);hwin=eval(win);
n0=fix(n/2);n1=1:n0;if rem(n,2)==0,n1=n1-.5;end

wx=2*pi*fc;while wx~=0
h=2*fc*cos(2*pi*n1*fc)./n1 - sin(2*pi*n1*fc)/pi./n1./n1;
if rem(n,2)==1,h=[-fliplr(h) 0 h];else,h=[-fliplr(h) h];end
hw=h.*hwin;

cut=num2str(fc);hh=polyval(hw,exp(sqrt(-1)*w));
if ml==3,plot(f,abs(hh),[0 0.5],[0 pi]),grid,end
if ml==4,xd='xdata';yd='ydata';eval('set(plt,xd,f,yd,abs(hh));drawnow');end
xlabel(g),ylabel('Magnitude')
title(['FIR Differentiator: Window=' wind '  n=' num '  fc=' cut ])
if ml==3,pause,end
wx1=input('enter new fc or press ENTER to change order :');
if isempty(wx1),break,else,fc=wx1;end
end
n2=input('enter new n or press ENTER to stop :');
if isempty(n2),hold off,break,end
if rem(n2,2)==0,if odd==1,n2=n2+1;end,end
n=n2;
end
return