fftsig.m

结构力学中的有限元例子,包含了7个分类文件夹

function [FT_T,omega] = FFTsig (T,T2,dt)

% Fourier spectrum of a signal (plot of fft, signal)
% INPUT:
%       T,T2  -  input, response signals
%       dt    - sampling period
% OUTPUT:
%       FT_T  - fft
%       omega - frequency axis (rad/s)

if nargin==0 T=rand(2000,1); T2=rand(2000,1);  dt=0.005; end;
z = length(T);
omega = [1:z].*(2*pi)./(z*dt); mF=floor(z/3); % freq axis

FT_T = fft(T); FT_T2 = fft(T2);
fm=(2*pi)/omega(  find(FT_T(:)==max(FT_T)));
fm2=(2*pi)/omega(  find(FT_T2(:)==max(FT_T2)));
figure(390)
subplot(211)
loglog(omega(1:mF)./(2*pi),abs(FT_T(1:mF)),'b-'); hold on;
loglog(omega(1:mF)./(2*pi),abs(FT_T2(1:mF)),'r-'); ylabel('FFT');
xlabel('Frequency, Hz');
title(['max FFT excit./resp.:  ' num2str(fm) '/' num2str(fm2) ' sec']);
hold off;

t=[1:z].*dt;
subplot(212)
plot(t,T2,'r-',t,T,'b-'); legend('response','excitation');
hold on; xlabel('Time, s'); ylabel('Amplitude'); 

index1=find(T==max(T)); index2=find(T==min(T));
PGA=0;
if abs(index1)>abs(index2) index=index1; PGA=max(T); else index=index2; PGA=min(T); end;
plot(t(index),T(index),'ko','MarkerSize',10);
text(t(index)+2,T(index),['peak excit. =  ' num2str(PGA)])
hold off;