//Datei FFT.sci
//geschrieben für das Open Source Programm scilab 5.2, www.scilab.org

function drawFFT(Datei, Titel, Abtastzeit, Pfad, Einheit, varargin)
//Routine zum Erstellen eines Frequenzspektrum-Diagramms
//Übergabewerte:
//Datei: Pfad zur Datei mit den Messwerten (in Anführungszeichen)
//Titel: Name der Messung (in "), darf keine . enthalten
//Abtastzeit: zeitlicher Abstand zwischen zwei Messungen in ms
//Pfad: String für den Speicherort der Ausgabe, z.B. 'c:\Daten\'
//Einheit: String mit der Masseinheit der aufgenommenen Werte
//Xmax (optional): Frequenz in Hz, bis zu welcher das Diagramm
//gezeichnet wird.

//Die Routine speichert das Diagramm in
//[Pfad]\[Titel] als .pdf und .gif

//.pdf-Dateien können direkt in pdfTeX eingebunden werden oder mit
//einem PDF-Reader in voller Qualität gedruckt werden.
	
//lese Vektor aus Pfad Datei ein, merke seine Länge
f = read(Datei, -1, 1);
laengeVonF = length(f);
	
//berechne ein Kaiser-Fenster mit alpha=6, ähnlich Hanning
Fenster = window ('kr', length(f), 6);
  
//appliziere das Fenster
for i = 1:length(f), f(i) = f(i) * Fenster(i); end

//fülle den Vektor bis zur nächsten Zweierpotenz mit
//Nullen auf (zero padding). Dadurch verbessert sich die spektrale
//Auflösung und es kann der echte FFT-Algorithmus genommen werden.
NaechsteZweierp = 1;
while 2^(NaechsteZweierp - 1) < length(f),
	NaechsteZweierp = NaechsteZweierp + 1; end
for i = length(f) + 1:2^NaechsteZweierp, f(i) = 0; end

//berechne die Fouriertrafo
c = fft(f, -1);

//skaliere die x-Achse, Einheit Hz
n = 0:length(c) - 1;
n = n / length(c) * 1000 / Abtastzeit;

//Resultat unabhängig von der Anzahl Messungen und vom Fenster machen
c = c * 2 / laengeVonF;

//nimm Betrag, Phasenverschiebungen interessieren nicht
c = abs(c);
	
//begrenze die Frequenz, bis zu welcher das Diagramm gezeichnet wird,
//auf die halbe Maschinenfrequenz, lasse schmalbandigere Diagramme zu
if isempty(varargin) | varargin(1) > 1000 / Abtastzeit / 2 then
	varargin(1) = 1000 / Abtastzeit / 2;
end
	
//ein brachiales Tiefpassfilter...
laenge = length(c); [maximum] = max(c);
for i = 1:0.002 * varargin(1) * Abtastzeit / 1000 * length(c) + 1,
	c(i) = maximum / 1000000;
	c(laenge - i + 1) = maximum / 1000000;
end
	
//lösche alle Koeffizienten ausserhalb des darzustellenden Bereichs,
//damit die y-Achse richtig skaliert wird
for j = varargin(1) * Abtastzeit / 1000 * length(c) + 2:length(c),
	c(j) = 0;
end,

//Säulendiagramm erstellen und formatieren
plot2d3 (n, c, rect = [0, 0, varargin(1)*1.005, max(c)*1.04]);
ordina = 'Freq. / Hz';
xtitle ('FFT ' + Titel, ordina, 'Amplitude / ' + Einheit);
t = get('hdl');
t.font_size = 3;
t.title.font_size = 4;
t.x_label.font_size = 3;
t.y_label.font_size = 3;
t.sub_ticks = [9, 3];

//speichere Plot im angegebenen Verzeichnis als .pdf und .gif
//entferne zuerst Leerzeichen und Punkte aus dem Dateinamen
Titel = strsubst(strsubst(Titel, ' ', '_'), '.', '_') ;
xs2gif(0, Pfad + Titel + '.gif');
xs2pdf(0, Pfad + Titel + '.pdf');
	
endfunction

function y = arg(x)
//berechnet aus einem komplexen Spaltenvektor x einen
//Spaltenvektor mit Winkeln in Grad als Rückgabewert
	for j = 1:length(x)
	 if real(x(j)) == 0  then
	   if imag(x(j)) == 0 then y(j) = 0;
	   elseif imag(x(j)) < 0 then y(j) = -90;
	   elseif imag (x(j)) > 0 then y(j) = 90; end,
	 else
	   y(j) = atan(imag(x(j)) / real(x(j))) / (2*%pi) * 360;
	 end,
	end
endfunction

function x = Frequenzspitze(Datei, Abtastzeit, varargin)
//Diese Funktion gibt für eine Datei, welche zeilenetrennte Messwerte
//enthält, den Betrag, die Phasenlage in °, die Frequenz, die Fre-
//quenz als Vielfaches der Abtastfrequenz und den Index des Fourier-
//koeffizienten mit dem grössten Betrag im Frequenzspektrum der
//Messung zurück.
//Der Index beginnt bei 0 (arithmetisches Mittel)
//Die ersten Koeffizienten werden nicht berücksichtigt (DC-Anteil)
//optional kann ein Parameter minFrequ und maxFrequ angegeben werden,
//um nach Nebenfrequenzspitzen zu suchen.

//Die Rückgabe erfolgt als Zeilenvektor in obgenannter Reihenfolge

//lese Vektor aus Pfad Datei ein, merke seine Länge
f = read(Datei, -1, 1);
laengeVonF = length(f);
	
//berechne ein Kaiser-Fenster mit alpha=6, ähnlich Hanning
Fenster = window ('kr', length(f), 6);
  
//appliziere das Fenster
for j = 1:length(f), f(j) = f(j) * Fenster(j); end

//fülle den Vektor bis zur nächsten Zweierpotenz mit Nullen auf
//(zero padding). Dadurch verbessert sich die spektrale Auflösung
//und es kann der echte FFT-Algorithmus genommen werden.
NaechsteZweierp = 1;
while 2^(NaechsteZweierp - 1) < length(f),
 NaechsteZweierp = NaechsteZweierp + 1;
end
for i = length(f) + 1:2^NaechsteZweierp, f(i) = 0; end

//berechne die Fouriertrafo
c = fft(f, -1);

//Resultat unabhängig von der Anzahl Messungen und vom Fenster machen
c = c * 2 / laengeVonF;

//ein brachiales Tiefpassfilter...
laenge = length(c); maximum = 1;
for i = 1:round(0.002 * laenge),
 c(i) = maximum / 10000;
 c(laenge - i + 1) = maximum / 10000;
end

//sind die optionalen Parameter Xmin und Xmax gesetzt, setze alle
//Koeffizienten ausserhalb dieses Intervalles = 0
if ~isempty (varargin) then
 for j = 1:varargin(1) * Abtastzeit / 1000 * length(c) + 1,
  c(j) = 0;
 end,
 for j = varargin(2) * Abtastzeit / 1000 * length(c) + 1:length(c),
  c(j) = 0;
 end,
else
 //sind die Parameter nicht gesetzt, lösche die rechte Seite der FFT
 for j = round(length(c) / 2) + 2:length(c), c(j) = 0; end,
end
  
//suche grössten Koeffizienten, setze Rückgabewerte
  [Betrag, i] = max( abs(c));
  x(1,1) = Betrag;
  x(1,2) = arg(c(i));
  x(1,3) = (i - 1) / length(c) * 1000 / Abtastzeit;
  x(1,4) = (i - 1) / length(c);
  x(1,5) = i - 1;
endfunction