source.m 1.24 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69
% The function computes amplitude spectrum of a Ricker, Fumue
% or extern source wavelet (central frequency: fs, spectral sampling: df)

function [amp,f]=source(s,fp1,fp2,df,fs)

if ~strcmp(s,'from_file'),
   tsour=1/fs;

   dt=tsour/256;

  k=5;
   tbeg=0;
   tend=k*tsour;

   t=[tbeg:dt:tend];
   n=length(t);
   lsour=tsour/dt;

   if strcmp(s,'fumue'),
      % FUCHS-MUELLER-SIGNAL
      w1=2*pi/tsour;
      ft=sin(w1*(t-tsour))-0.5*sin(2*w1*(t-tsour));
      ft(1:(k/2-0.5)*n/k)=0.0;
      ft((k/2+0.5)*n/k:n)=0.0;
   elseif strcmp(s,'ricker'),
      % RICKER-SIGNAL:
      t0=tsour*1.5;
      T0=tsour*1.5;
      tau=pi*(t-t0)/T0;
      a=4;
      ft=(1-a*tau.*tau).*exp(-2*tau.*tau);
   end

else
   load wavelet.dat
   ft=wavelet(2:size(wavelet,1))';
   dt=1/fs;
end

ft=ft-mean(ft);

% Graphische Darstellung der eingelesenen Zeitreihe
figure;
plot(t,ft)
title(['Eingelesene Zeitreihe']);
xlabel('Zeit [s]');
ylabel('Amplitude');



fnyq=1/(2*dt);
nfft=2^nextpow2(length(ft))

%FFT
y=fftshift(fft(ft,nfft));
amp=abs(y);
amp=amp/max(amp);
f=fnyq*(-nfft/2:nfft/2-1)/(nfft/2);

nn1=(nfft/2+1)+(fp1/fnyq)*nfft/2 ;
nn2=(nfft/2+1)+(fp2/fnyq)*nfft/2;

figure;
plot(f(nn1:nn2),amp(nn1:nn2))
title('Amplitudenspektrum')
xlabel('Frequenz [Hz]')
ylabel('Amplitude');