psdwelch.m

很多matlab的源代码

function [psd,sl] = psdwelch(x,m,ol,win,p1)
% PSDWELCH Welch method for PSD estimate.
%
%	[P, sl] = PSDWELCH(X,M,OL,WIN,P1) PSD estimate using Welch method
%	X = signal vector/array
%	M =  number of points per section [Default: 0.1*LENGTH(X)], 
%	OL = percent overlap between sections [Default: 50]
%	WIN = string variable for window type [Default: 'hamming'] 
%	P1: Optional parameter for some windows. Type HELP WINDOW for more 
%	P returns the PSD estimate in a length M vector.
%	sl returns the FFT spectrum of each windowed section.
%
%	NOTE: Units of P are such that sum(X.*X)/length(x) = sum(P)/length(P)
%
%       If no output argument specified, PSD vs digital frequency is  plotted.
%
%       PSDWELCH (with no input arguments) invokes the following example:
%
%	% Plot the PSD of 1000 samples of a chirp between F=0.2 and F=0.4 
%       % using 64 points per section, 50% overlap and a vonHann window
%        >>z = chirp([0.2 0.4],1000);
%        >>psdwelch(z,64,50,'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 psdwelch,disp('Strike a key to see results of the example')
pause,z=chirp([0.2 0.4],1000);
 vx=matverch;
 if vx < 4, eval('clg');else,eval('clf');end
psdwelch(z,64,50,'vonhann');return,end

[mx,nx]=size(x);
if mx==1,x=x(:);n=nx;else,n=mx;end
if nargin<4,win='hamming';end,
if nargin<3,ol=50;end,
if nargin<2,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==5,w=window(win,m,p1);
else w=window(win,m);
end,
w=w(:);
%FOR PSD
p=zeros(m,1);
for i=1:k,xw=x(l).*w;l=l+(m-s);xf=abs(fft(xw));sl(i,:)=xf.';p=p+xf.*xf;end
u=sum(w.*w);
%psd=p/(k*m*u); %%OLD
psd=p/(k*u);%%MATCHES SPECTRUM
%psd=2*sqrt(2*psd/m);psd(1)=psd(1)/2; %Normalized 
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,dtplot(f,ps,'.'),hold on,plot(f,ps,':'),hold off
else,plot(f,ps,f,ps,'.'),end,title('PSD vs digital frequency F'),
end