%% 6.5.2 Signalanalyse und digitale Filterung

% Der Datensatz wurde mit den folgenden Zeilen erzeugt und gespeichert:

t = 0.001:0.001:1;                                   % Zeitvektor Fs = 1000 Hz
xx = 12 * sin (2*pi*4*t) + ...                       % Sollsignal      4 Hz
      4 * sin (2*pi*12*t);                           %                12 Hz
x =  xx + ...
     20 * cos (2*pi*63*t) + ...                      % Störsignal     63 Hz
     18 * cos (2*pi*137*t) .* sin (2*pi*29*t) +  ... % 137-29 Hz, 137+29 Hz
      4 + randn (size(t));                           % Rauschen

datei = fopen ('ueb_signal.dat', 'w');               % Signal speichern (ASCII)
fprintf (datei, '%5.3f  %6.2f\r\n', [t' x']');       % \r\n = CRLF im DOS-Format
fclose (datei);

save ueb_signal t x                                  % Signal speichern (MAT)

% Hinweis: Die Multiplikation eines 137-Hz mit einem 29-Hz-Signal
% (d.h. Amplitudenmodulation) entspricht einer Faltung im Frequenzbereich.
% Dadurch werden Anteile mit 108 und 166 Hz erzeugt (=137 +/- 29 Hz).

clear

% Lösung der Aufgabe

load ueb_signal.dat                          % Signal laden
t = ueb_signal (:,1)';                       % Zeitvektor
x = ueb_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')