firhilb.m

很多matlab的源代码

function hw = firhilb(n,fc,wind,p1)
% FIRHILB FIR Hilbert Transformer design with windows.
%
%	HW = FIRHILB(N,FC,WIN,P1) Designs a FIR Hilbert Transformer
%       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 = Coefficients of the designed filter.
%
%	Program allows user to change N and FC when prompted.
%
%	FIRHILB (with no input arguments) invokes the following example:
%
%	% Design a 14-point Hilbert transformer with Fc=0.5 and Hamming window
%    	  >>hhilb=firhilb(14,0.5,'hamm')


% 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 firhilb,disp('Strike a key to see results of the example')
disp('NOTE: This program returns control to you. To make no changes,')
disp('Hit the enter key at each prompt')
pause,hhilb=firhilb(14,0.5,'hamm'),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],[1 1],'--'),axis([0 0.5 0 1.2]);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=(1-cos(2*pi*n1*fc))./n1/pi;
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],[1 1],'--'),grid,end
if ml==4,xd='xdata';yd='ydata';eval('set(plt,xd,f,yd,abs(hh));drawnow');end
xlabel(g),ylabel('Magnitude')
title(['Hilbert Transformer: 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