tfwelch.m

很多matlab的源代码

function tf = tfwelch(x,y,m,ol,win,p1,p2)
%TFWELCH Transfer function estimate using Welch method
%
%	TF=TFWELCH(x,y,m,ol,win,p1,p2) m-point TF estimate using Welch method.
%	x and y are filter input and output. (Must be same length).
%	m = Number of points per windowed section [Def=0.1*length(x)]
%	ol= Percent overlap (Def=50)[ol=50 is a good starting point]
%	win = window as a string [Def= 'hamming']  
%	p1,p2= additional parameters for some windows. Type HELP WINDOW for more
%	TF = TF estimate with length = 1+fix(m/2)
%
%	abs(TF) is plotted vs digital frequency F if no output argument given
%
%	TF (with no input arguments) invokes the following example: 
%
%	Apply a random noise input to a known filter, find the output and 
%	compare the estimated TF from input/output data with actual TF
%
%	h=0.1*ones(1,10),	%a 10 point MA filter. 
%	xx=randist(1000),	%1000 point gaussian random noise input
%	yy=filter(h,1,xx);	%Filter output
%	he=tfwelch(xx,yy); 	%TF estimate using built in defaults	
%	tfplot('z',h,1,[0 0.5],0,1);l=length(he);  	%Plot actual TF
%	hold on,plot((0:l-1)/2/l,abs(he),'g')		%Overplot estimate


% 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


%NOTE: To find the 2l-point impulse response h, we could, for example, use:
%lh=length(he);h=real(ifft([tf;0;fliplr(conj(tf(2:lh)))]));


if nargin==0,help tfwelch,disp('Strike a key to see results of example'),pause
h=0.1*ones(1,10);xx=randist(1000);yy=filter(h,1,xx);tfplot('z',h,1,[0 0.5],0,1);
he=tfwelch(xx,yy);l=length(he);hold on,plot((0:l-1)/2/l,abs(he),'g'),hold off
title('Actual and estimated (green) transfer function magnitude'),return,end

[mx,nx]=size(x);if mx==1,x=x(:);n=nx;else,n=mx;end,y=y(:);
if length(x)~=length(y),error('x and y must be same length'),return,end
if nargin<5,win='hamming';end,if nargin<4,ol=50;end,
if nargin<3,m=max([16,fix(0.1*n)]);end,if m>n,m=n;end
s=fix(ol*m/100);k=fix((n-s)/(m-s));l=1:m;
if nargin==7,w=window(win,m,p1,p2);
elseif nargin==6,w=window(win,m,p1);
else w=window(win,m);end,w=w(:);
%FOR PSD
pxx=zeros(m,1);pxy=pxx;
for i=1:k,
xw=x(l).*w;yw=y(l).*w;
l=l+(m-s);xf=fft(xw);yf=fft(yw);
%slx(i,:)=abs(xf.');sly(i,:)=abs(yf.');
pxx=pxx+xf.*conj(xf);pxy=pxy+yf.*conj(xf);
end

tf=pxy./pxx;tf=tf(1:fix((m+1)/2));
%tf=pxy./pxx;tf=tf(1:1+fix(m/2));  %%Previous gave error in v5 for even m


if nargout==0,ll=length(tf);plot((0:ll-1)/2/ll,abs(tf))
title('Estimated TF magnitude vs Digital frequency F'),end

%%FOR PSD ONLY
%u=sum(w.*w);psd=pxx./(k*u);%%MATCHES SPECTRUM
%n=m/2;f=((0:m-1)-n)/m;if mx==1,psd=psd.';end
%if nargout==0,ps=fftshift(psd);
%if length(ps)<50,lines(f,ps,'.'),hold on,plot(f,ps,':'),hold off
%else,plot(f,ps,f,ps,'.'),end,title('PSD vs digital frequency F'),end