% =========================================================
% Demo: Korrelation
% Modul DSV1
% Hanspeter Hochreutener, 31.3.2007
% =========================================================
% Wenn nach der erten Grafik die Tasten "j" und "ENTER" gedrückt werden,
% wird der Korrelationsvorgang Schritt für Schritt gezeigt.
% Zum nächsten Schritt: drücken der "LEER"- oder der "ETER"-Taste
% Abbrechen der Demo: auf der Matlab-Kommandozeile CTRL+C drücken

% https://home.zhaw.ch/.../dsv1kap4_digisys_korrelation.m


clear all; clc; close all


% Für die verschiedenen Anwendungen und Beispiele
% jeweils beim entsprechenden Block die Kommentare entfernen

% titel = 'Test-Signale';
% muster = [2 1.75 1.5 1.25 1 0.75 0.5 0.25 0];     % 1 Sägezahn
% signal = [1 1 1 1 1 1];         % 1 Rechteck

titel = 'Ähnlichkeit von zwei Signalen';
muster = [0 1 1 1 1 0 -1 -1 -1 -1];     % 1 Rechteck
signal = sin(2*pi*(0:0.1:2.5));         % 2.5 Sinus-Schwingungen

% titel = 'Phasen- oder Zeit-Verschiebung von zwei Signalen messen';
% muster = sin(2*pi*(0:0.1:1));           % 1 Sinus-Schwingungen
% signal = sin(2*pi*((0:0.1:3)-0.4));     % 3 Sinus um 2/5 Periode Phasen-verschoben
% 
% titel = 'Die Signale sind sich nicht ähnlich';
% muster = [1 1 1 1 1 -1 -1 1 1 -1 1 -1 1];   % Synchronisations-Sequenz
% signal = sin(2*pi*(0:0.2:2.5));         % 2.5 Sinus-Schwingungen
% 
% titel = 'Start-Sequenz finden in einem Bit-Strom';
% muster = [1 1 1 1 1 -1 -1 1 1 -1 1 -1 1];   % Synchronisations-Sequenz
% signal = sign([randn(1,5) muster randn(1,10) muster randn(1,5)]);    % Zufall Synchronisations-Sequenz Zufall
% 
% % randn macht weisses Rauschen (Gauss-verteilt)
% % rand macht normalverteiltes Rauschen, alles Werte kommen gleich häufig
% 
% titel = 'Signal-Sequenz finden in einem verrauschten Signal-Strom';
% muster = [1 1 1 1 1 -1 -1 1 1 -1 1 -1 1];   % Signal-Sequenz
% signal = randn(1,40);                     % Rauschen
% signal(11:23) = signal(11:23)+muster;   % Rauschen und Signal-Sequenz
% 
% titel = 'Autokorrelation eines Rauschsignals = Diracstoss';
% muster = randn(1,500);                   % Weisses Rauschen
% signal = muster;                        % gleiche Sequenz = Autokorrelation
% 
% titel = 'Impulsantwort eines FIR-Bandpass-Filters mit Rauschen messen';
% signal = randn(1,1000);                 % Weisses Rauschen
% b = [1 1 1 1]/4;             % Filter-Koeffizienten: gleitender Mittelwert
% b = [-0.1 0.3 0.7 0.3 -0.1]; % Filter-Koeffizienten: Tiefpass
% b = [0.1 -0.3 0.3 -0.3 0.1]; % Filter-Koeffizienten: Hochpass
% b = [0.0671 -0.3394 0.3250 0.3250 -0.3394 0.0671]; % Bandpass
% muster = filter(b,1,signal);            % Filter-Ausgang berechnen
% Die Korrelation ergibt die Filterkoeffizienten b (= Impulsantwort h),
% falls die Messdauer genügend lange ist.

% titel = 'Impulsantwort eines RC-Tiefpass-Filters mit Rauschen messen';
% signal = randn(1,1000);                 % Weisses Rauschen
% b = 0.2; a = [1 -0.8];                  % Filter-Koeffizienten: RC-Tiefpass-Filter
% muster = filter(b,a,signal);            % Filter-Ausgang berechnen
% % Die Korrelation ergibt die Impulsantwort h (Kontrolle mit impz(b,a)),
% % falls die Messdauer genügend lange ist.

if (length(muster) == length(signal))   % Korrelation berechnen
    korr = xcorr(muster,signal,'biased');   % "normieren", wenn Signale gleich lang sind
else
    korr = xcorr(muster,signal);
end;

xmax = (length(korr)-1)/2;              % Achsen für die Grafiken bestimmen
ymax = 1.2*max([muster signal -muster -signal]);
[ykorrmax, yindexmax] = max(korr);

figure(1);                              % Schlussresultat zeigen
subplot(3,1,1);
plot(0:length(signal)-1, signal, '-x');
title(titel); grid;
axis([-xmax, xmax, -ymax, ymax]);
subplot(3,1,2);
plot(0:length(muster)-1, muster,'-o');
title('Signal x (oben) und Muster m (unten)'); grid;
axis([-xmax, xmax, -ymax, ymax]);
subplot(3,1,3);
n = -xmax:xmax;
plot(n, korr, '-o', n, korr, ':');
title(['Korrelation mit Maximum bei \tau = ' num2str(yindexmax-1-xmax) '\cdotTs']); grid;
axis([-xmax, xmax, -ykorrmax, ykorrmax]);

detail = strtrim(input('Schritt für Schritt ablaufen lassen? j/n ', 's'));

if (detail == 'j' || detail == 'J')     % Korrelation Schritt für Schritt zeigen
    figure(2);
    for k = 0:length(korr)-1
        subplot(2,1,1);
        plot((0:length(muster)-1)+xmax-k, muster,'-o', 0:length(signal)-1, signal, '-x');
        title(titel); grid;
        axis([-xmax, xmax, -ymax, ymax]);
        legend(['Muster mit \tau = ' num2str(k-xmax) '\cdotTs'], 'Signal', 'Location', 'SouthWest');
        subplot(2,1,2);
        n = -xmax:xmax;
        plot(n(1:k+1), korr(1:k+1), '-o', n, korr, ':');
        title('Korrelation'); grid;
        axis([-xmax, xmax, -ykorrmax, ykorrmax]);
        if (k == 0)
            text(1,0,'Weiter mit der Leer-Taste');
        end;
        pause()
    end;
end;