% 0. Lesen Sie sich das folgende Matlab-Programm durch... danach: 

% 1. Gehen Sie auf die Seite  http://octave-online.net/
% 2. "Dismissen" Sie die Kommentar-Fenster (kein "Sign in" nötig)
% 3. "Copy&Pasten" Sie den Text dieser Datei in die Kommandozeile
%    am unteren Ende der Seite 
% 4. Drücken Sie "Enter" und warten Sie die Ergebnisse ab
% 5. Scrolen Sie durch die Ergebnisse (angefangen bei dem blauen Plot)



clear; %Neustart
% Bei einer NMR-Messung wird ein Signal gemessen, das eine 
% Überlagerung von zwei Schwingungen mit verschiedenen 
% Frequenzen und Phasenverschiebungen ist

frequenz1  =  8; %ganzzahlige Frequenzen!
frequenz2  =  3;
phasenversch1 = -0.3; % die Signale setzen zeitversetzt ein
phasenversch2 = 1.0;

% Wichtig ist nun, dass die beiden Schwingungen nicht
% gleichgewichtet auftreten, da z.B. von der einen Substanz
% bei dem NMR-Experiment mehr vorliegt als von der anderen.

anteil1 = 1; 
anteil2 = 3;

% Bei dem NMR-Experiment wird regelmäßig im Zeit-Interval 
% t=[0, 2*pi] in n Zeitschritten abgetastet

nSchritte = 80;
t = ((1/nSchritte:1/nSchritte:1)*2*pi)'; 
% Matlab-Schreibweise für die Abtastpunkte

% Es ergibt sich als Zeitsignal b(t) insgesamt
b =  anteil1 * cos(frequenz1*(t + phasenversch1)) ...
    + anteil2 * cos(frequenz2*(t + phasenversch2)); 

% Plot des "exakten" Experiments erstellen für die Ergebnisse 
% der Abtastungen: "Zeit" vs "Stärke des Signals (und Richtung) im NMR"
figure(1); hold on;
plot(t, b,'-'); % Die blaue Linie zeigt die exakte Messung


% Messungenauigkeiten als Gauß-Verteilung aufaddieren. 
% Varianz in Signalstärke = 1% vom Maximal-Ausschlag  
b = b + randn(nSchritte,1)*max(b)*sqrt(0.01)

% Varianz in der Zeitachse = 0.01%
t = t + randn(nSchritte,1)*sqrt(0.0001);


% Plot des "ungenauen" Experiments erstellen mit roten Sternchen
plot(t, b,'r*'); % Diese Messpunkte gelten jetzt als Grundlage 


%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% EIGENTLICHE AUFGABE!!!

% NUN SOLLEN AUS DIESEM SIGNAL DIE ANTEILE DER SUBSTANZEN
% ERMITTELT WERDEN... NUR DIE FREQUENZEN SEIEN BEKANNT
% NICHT JEDOCH DIE PHASENVERSCHIEBUNGEN

% Mathematsche Überlegungen:

% b ist der Realteil von 
% anteil1 * exp(i*frequenz1*(t+phasenversch1)) +
% anteil2 * exp(i*frequenz2*(t+phasenversch2))
% wegen der Eulerformel cos(a) = Realteil(exp(i*a))

% also ist b der Realteil von
% anteil1 * exp(i*phasenversch1*frequenz1) * exp(i*frequenz1*t) +
% anteil2 * exp(i*phasenversch2*frequenz2) * exp(i*frequenz2*t)
% wegen Potenzgesetzen exp(a+b) = exp(a) * exp(b)

% Jetzt fasst man die t-unabhängigen Faktoren zusammen:
% w1 = anteil1 * exp(i*phasenversch1*frequenz1)
% w2 = anteil2 * exp(i*phasenversch2*frequenz2)

% und trennt diese von den t-abhängigen Faktoren...
% a1 = exp(i*frequenz1*t) und a2 = exp(i*frequenz2*t)
% In Matrixschreibweise ist b der Realteil von A(t)*w, wobei
% für A(t) gilt:

A=[exp(i*frequenz1*t), exp(i*frequenz2*t)]

% Da b bereits reell ist und da eine ganze Periode [0, 2*pi] abgetastet
% wird (Ausgleich der Abweichungen bei ganzzahligen Frequenzen), kann man 
% w approximieren, indem man das überbestimmte Gleichungssystem A*w = b 
% näherungsweise löst

w = A\b % Matlab-Schreibweise für Lösung von A'*A*w=A'*b

% Da der anteil1 eine reelle positive Zahl ist, gilt
% antei1 = abs(anteil1) ... Absolutbetrag

% Da exp(i*phasenversch1*frequenz1) eine komplexe Zahl mit Betrag 1 
% ist, gilt ferner: anteil1 = abs(anteil1) = abs(w1)

% Also bekommt man die approx. Anteile, indem man den Absolutbetrag des 
% w-Vektors rechnet

approx_Anteil1 = abs(w(1))
approx_Anteil2 = abs(w(2))

% Das sind die richtigen Anteile, bis auf einen Skalierungsfaktor. 
% Skalieren muss man noch mit dem maximalen Ausschlag des Zeitsignals
% approx_Anteile = abs(w) * (max(b)/max(real(A*w)))
