fftpsd.m

一个利用元胞自动机进行生态仿真的matlab程序

function [f,px]=fftpsd(x,semi,m,dt,tit)
%FFTPSD	PSD by FFT
%
%	[f,px]=fftpsd(x,m,semi,dt,tit)
%	f    : frequency
%	px    : psd for x
%
%	x    : signal matrix (each column indicates time-series)
%    semi  : Display semilog-graph when yes otherwise loglog
%	m   : Number of points to be used FFT; this must be power of 2. (default=128)
%	dt   : Sampling priod (Nyquist frequency is 1/(2*dt). If missing dt=0
%	tit	: title for the graphs
%
%	Refer to SPECTRUM
%
%	First Code		February     22 1995          
%	Last Modified	Octber       27 1998
%	
%     Copyright	JAERI-CAIL		T. Suzudo

% Evaluate input
if nargin<4	dt=1;	end
sr = 1/dt;
[n l]= size(x);
if n<l x=x'; [n l]= size(x); end;
if n<128 error(' *** Too short time-series data ***'); end
if nargin<3 m=128; elseif m<128 m=128; end

f=(0:sr/m:sr/2)';
if l>4 error('Too many columns!');  end

% Spectrums
for k=1:l
	p=spectrum(x(:,k),m);
	px(:,k)=p(:,1);
end

% Graphics
fig=figure(gcf+1);
clf;

if ~semi
	if l==1
		loglog(f,px(:,1));
	elseif l==2
		loglog(f,px(:,1),f,px(:,2),'--');
	elseif l==3
		loglog(f,px(:,1),f,px(:,2),'--',f,px(:,3),':');
	else
		loglog(f,px(:,1),f,px(:,2),'--',f,px(:,3),':',f,px(:,4),'-.');
	end
	alter='semilog';
	alt=1;
else
	if l==1
		semilogy(f,px(:,1));
	elseif l==2
		semilogy(f,px(:,1),f,px(:,2),'--');
	elseif l==3
		semilogy(f,px(:,1),f,px(:,2),'--',f,px(:,3),':');
	else
		semilogy(f,px(:,1),f,px(:,2),'--',f,px(:,3),':',f,px(:,4),'-.');
	end
	alter='loglog';
	alt=0;
end

xlabel('Frequency (Hz)');
ylabel('PSD');
if nargin>=5
	if ~isstr(tit) error('Invalid argument ''tit'' must be string'); end
	title(tit);
end

% Close Grpah
clphandle = uicontrol('position',[380 5 50 25],'string','close',...
    'callback','close(gcf);');