% Skript sgolay_2.m, in dem die Glättung eines Spektrums mit
% dem Savitzky-Golay Filter untersucht wird

clear;
%-------- Signal in Form eines Spektogramms
nl = 1000;
x = 0:nl-1;
x0 = [200, 500, 650, 750, 850, 950];     % Stellen der Spektrallinien
alpha = [0.0002, 0.0006, 0.0018, 0.054, 0.108, 0.324]; 
                       % Steilheit der Linien
nk = length(x0);
y = zeros(1,nl);
for k = 1:nk,          % Bildung des Spektrums
   y = y + 1./(1+((x-x0(k)).^2)*alpha(k));
end;
noise = 0.02;             % Varianz des Rauschens
yr = y + sqrt(noise)*randn(1,nl);  % Spektrum mit Rauschanteil
figure(1);    clf;   % -------------
subplot(221), plot(x,y);
title('Ideales Spektrum');
xlabel('n');  grid;
subplot(222), plot(x,yr);
title('Verrauschtes Spektrum');
xlabel('n');  grid;
%------- Savitzky-Golay-Filter
k = 3;      % Grad des Interpolationspolynom
f = 17;     % MR+ML+1 Länge des Filters (ungerade Zahl)
% k = 4;      % Grad des Interpolationspolynom
% f = 33;     % MR+ML+1 Länge des Filters (ungerade Zahl)
yf = sgolayfilt(yr,k,f);  % Filterung mit Savitzky-Golay Filter

subplot(212), plot(x,[yf' y']);
title('Spektrum ohne Rauschanteil und gefiltertes Spektrum');
xlabel('n');  grid;   axis tight;

% ------- Einheitspulsantwort und Frequenzgang des Filters
[b, g] = sgolay(k,f); % Filtermatrix in b und Ableitungsfilter in g
hm = b((f-1)/2+1,:);    % Filter für die mittlere Sequenz
figure(2);    clf;
subplot(211), stem(0:length(hm)-1, hm);
title('Einheitspulsantwort des mittleren Filters');
xlabel('n');  grid;    axis tight;

Nfft = 256;
H = fft(hm,Nfft);
subplot(212), plot((0:Nfft-1)/Nfft, 20*log10(abs(H)));
title(' Amplitudengang des mittleren Savitzky-Golay Filters in dB');
xlabel('f/fs');  grid;
La = axis;    axis([La(1:2), -50, La(4)]);

% -------  Einheitspulsantworten der Ableitungsfilter
figure(3);    clf;
subplot(211), stem(0:f-1, g(:,2));
title('Einheitspulsantwort des Filters für die erste Ableitung');
xlabel('n');    grid on;    axis tight;    
subplot(212), stem(0:f-1, g(:,3));
title('Einheitspulsantwort des Filters für die zweite Ableitung');
xlabel('n');    grid on;    axis tight;

% ------- Geschätzte Ableitung des Spektrums
% Erste und zweite Ableitung mit Filter g(:,2) und g(:,3);
yabl1 = filter(g(:,2),1,y);  % Erste Ableitung des Spektrums ohne Rauschen
yabl1r = filter(g(:,2),1,yr);% Erste Ableitung des Spektrums mit Rauschen
yabl2 = filter(g(:,3),1,y);  % Zweite Ableitung des Spektrums ohne Rauschen
yabl2r = filter(g(:,3),1,yr);% Zweite Ableitung des Spektrums mit Rauschen
figure(4);    clf; % -------------
subplot(211), plot(0:nl-1, yabl1);
title('Erste Ableitung des Spektrums mit und ohne Rauschen');
xlabel('n');  grid;
hold on;  plot(0:nl-1, yabl1r);
hold off;
subplot(212), plot(yabl2);
title('Zweite Ableitung des Spektrums mit und ohne Rauschen');
xlabel('n');  grid;
hold on;  plot(0:nl-1, yabl2r);
hold off;



