load Signal.dat                              % Signal laden
t = Signal (:,1)';                           % Zeitvektor
x = Signal (:,2)';                           % Messdaten-Vektor

T = diff (t(1:2));                           % Abtastzeit
N = length (x);                              % Länge des Datenvektors
f = [0:floor((N-1)/2)] / (N*T);              % Frequenzvektor für Plot

X = fft (x);                                 % Fouriertransformation
X = [X(1) 2*X(2:floor((N-1)/2)+1)] / N;      % Begrenzung und Normierung

                                             % FIR-Filter auslegen
ord    = 40;                                 % Ordnung des FIR-Filters
B      = fir1 (ord, 0.04);                   % FIR-Filter berechnen f_g = 20 Hz
x_filt = filter (B, 1, x);                   % Daten filtern
x_filt = x_filt - mean (x_filt);             % Gleichanteil ausblenden

figure
subplot (221)
    plot (t, x, 'b-')
    xlabel ('Zeit [s]')
    title ('Ungefiltertes Signal')
subplot (222)
    stem (f, abs (X), 'bo')
    axis ([0 500 0 25])
    xlabel ('Frequenz [Hz]')
    title ('Spektrum')
subplot (223)
    plot (t, x_filt, 'r-')
    xlabel ('Zeit [s]')
    title ('Gefiltertes Signal')

% Die Übertragungsfunktion des FIR-Filters kann angezeigt werden mit:

F   = [-ord/2:0.1:ord/2] / (ord*T);
amp = zeros (size (F));
for i = -ord/2:ord/2
   amp = amp + B(i+ord/2+1) * cos (2*pi*F*i*T);
end
subplot (224)
    plot (F, abs (amp))
    axis ([0 500 0 1])
    xlabel ('Frequenz [Hz]')
    title ('Filter-Übertragungsfunktion')