summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--epost.lpi10
-rw-r--r--epost.lps172
-rw-r--r--epostunit.pas320
-rw-r--r--typenunit.pas1
-rw-r--r--werteunit.pas362
5 files changed, 273 insertions, 592 deletions
diff --git a/epost.lpi b/epost.lpi
index dd933d7..3c53718 100644
--- a/epost.lpi
+++ b/epost.lpi
@@ -39,7 +39,7 @@
<PackageName Value="LazUtils"/>
</Item1>
</RequiredPackages>
- <Units Count="6">
+ <Units Count="8">
<Unit0>
<Filename Value="epost.lpr"/>
<IsPartOfProject Value="True"/>
@@ -64,6 +64,14 @@
<Filename Value="typenunit.pas"/>
<IsPartOfProject Value="True"/>
</Unit5>
+ <Unit6>
+ <Filename Value="../units/fftunit.pas"/>
+ <IsPartOfProject Value="True"/>
+ </Unit6>
+ <Unit7>
+ <Filename Value="../units/fftunit.inc"/>
+ <IsPartOfProject Value="True"/>
+ </Unit7>
</Units>
</ProjectOptions>
<CompilerOptions>
diff --git a/epost.lps b/epost.lps
index 79154bf..8214172 100644
--- a/epost.lps
+++ b/epost.lps
@@ -3,12 +3,12 @@
<ProjectSession>
<Version Value="9"/>
<BuildModes Active="Default"/>
- <Units Count="15">
+ <Units Count="17">
<Unit0>
<Filename Value="epost.lpr"/>
<IsPartOfProject Value="True"/>
- <TopLine Value="66"/>
- <CursorPos X="46" Y="76"/>
+ <TopLine Value="28"/>
+ <CursorPos X="54" Y="17"/>
<UsageCount Value="202"/>
<Loaded Value="True"/>
</Unit0>
@@ -24,8 +24,8 @@
<IsPartOfProject Value="True"/>
<IsVisibleTab Value="True"/>
<EditorIndex Value="1"/>
- <TopLine Value="3084"/>
- <CursorPos X="45" Y="3088"/>
+ <TopLine Value="2479"/>
+ <CursorPos X="59" Y="2499"/>
<UsageCount Value="201"/>
<Loaded Value="True"/>
</Unit2>
@@ -40,203 +40,215 @@
<Filename Value="werteunit.pas"/>
<IsPartOfProject Value="True"/>
<EditorIndex Value="2"/>
- <TopLine Value="695"/>
- <CursorPos X="48" Y="706"/>
+ <TopLine Value="606"/>
+ <CursorPos Y="638"/>
<UsageCount Value="200"/>
<Loaded Value="True"/>
</Unit4>
<Unit5>
<Filename Value="typenunit.pas"/>
<IsPartOfProject Value="True"/>
- <EditorIndex Value="3"/>
- <TopLine Value="836"/>
- <CursorPos Y="868"/>
+ <EditorIndex Value="5"/>
+ <TopLine Value="211"/>
<UsageCount Value="200"/>
<Loaded Value="True"/>
</Unit5>
<Unit6>
+ <Filename Value="../units/fftunit.pas"/>
+ <IsPartOfProject Value="True"/>
+ <EditorIndex Value="3"/>
+ <UsageCount Value="39"/>
+ <Loaded Value="True"/>
+ </Unit6>
+ <Unit7>
+ <Filename Value="../units/fftunit.inc"/>
+ <IsPartOfProject Value="True"/>
+ <EditorIndex Value="4"/>
+ <CursorPos X="20" Y="20"/>
+ <UsageCount Value="36"/>
+ <Loaded Value="True"/>
+ </Unit7>
+ <Unit8>
<Filename Value="../fpGUI/src/corelib/render/software/agg_scanline_storage_aa.pas"/>
<EditorIndex Value="-1"/>
<TopLine Value="1612"/>
<CursorPos X="2" Y="1675"/>
- <UsageCount Value="4"/>
- </Unit6>
- <Unit7>
+ <UsageCount Value="2"/>
+ </Unit8>
+ <Unit9>
<Filename Value="../units/mystringlistunit.pas"/>
<EditorIndex Value="-1"/>
<TopLine Value="415"/>
<CursorPos Y="435"/>
- <UsageCount Value="16"/>
- </Unit7>
- <Unit8>
+ <UsageCount Value="14"/>
+ </Unit9>
+ <Unit10>
<Filename Value="../units/lowlevelunit.pas"/>
<EditorIndex Value="-1"/>
- <TopLine Value="852"/>
- <CursorPos Y="892"/>
- <UsageCount Value="26"/>
- </Unit8>
- <Unit9>
+ <TopLine Value="2"/>
+ <CursorPos X="3" Y="22"/>
+ <UsageCount Value="24"/>
+ </Unit10>
+ <Unit11>
<Filename Value="../units/randomunit.pas"/>
<EditorIndex Value="-1"/>
- <UsageCount Value="4"/>
- </Unit9>
- <Unit10>
+ <UsageCount Value="2"/>
+ </Unit11>
+ <Unit12>
<Filename Value="../units/matheunit.pas"/>
<EditorIndex Value="-1"/>
<TopLine Value="544"/>
<CursorPos X="53" Y="567"/>
<FoldState Value=" T3q50{m012A"/>
- <UsageCount Value="16"/>
- </Unit10>
- <Unit11>
+ <UsageCount Value="14"/>
+ </Unit12>
+ <Unit13>
<Filename Value="../units/systemunit.pas"/>
<EditorIndex Value="-1"/>
<TopLine Value="65"/>
<CursorPos X="26" Y="80"/>
- <UsageCount Value="17"/>
- </Unit11>
- <Unit12>
+ <UsageCount Value="15"/>
+ </Unit13>
+ <Unit14>
<Filename Value="../fpGUI/src/corelib/render/software/agg_2D.pas"/>
<EditorIndex Value="-1"/>
<TopLine Value="807"/>
<CursorPos Y="818"/>
- <UsageCount Value="7"/>
- </Unit12>
- <Unit13>
+ <UsageCount Value="5"/>
+ </Unit14>
+ <Unit15>
<Filename Value="../units/protokollunit.pas"/>
<EditorIndex Value="-1"/>
<TopLine Value="82"/>
<CursorPos X="15" Y="30"/>
- <UsageCount Value="7"/>
- </Unit13>
- <Unit14>
+ <UsageCount Value="5"/>
+ </Unit15>
+ <Unit16>
<Filename Value="/usr/lib/fpc/src/rtl/inc/objpash.inc"/>
<EditorIndex Value="-1"/>
<TopLine Value="182"/>
<CursorPos X="21" Y="202"/>
- <UsageCount Value="9"/>
- </Unit14>
+ <UsageCount Value="7"/>
+ </Unit16>
</Units>
<JumpHistory Count="30" HistoryIndex="29">
<Position1>
<Filename Value="epostunit.pas"/>
- <Caret Line="1549" Column="35" TopLine="1538"/>
+ <Caret Line="3038" Column="13" TopLine="3006"/>
</Position1>
<Position2>
- <Filename Value="epostunit.pas"/>
- <Caret Line="78" Column="43" TopLine="58"/>
+ <Filename Value="werteunit.pas"/>
+ <Caret Line="652" TopLine="636"/>
</Position2>
<Position3>
- <Filename Value="epostunit.pas"/>
- <Caret Line="1752" Column="14" TopLine="1726"/>
+ <Filename Value="werteunit.pas"/>
</Position3>
<Position4>
- <Filename Value="epostunit.pas"/>
- <Caret Line="1476" TopLine="1456"/>
+ <Filename Value="werteunit.pas"/>
+ <Caret Line="46" Column="17" TopLine="13"/>
</Position4>
<Position5>
- <Filename Value="epostunit.pas"/>
- <Caret Line="1483" Column="15" TopLine="1468"/>
+ <Filename Value="werteunit.pas"/>
+ <Caret Line="47" Column="17" TopLine="14"/>
</Position5>
<Position6>
- <Filename Value="epostunit.pas"/>
- <Caret Line="1740" Column="27" TopLine="1720"/>
+ <Filename Value="werteunit.pas"/>
+ <Caret Line="639" Column="22" TopLine="648"/>
</Position6>
<Position7>
- <Filename Value="epostunit.pas"/>
- <Caret Line="1765" Column="21" TopLine="1734"/>
+ <Filename Value="werteunit.pas"/>
+ <Caret Line="642" Column="16" TopLine="634"/>
</Position7>
<Position8>
- <Filename Value="epostunit.pas"/>
- <Caret Line="1752" Column="26" TopLine="1749"/>
+ <Filename Value="werteunit.pas"/>
+ <Caret Line="644" Column="16" TopLine="634"/>
</Position8>
<Position9>
<Filename Value="epostunit.pas"/>
- <Caret Line="141" Column="25" TopLine="123"/>
+ <Caret Line="4593" Column="34" TopLine="4565"/>
</Position9>
<Position10>
<Filename Value="epostunit.pas"/>
- <Caret Line="2982" Column="23" TopLine="2969"/>
+ <Caret Line="253" Column="10" TopLine="232"/>
</Position10>
<Position11>
<Filename Value="epostunit.pas"/>
- <Caret Line="343" Column="20" TopLine="323"/>
+ <Caret Line="2462" TopLine="2445"/>
</Position11>
<Position12>
<Filename Value="epostunit.pas"/>
- <Caret Line="2982" Column="24" TopLine="2949"/>
+ <Caret Line="2460" Column="30" TopLine="2440"/>
</Position12>
<Position13>
<Filename Value="epostunit.pas"/>
- <Caret Line="3069" Column="24" TopLine="3036"/>
+ <Caret Line="2483" Column="28" TopLine="2459"/>
</Position13>
<Position14>
<Filename Value="epostunit.pas"/>
- <Caret Line="5917" Column="9" TopLine="5902"/>
+ <Caret Line="4570" Column="46" TopLine="4550"/>
</Position14>
<Position15>
<Filename Value="epostunit.pas"/>
- <Caret Line="260" Column="143" TopLine="240"/>
+ <Caret Line="260" Column="59" TopLine="240"/>
</Position15>
<Position16>
<Filename Value="epostunit.pas"/>
- <Caret Line="2567" Column="10" TopLine="2534"/>
+ <Caret Line="4575" Column="12" TopLine="4556"/>
</Position16>
<Position17>
<Filename Value="epostunit.pas"/>
- <Caret Line="2577" Column="10" TopLine="2544"/>
+ <Caret Line="4588" TopLine="4560"/>
</Position17>
<Position18>
<Filename Value="epostunit.pas"/>
- <Caret Line="2578" Column="10" TopLine="2545"/>
+ <Caret Line="4582" Column="46" TopLine="4562"/>
</Position18>
<Position19>
<Filename Value="epostunit.pas"/>
- <Caret Line="2579" Column="10" TopLine="2546"/>
+ <Caret Line="4604" Column="79" TopLine="4578"/>
</Position19>
<Position20>
<Filename Value="epostunit.pas"/>
- <Caret Line="2598" Column="27" TopLine="2566"/>
+ <Caret Line="4608" Column="47" TopLine="4584"/>
</Position20>
<Position21>
<Filename Value="epostunit.pas"/>
- <Caret Line="2599" Column="14" TopLine="2567"/>
+ <Caret Line="2212" TopLine="2192"/>
</Position21>
<Position22>
<Filename Value="epostunit.pas"/>
- <Caret Line="2600" Column="14" TopLine="2568"/>
+ <Caret Line="2233" TopLine="2214"/>
</Position22>
<Position23>
- <Filename Value="epostunit.pas"/>
- <Caret Line="2641" Column="14" TopLine="2609"/>
+ <Filename Value="werteunit.pas"/>
+ <Caret Line="646" Column="3" TopLine="625"/>
</Position23>
<Position24>
- <Filename Value="epostunit.pas"/>
- <Caret Line="2663" Column="39" TopLine="2631"/>
+ <Filename Value="werteunit.pas"/>
+ <Caret Line="47" Column="151" TopLine="29"/>
</Position24>
<Position25>
<Filename Value="epostunit.pas"/>
- <Caret Line="2686" Column="24" TopLine="2654"/>
+ <Caret Line="1779" Column="111" TopLine="1767"/>
</Position25>
<Position26>
<Filename Value="epostunit.pas"/>
- <Caret Line="2730" Column="30" TopLine="2697"/>
</Position26>
<Position27>
<Filename Value="epostunit.pas"/>
- <Caret Line="2964" Column="10" TopLine="2931"/>
+ <Caret Line="113" Column="29" TopLine="104"/>
</Position27>
<Position28>
- <Filename Value="typenunit.pas"/>
- <Caret Line="1443" Column="41" TopLine="1433"/>
+ <Filename Value="epostunit.pas"/>
+ <Caret Line="2313" Column="21" TopLine="2281"/>
</Position28>
<Position29>
- <Filename Value="typenunit.pas"/>
- <Caret Line="172" Column="11" TopLine="152"/>
+ <Filename Value="epostunit.pas"/>
+ <Caret Line="2499" Column="104" TopLine="2467"/>
</Position29>
<Position30>
- <Filename Value="werteunit.pas"/>
- <Caret Line="431" Column="34" TopLine="431"/>
+ <Filename Value="../units/fftunit.pas"/>
+ <Caret Line="174" Column="32" TopLine="84"/>
</Position30>
</JumpHistory>
</ProjectSession>
diff --git a/epostunit.pas b/epostunit.pas
index da18189..9102396 100644
--- a/epostunit.pas
+++ b/epostunit.pas
@@ -5,7 +5,7 @@ unit epostunit;
interface
uses
- Classes, SysUtils, mystringlistunit, werteunit, typenunit, process, lowlevelunit, matheunit;
+ Classes, SysUtils, mystringlistunit, werteunit, typenunit, process, lowlevelunit, matheunit, fftunit;
type
TBmpHeader = packed record
@@ -134,7 +134,6 @@ type
procedure ermittleMinMaxDichten(st: boolean; threads,xmin,xmax,tmin,tmax: longint; symmetrisch: boolean); overload;
procedure gleicheMinMaxDichtenAn(st: boolean; var f: tMyStringlist; symmetrisch: boolean);
function fft(threads: longint; senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fen: tFenster; out pvFehler: extended; Warn: tWarnstufe): boolean; overload;
- function fft(threads,xmin,xmax,tmin,tmax: longint; senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fen: tFenster; out pvFehler: extended; Warn: tWarnstufe): boolean; overload;
function berechneZeitfrequenzanalyse(st: boolean; var f: tMyStringlist; threads: longint; const quelle: tWerte; Warn: tWarnstufe): boolean;
function berechneVerzerrung(st: boolean; var f: tMyStringlist; threads: longint; const quelle: tWerte; Warn: tWarnstufe): boolean;
function berechneIntegral(st: boolean; var f: tMyStringlist; threads: longint; const quelle: tWerte): boolean;
@@ -251,22 +250,24 @@ type
procedure stExecute; override;
end;
tFFTThread = class(tLogThread)
- xMi,xMa,tMi,tMa: longint;
- vo,na: tFFTDatenordnung;
+ smi,sma: longint;
fen: tFenster;
erfolg,sen,inv: boolean;
+ algo: tFFTAlgorithmus;
pW: tWerte;
pvFehler: extended;
- constructor create(werte: tWerte; xMin,xMax,tMin,tMax: longint; senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fenster: tFenster);
+ constructor create(werte: tWerte; smin,smax: longint; senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fenster: tFenster); overload;
+ constructor create(werte: tWerte; smin,smax: longint; senkrecht,invers: boolean; algorithmus: tFFTAlgorithmus; const fenster: tFenster); overload;
destructor destroy; override;
procedure stExecute; override;
end;
tKorrelThread = class(tLogThread)
- wl: tWavelet;
- xMi,xMa,tMi,tMa,tOf,xOf: longint;
- qu,zi: tWerte;
- pvFehler: extended;
- constructor create(quelle,ziel: tWerte; xMin,xMax,tMin,tMax,xOff,tOff: longint; wavelet: tWavelet);
+ wl: tWavelet;
+ xMi,xMa: longint;
+ qu,zi: tWerte;
+ pvFehler: extended;
+ algo: tFFTAlgorithmus;
+ constructor create(quelle,ziel: tWerte; xMin,xMax: longint; wavelet: tWavelet);
destructor destroy; override;
procedure stExecute; override;
end;
@@ -2204,11 +2205,10 @@ end;
function tWerte.berechneKorrelation(st: boolean; var f: tMyStringlist; threads: longint; const quelle: tWerte): boolean;
var
- i,xmin,xmax,tmin,tmax: longint;
+ i: longint;
s: string;
wavelet: tWavelet;
fertig: boolean;
- ausrichtung: word; // 0 = zentriert, 1 = anfangsbündig, 2 = endbündig
korrelThreads: array of tKorrelThread;
Zeit,pvFehler: extended;
pSi: pTLLWerteSingle;
@@ -2227,36 +2227,15 @@ begin
wavelet.tfwhm:=1;
wavelet.typ:=wtSin2;
_xsteps:=quelle._xsteps;
- xmin:=0;
- xmax:=_xsteps-1;
_tsiz:=quelle._tsiz;
- tmin:=0;
- tmax:=_tsiz-1;
_np:=quelle._np;
_beta:=quelle._beta;
- ausrichtung:=0;
repeat
if not f.metaReadln(s,true) then begin
gibAus('Unerwartetes Dateiende!',3);
exit;
end;
if s='Ende' then break;
- if startetMit('xmin:',s) then begin
- xmin:=kont2disk('x',exprtofloat(st,s));
- continue;
- end;
- if startetMit('xmax:',s) then begin
- xmax:=kont2disk('x',exprtofloat(st,s));
- continue;
- end;
- if startetMit('tmin:',s) then begin
- tmin:=kont2disk('t',exprtofloat(st,s));
- continue;
- end;
- if startetMit('tmax:',s) then begin
- tmax:=kont2disk('t',exprtofloat(st,s));
- continue;
- end;
if startetMit('freq:',s) then begin
wavelet.freq:=1/kont2diskFak('t',1/exprtofloat(st,s));
continue;
@@ -2269,23 +2248,9 @@ begin
if not wavelet.setzeTyp(s) then exit;
continue;
end;
- if startetMit('mit FFT',s) then begin
+ if s='mit FFT' then begin
wavelet.mitFFT:=true;
- if s='' then continue;
- if s='zentriert' then begin
- ausrichtung:=0;
- continue;
- end;
- if s='anfangsbündig' then begin
- ausrichtung:=1;
- continue;
- end;
- if s='endbündig' then begin
- ausrichtung:=2;
- continue;
- end;
- gibAus('Kenne Ausrichtung '''+s+''' nicht für die FFT bei einer Korrelationsberechnung!',3);
- exit;
+ continue;
end;
if s='ohne FFT' then begin
wavelet.mitFFT:=false;
@@ -2295,28 +2260,6 @@ begin
exit;
until false;
- Transformationen:=tKoordinatenAusschnitt.create(Transformationen,xmin,xmax,tmin,tmax);
- _xsteps:=xmax-xmin+1;
- if wavelet.mitFFT then begin
- i:=1;
- while 2*i<=tmax-tmin+1 do
- i:=i*2;
- case ausrichtung of
- 0: begin
- tmin:=(tmax+1+tmin-i) div 2;
- tmax:=tmin+i-1;
- end;
- 1: tmax:=tmin+i-1;
- 2: tmin:=tmax-i+1;
- end{of case};
- if (tmin<0) or (tmax>=_tsiz) then begin
- gibAus('Der ausgewählte Bereich liegt irgendwie außerhalb der Grenzen, das sollte nicht passieren! ('+inttostr(tmin)+'-'+inttostr(tmax)+' aus 0-'+inttostr(_tsiz-1)+')',3);
- gibAus('selbst: '+paramsDump,3);
- gibAus('Quelle: '+quelle.paramsDump,3);
- exit;
- end;
- _tsiz:=tmax-tmin+1;
- end;
if st then begin
result:=true;
exit;
@@ -2328,22 +2271,21 @@ begin
case quelle.genauigkeit of
gSingle: begin
pSi:=@(quelle.sWerte);
- eWerte.kopiereVon(st,pSi,xmin,xmax,tmin,tmax);
+ eWerte.kopiereVon(st,pSi);
end;
gDouble: begin
pDo:=@(quelle.dWerte);
- dWerte.kopiereVon(st,pDo,xmin,xmax,tmin,tmax);
+ dWerte.kopiereVon(st,pDo);
end;
gExtended: begin
pEx:=@(quelle.eWerte);
- eWerte.kopiereVon(st,pEx,xmin,xmax,tmin,tmax);
+ eWerte.kopiereVon(st,pEx);
end;
end{of case};
gibAus('... fertig '+timetostr(now-Zeit)+', berechne ...',3);
end
else begin
genauigkeit:=gExtended;
- _tsiz:=tmax+1-tmin;
eWerte.holeRam(3);
gibAus('Berechne ...',3);
end;
@@ -2356,7 +2298,7 @@ begin
end;
setlength(korrelThreads,threads);
for i:=0 to length(korrelThreads)-1 do
- korrelThreads[i]:=tKorrelThread.create(quelle,self,round(i*_xsteps/threads),round((i+1)*_xsteps/threads-1),0,_tsiz-1,xmin,tmin,wavelet);
+ korrelThreads[i]:=tKorrelThread.create(quelle,self,round(i*_xsteps/threads),round((i+1)*_xsteps/threads-1),wavelet);
repeat
sleep(10);
fertig:=true;
@@ -2476,69 +2418,67 @@ end;
function tWerte.fft(threads: longint; senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fen: tFenster; out pvFehler: extended; Warn: tWarnstufe): boolean;
var
- len: longint;
-begin
- len:=1;
- if senkrecht then begin
- while 2*len<=_tsiz do
- len:=len*2;
- result:=fft(threads,0,_xsteps-1,(_tsiz-len) div 2,((_tsiz+len) div 2) - 1,true,invers,vor,nach,fen,pvFehler,Warn);
- end
- else begin
- while 2*len<=_xsteps do
- len:=len*2;
- result:=fft(threads,(_xsteps-len) div 2,((_xsteps+len) div 2) - 1,0,_tsiz-1,false,invers,vor,nach,fen,pvFehler,Warn);
- end;
-end;
-
-function tWerte.fft(threads,xmin,xmax,tmin,tmax: longint; senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fen: tFenster; out pvFehler: extended; Warn: tWarnstufe): boolean;
-var
- fftThreads: array of tFFTThread;
- i: longint;
- fertig: boolean;
+ fftThreads: array of tFFTThread;
+ i: longint;
+ fertig: boolean;
begin
result:=false;
if senkrecht then begin
- if length(fen.Werte)<>tmax+1-tmin then
- fen.berechneWerte(tmax+1-tmin);
+ if length(fen.Werte)<>_tsiz then
+ fen.berechneWerte(_tsiz);
if threads>_xsteps then
threads:=_xsteps;
end
else begin
- if length(fen.Werte)<>xmax+1-xmin then
- fen.berechneWerte(xmax+1-xmin);
+ if length(fen.Werte)<>_xsteps then
+ fen.berechneWerte(_xsteps);
if threads>_tsiz then
threads:=_tsiz;
end;
+
setlength(fftThreads,threads);
if senkrecht then begin
- for i:=0 to threads-1 do
+ fftThreads[0]:=
+ tFFTThread.create(
+ self,
+ 0,
+ round(_xsteps/threads-1),
+ senkrecht,
+ invers,
+ vor,
+ nach,
+ fen);
+ for i:=1 to threads-1 do
fftThreads[i]:=
tFFTThread.create(
self,
- xmin+round((xmax+1-xmin)/threads*i),
- xmin+round((xmax+1-xmin)/threads*(i+1)-1),
- tmin,
- tmax,
+ round(_xsteps/threads*i),
+ round(_xsteps/threads*(i+1)-1),
senkrecht,
invers,
- vor,
- nach,
+ fftThreads[0].algo,
fen);
end
else begin
- for i:=0 to threads-1 do
+ fftThreads[0]:=
+ tFFTThread.create(
+ self,
+ 0,
+ round(_tsiz/threads-1),
+ senkrecht,
+ invers,
+ vor,
+ nach,
+ fen);
+ for i:=1 to threads-1 do
fftThreads[i]:=
tFFTThread.create(
self,
- xmin,
- xmax,
- tmin+round((tmax+1-tmin)/threads*i),
- tmin+round((tmax+1-tmin)/threads*(i+1)-1),
+ round(_tsiz/threads*i),
+ round(_tsiz/threads*(i+1)-1),
senkrecht,
invers,
- vor,
- nach,
+ fftThreads[0].algo,
fen);
end;
repeat
@@ -2958,7 +2898,6 @@ end;
function tWerte.berechneFFT(st: boolean; var f: tMyStringlist; threads: longint; Warn: tWarnstufe): boolean;
var
- i: longint;
Zeit,pvFehler: extended;
NB: tFFTDatenordnung;
Fenster: tFenster;
@@ -3013,26 +2952,6 @@ begin
end;
gibAus(' (Parseval-Fehler = '+floattostr(pvFehler)+')',3);
end;
- if senkrecht then begin
- i:=1;
- while 2*i<=_tsiz do
- i:=i*2;
- Transformationen:=tKoordinatenAusschnitt.create(Transformationen,0,_xsteps-1,(_tsiz-i) div 2,((_tsiz+i) div 2) - 1);
- if not st then begin
- gibAus('Die Länge wird von '+inttostr(_tsiz)+' auf '+inttostr(i)+' Zeitschritte gekürzt!',3);
- _tsiz:=i;
- end;
- end
- else begin
- i:=1;
- while 2*i<=_xsteps do
- i:=i*2;
- Transformationen:=tKoordinatenAusschnitt.create(Transformationen,(_xsteps-i) div 2,((_xsteps+i) div 2) - 1,0,_tsiz-1);
- if not st then begin
- gibAus('Die Länge wird von '+inttostr(_xsteps)+' auf '+inttostr(i)+' Ortsschritte gekürzt!',3);
- _xsteps:=i;
- end;
- end;
Transformationen:=tFFTTransformation.create(Transformationen,not senkrecht,senkrecht);
if not st then begin
eWerte.holeRam(0);
@@ -3043,12 +2962,11 @@ end;
function tWerte.berechneFFT2d(st: boolean; var f: tMyStringlist; threads: longint; Warn: tWarnstufe): boolean;
var
- i,k: longint;
- Zeit,pvFehler: extended;
- NB,preOrd: tFFTDatenordnung;
- Fensters: array[boolean] of tFenster;
- s: string;
- b,spiegeln: boolean;
+ Zeit,pvFehler: extended;
+ NB,preOrd: tFFTDatenordnung;
+ Fensters: array[boolean] of tFenster;
+ s: string;
+ b,spiegeln: boolean;
begin
result:=false;
warteaufBeendigungDesLeseThreads;
@@ -3095,26 +3013,6 @@ begin
exit;
end;
- i:=1;
- while 2*i<=_tsiz do
- i:=i*2;
- if _tsiz>i then begin
- gibAus('Die Länge wird von '+inttostr(_tsiz)+' auf '+inttostr(i)+' Zeitschritte gekürzt!',3);
- Transformationen:=tKoordinatenAusschnitt.create(Transformationen,0,_xsteps-1,(_tsiz-i) div 2,((_tsiz+i) div 2) - 1);
- _tsiz:=i;
- eWerte.holeRam(0);
- end;
- i:=1;
- while 2*i<=_xsteps do
- i:=i*2;
- if _xsteps>i then begin
- gibAus('Die Breite wird von '+inttostr(_xsteps)+' auf '+inttostr(i)+' Ortsschritte gekürzt!',3);
- Transformationen:=tKoordinatenAusschnitt.create(Transformationen,(_xsteps-i) div 2,((_xsteps+i) div 2) - 1,0,_tsiz-1);
- for k:=1 to _tsiz-1 do
- Move(eWerte.werte[k*_xsteps],eWerte.werte[k*i],i*sizeof(extended));
- _xsteps:=i;
- eWerte.holeRam(0);
- end;
Transformationen:=tFFTTransformation.create(Transformationen,true,true);
gibAus('... fertig! '+timetostr(now-Zeit),3);
if spiegeln then begin
@@ -3123,14 +3021,14 @@ begin
gibAus('... fertig! '+timetostr(now-Zeit),3);
end;
gibAus('berechne t-FFT ...',3);
- if not fft(threads,0,_xsteps-1,0,_tsiz-1,true,false,doRes,preOrd,Fensters[true],pvFehler,Warn) then begin
+ if not fft(threads,true,false,doRes,preOrd,Fensters[true],pvFehler,Warn) then begin
gibAus('Es traten Fehler auf!',3);
exit;
end;
gibAus(' (Parseval-Fehler = '+floattostr(pvFehler)+')',3);
gibAus('... fertig! '+timetostr(now-Zeit),3);
gibAus('berechne x-FFT ...',3);
- if not fft(threads,0,_xsteps-1,0,_tsiz-1,false,false,doRes,preOrd,Fensters[false],pvFehler,Warn) then begin
+ if not fft(threads,false,false,doRes,preOrd,Fensters[false],pvFehler,Warn) then begin
gibAus('Es traten Fehler auf!',3);
exit;
end;
@@ -4653,21 +4551,30 @@ end;
// tFFTThread ******************************************************************
-constructor tFFTThread.create(werte: tWerte; xMin,xMax,tMin,tMax: longint; senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fenster: tFenster);
+constructor tFFTThread.create(werte: tWerte; smin,smax: longint; senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fenster: tFenster);
+var
+ tmpFFTAlgo: tFFTAlgorithmus;
+begin
+ if senkrecht then
+ tmpFFTAlgo:=createFFTAlgorithmus(werte._tsiz,vor,nach)
+ else
+ tmpFFTAlgo:=createFFTAlgorithmus(werte._xsteps,vor,nach);
+ create(werte,smin,smax,senkrecht,invers,tmpFFTAlgo,fenster);
+ tmpFFTAlgo.free;
+end;
+
+constructor tFFTThread.create(werte: tWerte; smin,smax: longint; senkrecht,invers: boolean; algorithmus: tFFTAlgorithmus; const fenster: tFenster);
begin
inherited create;
pW:=werte;
- xMi:=xMin;
- xMa:=xMax;
- tMi:=tMin;
- tMa:=tMax;
+ smi:=smin;
+ sma:=smax;
sen:=senkrecht;
inv:=invers;
- vo:=vor;
- na:=nach;
+ algo:=createFFTAlgorithmus(algorithmus);
fen:=fenster;
erfolg:=false;
- gibAus('FFTthread kreiert: '+inttostr(xmin)+'-'+inttostr(xmax)+' '+inttostr(tmin)+'-'+inttostr(tmax),1);
+ gibAus('FFTthread kreiert ('+inttostr(pW._xsteps)+' x '+inttostr(pW._tsiz)+': '+inttostr(smi)+' .. '+inttostr(sma)+')',3);
suspended:=false;
end;
@@ -4678,11 +4585,11 @@ end;
procedure tFFTThread.stExecute;
begin
- gibAus('FFTthread gestartet: '+inttostr(xMi)+'-'+inttostr(xMa)+' '+inttostr(tMi)+'-'+inttostr(tMa)+' ...',1);
+ gibAus('FFTthread gestartet ('+inttostr(pW._xsteps)+' x '+inttostr(pW._tsiz)+': '+inttostr(smi)+' .. '+inttostr(sma)+') ...',1);
case pW.Genauigkeit of
- gSingle: erfolg:=pW.sWerte.fft(xMi,xMa,tMi,tMa,sen,inv,vo,na,fen,pvFehler);
- gDouble: erfolg:=pW.dWerte.fft(xMi,xMa,tMi,tMa,sen,inv,vo,na,fen,pvFehler);
- gExtended: erfolg:=pW.eWerte.fft(xMi,xMa,tMi,tMa,sen,inv,vo,na,fen,pvFehler);
+ gSingle: erfolg:=pW.sWerte.fft(smi,sma,sen,inv,algo,fen,pvFehler);
+ gDouble: erfolg:=pW.dWerte.fft(smi,sma,sen,inv,algo,fen,pvFehler);
+ gExtended: erfolg:=pW.eWerte.fft(smi,sma,sen,inv,algo,fen,pvFehler);
end{of case};
gibAus('... und fertig! ',1);
fertig:=true;
@@ -4749,19 +4656,15 @@ end;
// tKorrelThread ***************************************************************
-constructor tKorrelThread.create(quelle,ziel: tWerte; xMin,xMax,tMin,tMax,xOff,tOff: longint; wavelet: tWavelet);
+constructor tKorrelThread.create(quelle,ziel: tWerte; xMin,xMax: longint; wavelet: tWavelet);
begin
inherited create;
qu:=quelle;
zi:=ziel;
xMi:=xMin;
xMa:=xMax;
- tMi:=tMin;
- tMa:=tMax;
- tOf:=tOff;
- xOf:=xOff;
wl:=wavelet;
- gibAus('Korrelationsthread kreiert: '+inttostr(xmin)+'-'+inttostr(xmax)+' '+inttostr(tmin)+'-'+inttostr(tmax),1);
+ gibAus('Korrelationsthread kreiert: '+inttostr(xmin)+'-'+inttostr(xmax),1);
suspended:=false;
end;
@@ -4777,27 +4680,26 @@ var
in0,out0: boolean;
fenster: tFenster;
tmpW: tWerte;
+ tmpFFTAlgo: tFFTAlgorithmus;
begin
gibAus('Korrelationsberechnungsthread gestartet ...',1);
- gibAus('('+inttostr(xmi)+'-'+inttostr(xma)+' x '+inttostr(tmi)+'-'+inttostr(tma)+'), '+inttostr(wl.werte.params.tsiz),1);
+ gibAus('('+inttostr(xmi)+'-'+inttostr(xma)+' x '+inttostr(qu._tsiz)+'), '+inttostr(wl.werte.params.tsiz),1);
in0:=true;
out0:=true;
pvFehler:=0;
if wl.mitFFT then begin
- tma:=tma-tmi; // wenn FFT, dann wurde das schon außerhalb von "kopiereVon" geändert
- tmi:=0;
for i:=xmi to xma do
- for j:=tmi to tma do in0:=in0 and (zi.eWerte.werte[i+j*zi._xsteps]=0);
+ for j:=0 to qu._tsiz-1 do in0:=in0 and (zi.eWerte.werte[i+j*zi._xsteps]=0);
fenster.aktiv:=false;
-// fenster.Breite:=0;
-// fenster.berechneWerte(tma+1-tmi);
gibAus('FFT berechnen ...',1);
- zi.eWerte.fft(xmi,xma,tmi,tma,true,false,doRes,doResIms,fenster,pvF);
+ tmpFFTAlgo:=createFFTAlgorithmus(zi._tsiz,doRes,doResIms);
+ zi.eWerte.fft(true,false,tmpFFTAlgo,fenster,pvF);
+ tmpFFTAlgo.free;
pvFehler:=pvF+pvFehler;
if wl.typ=wtSin2 then // Das Sin²-Wavelet besteht eigntlich aus zwei!
tmpW:=tWerte.create(zi,xmi,xma);
gibAus('... fertig, punktweise Produkte berechnen ...',1);
- hl:=(tma+1-tmi) div 2; // halbe Länge
+ hl:=qu._tsiz div 2; // halbe Länge
for i:=xmi to xma do begin
zi.eWerte.werte[i]:=zi.eWerte.werte[i]*wl.werte.werte[0]; // f_0
zi.eWerte.werte[i+hl*zi._xsteps]:=zi.eWerte.werte[i+hl*zi._xsteps]*wl.werte.werte[2*hl]; // f_n/2
@@ -4823,14 +4725,15 @@ begin
end;
end;
gibAus('... fertig, iFFT berechnen ...',1);
- zi.eWerte.fft(xmi,xma,0,tma-tmi,true,true,doResIms,doBetrQdr,fenster,pvF);
+ tmpFFTAlgo:=createFFTAlgorithmus(zi._tsiz,doResIms,doBetrQdr);
+ zi.eWerte.fft(xmi,xma,true,true,tmpFFTAlgo,fenster,pvF);
pvFehler:=pvF+pvFehler;
case wl.typ of
wtSin2: begin // Das Sin²-Wavelet besteht eigntlich aus zwei!
- tmpW.eWerte.fft(0,xma-xmi,0,tma-tmi,true,true,doResIms,doBetrQdr,fenster,pvF);
+ tmpW.eWerte.fft(true,true,tmpFFTAlgo,fenster,pvF);
pvFehler:=(pvF+pvFehler)/3;
for i:=xmi to xma do
- for j:=0 to tma-tmi do begin
+ for j:=0 to zi._tsiz-1 do begin
zi.eWerte.werte[i+j*zi._xsteps]:=zi.eWerte.werte[i+j*zi._xsteps]+tmpW.eWerte.werte[i-xmi+j*tmpW._xsteps];
out0:=out0 and (zi.eWerte.werte[i+j*zi._xsteps]=0);
end;
@@ -4839,10 +4742,11 @@ begin
wtFrequenzfenster: begin // Das Frequenzfenster-Wavelet ist nur eines!
pvFehler:=pvFehler/2;
for i:=xmi to xma do
- for j:=0 to tma-tmi do
+ for j:=0 to zi._tsiz-1 do
out0:=out0 and (zi.eWerte.werte[i+j*zi._xsteps]=0);
end;
end{of case};
+ tmpFFTAlgo.free;
gibAus(' (Parseval-Fehler = '+floattostr(pvFehler)+')',1);
gibAus('... fertig',1);
end
@@ -4851,45 +4755,45 @@ begin
gSingle:
for i:=xmi to xma do begin
if (xma-i) mod ((xma-xmi) div 10) = 0 then gibAus(inttostr(i)+'/'+inttostr(xmi)+'-'+inttostr(xma),1);
- for j:=tmi to tma do begin
+ for j:=0 to zi._tsiz-1 do begin
sus:=0;
suc:=0;
- for k:=max(-wl.hlen,-tOf-j) to min(wl.hlen,qu._tsiz-j-tOf-1) do begin
- suc:=suc + qu.sWerte.werte[i+xOf+(j+k+tOf)*qu._xsteps]*wl.werte.werte[(k+wl.hlen)*2];
- sus:=sus + qu.sWerte.werte[i+xOf+(j+k+tOf)*qu._xsteps]*wl.werte.werte[(k+wl.hlen)*2+1];
+ for k:=max(-wl.hlen,-j) to min(wl.hlen,qu._tsiz-j-1) do begin
+ suc:=suc + qu.sWerte.werte[i+(j+k)*qu._xsteps]*wl.werte.werte[(k+wl.hlen)*2];
+ sus:=sus + qu.sWerte.werte[i+(j+k)*qu._xsteps]*wl.werte.werte[(k+wl.hlen)*2+1];
end;
zi.eWerte.werte[i+j*zi._xsteps]:=(sqr(sus)+sqr(suc))/sqr(1+2*wl.hlen);
- in0:=in0 and (qu.sWerte.werte[i+xOf+(j+tOf)*qu._xsteps]=0);
+ in0:=in0 and (qu.sWerte.werte[i+j*qu._xsteps]=0);
out0:=out0 and (zi.eWerte.werte[i+j*zi._xsteps]=0);
end;
end;
gDouble:
for i:=xmi to xma do begin
if (xma-i) mod ((xma-xmi) div 10) = 0 then gibAus(inttostr(i)+'/'+inttostr(xmi)+'-'+inttostr(xma),1);
- for j:=tmi to tma do begin
+ for j:=0 to zi._tsiz-1 do begin
sus:=0;
suc:=0;
- for k:=max(-wl.hlen,-tOf-j) to min(wl.hlen,qu._tsiz-j-tOf-1) do begin
- suc:=suc + qu.dWerte.werte[i+xOf+(j+k+tOf)*qu._xsteps]*wl.werte.werte[(k+wl.hlen)*2];
- sus:=sus + qu.dWerte.werte[i+xOf+(j+k+tOf)*qu._xsteps]*wl.werte.werte[(k+wl.hlen)*2+1];
+ for k:=max(-wl.hlen,-j) to min(wl.hlen,qu._tsiz-j-1) do begin
+ suc:=suc + qu.dWerte.werte[i+(j+k)*qu._xsteps]*wl.werte.werte[(k+wl.hlen)*2];
+ sus:=sus + qu.dWerte.werte[i+(j+k)*qu._xsteps]*wl.werte.werte[(k+wl.hlen)*2+1];
end;
zi.eWerte.werte[i+j*zi._xsteps]:=(sqr(sus)+sqr(suc))/sqr(1+2*wl.hlen);
- in0:=in0 and (qu.dWerte.werte[i+xOf+(j+tOf)*qu._xsteps]=0);
+ in0:=in0 and (qu.dWerte.werte[i+j*qu._xsteps]=0);
out0:=out0 and (zi.eWerte.werte[i+j*zi._xsteps]=0);
end;
end;
gExtended:
for i:=xmi to xma do begin
if (xma-i) mod ((xma-xmi) div 10) = 0 then gibAus(inttostr(i)+'/'+inttostr(xmi)+'-'+inttostr(xma),1);
- for j:=tmi to tma do begin
+ for j:=0 to zi._tsiz-1 do begin
sus:=0;
suc:=0;
- for k:=max(-wl.hlen,-tOf-j) to min(wl.hlen,qu._tsiz-j-tOf-1) do begin
- suc:=suc + qu.eWerte.werte[i+xOf+(j+k+tOf)*qu._xsteps]*wl.werte.werte[(k+wl.hlen)*2];
- sus:=sus + qu.eWerte.werte[i+xOf+(j+k+tOf)*qu._xsteps]*wl.werte.werte[(k+wl.hlen)*2+1];
+ for k:=max(-wl.hlen,-j) to min(wl.hlen,qu._tsiz-j-1) do begin
+ suc:=suc + qu.eWerte.werte[i+(j+k)*qu._xsteps]*wl.werte.werte[(k+wl.hlen)*2];
+ sus:=sus + qu.eWerte.werte[i+(j+k)*qu._xsteps]*wl.werte.werte[(k+wl.hlen)*2+1];
end;
zi.eWerte.werte[i+j*zi._xsteps]:=(sqr(sus)+sqr(suc))/sqr(1+2*wl.hlen);
- in0:=in0 and (qu.eWerte.werte[i+xOf+(j+tOf)*qu._xsteps]=0);
+ in0:=in0 and (qu.eWerte.werte[i+j*qu._xsteps]=0);
out0:=out0 and (zi.eWerte.werte[i+j*zi._xsteps]=0);
end;
end;
diff --git a/typenunit.pas b/typenunit.pas
index a86de1c..b1dfdb4 100644
--- a/typenunit.pas
+++ b/typenunit.pas
@@ -168,7 +168,6 @@ type
end;
pTBeschriftungen = ^tBeschriftungen;
tLage = (lLinks,lOben,lRechts,lUnten);
- tFFTDatenordnung = (doResIms,doResSmi,doRes,doBetr,doBetrQdr);
tFenster = object
aktiv: boolean;
Breite,Rand: longint;
diff --git a/werteunit.pas b/werteunit.pas
index aea3462..115345b 100644
--- a/werteunit.pas
+++ b/werteunit.pas
@@ -5,7 +5,7 @@ unit werteunit;
interface
uses
- Classes, SysUtils, typenunit, math, process, lowlevelunit, matheunit;
+ Classes, SysUtils, typenunit, math, process, lowlevelunit, matheunit, fftunit;
type
// tLLWerte ********************************************************************
@@ -43,8 +43,8 @@ type
procedure kopiereVerzerrt(original: pTLLWerteExtended; ZPs: tIntPointArray; ZGs: tExtPointArray; ZAs: tExtendedArray; xmin,xmax,tmin,tmax: longint; vb,nb: tTransformation; va,na: longint); overload;
destructor destroy; override;
function liesDateien(dateien: tGenerischeInputDateiInfoArray): boolean;
- function fft(senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fen: tFenster; out pvFehler: extended): boolean; overload;
- function fft(xmin,xmax,tmin,tmax: longint; senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fen: tFenster; out pvFehler: extended): boolean; overload;
+ function fft(senkrecht,invers: boolean; const algo: tFFTAlgorithmus; const fen: tFenster; out pvFehler: extended): boolean; overload; inline;
+ function fft(smin,smax: longint; senkrecht,invers: boolean; const algo: tFFTAlgorithmus; const fen: tFenster; out pvFehler: extended): boolean; overload;
procedure spiegle; overload;
procedure spiegle(tmin,tmax: longint); overload;
procedure fft2dNachbearbeitungA(nb: tFFTDatenordnung);
@@ -636,329 +636,84 @@ begin
end;
end;
-function tLLWerte.fft(senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fen: tFenster; out pvFehler: extended): boolean;
-var
- len: longint;
+function tLLWerte.fft(senkrecht,invers: boolean; const algo: tFFTAlgorithmus; const fen: tFenster; out pvFehler: extended): boolean;
begin
- len:=1;
- if senkrecht then begin
- while 2*len<=params.tsiz do
- len:=len*2;
- result:=fft(0,params.xsteps-1,(params.tsiz-len) div 2,((params.tsiz+len) div 2) - 1,true,invers,vor,nach,fen,pvFehler);
- end
- else begin
- while 2*len<=params.xsteps do
- len:=len*2;
- result:=fft((params.xsteps-len) div 2,((params.xsteps+len) div 2) - 1,0,params.tsiz-1,false,invers,vor,nach,fen,pvFehler);
- end;
+ if senkrecht then
+ result:=fft(0,params.xsteps-1,senkrecht,invers,algo,fen,pvFehler)
+ else
+ result:=fft(0,params.tsiz-1,senkrecht,invers,algo,fen,pvFehler);
end;
-function tLLWerte.fft(xmin,xmax,tmin,tmax: longint; senkrecht,invers: boolean; const vor,nach: tFFTDatenordnung; const fen: tFenster; out pvFehler: extended): boolean;
+function tLLWerte.fft(smin,smax: longint; senkrecht,invers: boolean; const algo: tFFTAlgorithmus; const fen: tFenster; out pvFehler: extended): boolean;
var
- i,j,k,n,dist,absch,wnum,wstep,haL,
- pmax,pmin,smax,smin: longint;
- in0,out0: boolean;
- ims,wRe,wIm: tExtendedArray;
- t1,t2,vorher,nachher,fenavg: extended;
- umsortierung: tLongintArray;
-const faktoren: array[tFFTDatenordnung,0..2] of longint = // (doResIms,doResSmi,doRes,doBetr,doBetrQdr);
- ((2,1,1),(2,1,1),(1,1,1),(1,1,1),(1,1,1));
+ i,j,pmax,pstep,sstep: longint;
+ in0,out0: boolean;
+ vorher,nachher,fenavg: extended;
begin
result:=false;
if senkrecht then begin
- pmax:=tmax;
- pmin:=tmin;
- smax:=xmax;
- smin:=xmin;
+ pmax:=params.tsiz-1;
+ pstep:=params.xsteps;
+ sstep:=1;
end
else begin
- pmax:=xmax;
- pmin:=xmin;
- smax:=tmax;
- smin:=tmin;
- end;
- if length(fen.Werte)<>pmax+1-pmin then
- fen.berechneWerte(pmax+1-pmin);
- n:=round(ln(pmax+1-pmin)/ln(2));
- if round(power(2,n))<>pmax+1-pmin then begin
- gibAus(inttostr(pmax+1-pmin)+' ist keine Zweierpotenz, ich bin zu dumm um Werte dieser Anzahl fourierzutransformieren!',3);
- exit;
- faktoren[doRes,0]:=faktoren[doRes,1]; // um so eine dämliche Warnung loszuwerden ...
+ pmax:=params.xsteps-1;
+ pstep:=1;
+ sstep:=params.xsteps;
end;
+ if length(fen.Werte)<>pmax+1 then
+ fen.berechneWerte(pmax+1);
- haL:=(pmax+1-pmin) div 2;
if fen.aktiv then begin
- if senkrecht then begin
- for i:=xmin to xmax do // Werte fenstern
- for j:=0 to tmax-tmin do
- werte[i+(j+tmin)*params.xsteps]:=werte[i+(j+tmin)*params.xsteps]*fen.Werte[j];
- end
- else begin
- for i:=0 to xmax-xmin do
- for j:=tmin to tmax do // Werte fenstern
- werte[i+xmin+j*params.xsteps]:=werte[i+xmin+j*params.xsteps]*fen.Werte[i];
- end;
+ for i:=smin to smax do // Werte fenstern
+ for j:=0 to pmax do
+ werte[i*sstep+j*pstep]:=
+ werte[i*sstep+j*pstep]*fen.Werte[j];
fenavg:=0;
- for i:=0 to pmax-pmin do
+ for i:=0 to pmax do
fenavg:=fenavg+fen.Werte[i];
- fenavg:=fenavg/(pmax-pmin+1);
+ fenavg:=fenavg/(pmax+1);
if fenavg<0.5 then begin
- gibAus('Sehr geringer Fensterdurchschnitt: '+floattostr(fenavg)+' ('+inttostr(fen.Breite)+' '+inttostr(length(fen.werte))+' '+inttostr(pmax-pmin+1)+')!',3);
+ gibAus('Sehr geringer Fensterdurchschnitt: '+floattostr(fenavg)+' ('+inttostr(fen.Breite)+' '+inttostr(length(fen.werte))+' '+inttostr(pmax+1)+')!',3);
end;
end;
vorher:=0;
- if vor=doBetrQdr then begin
- for i:=xmin to xmax do
- for j:=tmin to tmax do
- vorher:=vorher + werte[i+j*params.xsteps];
- end
- else begin
- for i:=xmin to xmax do
- for j:=tmin to tmax do
- vorher:=vorher + werte[i+j*params.xsteps]*werte[i+j*params.xsteps]
- *faktoren[vor,
- Byte(((not senkrecht) and (i=xmin)) or (senkrecht and (j=tmin)))
- +2*Byte(((not senkrecht) and (2*i=xmin+xmax+1)) or (senkrecht and (2*j=tmin+tmax+1)))];
- end;
-
- setlength(umsortierung,pmax+1-pmin);
- for j:=0 to pmax-pmin do begin
- k:=0;
- for i:=0 to n-1 do // Bits umkehren
- k:=2*k + byte(odd(j shr i));
- umsortierung[j]:=k;
- end;
- setlength(ims,pmax+1-pmin);
- setlength(wRe,haL);
- setlength(wIm,haL);
- for i:=0 to haL-1 do begin
- wRe[i]:=cos(pi*2*i/(pmax+1-pmin));
- wIm[i]:=sin(pi*2*i/(pmax+1-pmin));
- end;
-
+ nachher:=0;
in0:=true;
out0:=true;
- for i:=smin to smax do begin
- k:=1-2*byte(invers); // Eingang komplex konjugieren
- case vor of
- doResIms: begin // Re_0 Re_1=Re_n-1, ... Re_n/2-1=Re_n/2+1 Re_n/2 Im_1=-Im_n-1, ... Im_n/2-1=-Im_n/2+1
- if pmin<>0 then begin // pmin != 0 wird bei dieser Eingangscodierung nicht geduldet (warum auch?)
- gibAus('Bei getrennten Imaginär- und Realteilen wird eine Verschiebung der Eingangsdaten in Transformationsrichtung nicht geduldet - wie kommt soetwas denn zustande?',3);
- exit;
- end;
- if senkrecht then begin
- for j:=1 to haL-1 do begin
- Ims[j]:=k*werte[i+(haL+j)*params.xsteps];
- Ims[2*haL-j]:=-k*werte[i+(haL+j)*params.xsteps];
- werte[i+(haL+j)*params.xsteps]:=werte[i+(haL-j)*params.xsteps];
- end;
- end
- else begin
- for j:=1 to haL-1 do begin
- Ims[j]:=k*werte[haL+j+i*params.xsteps];
- Ims[2*haL-j]:=-k*werte[haL+j+i*params.xsteps];
- werte[haL+j+i*params.xsteps]:=werte[haL-j+i*params.xsteps];
- end;
- end;
- ims[0]:=0;
- ims[haL]:=0;
- end;
- doResSmi: begin // Re_0 Re_1=Re_n-1, ... Re_n/2-1=Re_n/2+1 Re_n/2 -Im_n/2+1=Im_n/2-1, ... -Im_n-1=Im_1
- if pmin<>0 then begin // pmin != 0 wird bei dieser Eingangscodierung nicht geduldet (warum auch?)
- gibAus('Bei getrennten Imaginär- und Realteilen wird eine Verschiebung der Eingangsdaten in Transformationsrichtung nicht geduldet - wie kommt soetwas denn zustande?',3);
- exit;
- end;
- if senkrecht then begin
- for j:=1 to haL-1 do begin
- Ims[j]:=k*werte[i+(2*haL-j)*params.xsteps];
- Ims[2*haL-j]:=-Ims[j];
- werte[i+(2*haL-j)*params.xsteps]:=werte[i+(j)*params.xsteps];
- end;
- end
- else begin
- for j:=1 to haL-1 do begin
- Ims[j]:=k*werte[2*haL-j+i*params.xsteps];
- Ims[2*haL-j]:=-Ims[j];
- werte[2*haL-j+i*params.xsteps]:=werte[j+i*params.xsteps];
- end;
- end;
- ims[0]:=0;
- ims[haL]:=0;
- end;
- doRes: // im = 0
- for j:=0 to pmax-pmin do
- ims[j]:=0;
- doBetr,doBetrQdr: begin
- gibAus('Ich brauche mehr als Beträge oder Betragsquadrate um eine FFT auzurechnen!',3);
- exit;
- end;
- end{of case};
- if senkrecht then begin
- for j:=0 to tmax-tmin do begin
- if umsortierung[j]>j then begin
- t1:=werte[i+(j+tmin)*params.xsteps];
- werte[i+(j+tmin)*params.xsteps]:=werte[i+(umsortierung[j]+tmin)*params.xsteps];
- werte[i+(umsortierung[j]+tmin)*params.xsteps]:=t1;
- t1:=ims[j];
- ims[j]:=ims[umsortierung[j]];
- ims[umsortierung[j]]:=t1;
- end;
- in0:=in0 and (werte[i+(j+tmin)*params.xsteps]=0) and (ims[j]=0);
- end;
- end
- else begin
- for j:=0 to xmax-xmin do begin
- if umsortierung[j]>j then begin
- t1:=werte[j+xmin+i*params.xsteps];
- werte[j+xmin+i*params.xsteps]:=werte[umsortierung[j]+xmin+i*params.xsteps];
- werte[umsortierung[j]+xmin+i*params.xsteps]:=t1;
- t1:=ims[j];
- ims[j]:=ims[umsortierung[j]];
- ims[umsortierung[j]]:=t1;
- end;
- in0:=in0 and (werte[j+xmin+i*params.xsteps]=0) and (ims[j]=0);
- end;
- end;
- dist:=1;
- wstep:=haL;
- if senkrecht then begin
- while dist<tmax+1-tmin do begin
- absch:=0;
- while absch<tmax+1-tmin do begin
- wnum:=0;
- for j:=0 to dist-1 do begin
- // x_links: [absch+j]
- // x_rechts: [absch+j+dist]
- t1:=wRe[wnum]*werte[i+(absch+j+tmin+dist)*params.xsteps] - wIm[wnum]*ims[absch+j+dist];
- t2:=wRe[wnum]*ims[absch+j+dist] + wIm[wnum]*werte[i+(absch+j+tmin+dist)*params.xsteps];
- werte[i+(absch+j+tmin+dist)*params.xsteps]:=werte[i+(absch+j+tmin)*params.xsteps]-t1;
- ims[absch+j+dist]:=ims[absch+j]-t2;
- werte[i+(absch+j+tmin)*params.xsteps]:=werte[i+(absch+j+tmin)*params.xsteps]+t1;
- ims[absch+j]:=ims[absch+j]+t2;
- wnum:=wnum+wstep;
- end;
- absch:=absch+2*dist;
- end;
- wstep:=wstep div 2;
- dist:=dist*2;
+ case sizeof(wgen) of
+ sizeof(single):
+ for i:=smin to smax do begin
+ algo.laden(invers,pSingle(@(werte[i*sstep])),pstep);
+ algo.summen(vorher,in0);
+ algo.ausfuehren;
+ algo.summen(nachher,out0);
+ algo.speichern(invers,pSingle(@(werte[i*sstep])),pstep);
end;
- end
- else begin
- while dist<xmax+1-xmin do begin
- absch:=0;
- while absch<xmax+1-xmin do begin
- wnum:=0;
- for j:=0 to dist-1 do begin
- // x_links: [absch+j]
- // x_rechts: [absch+j+dist]
- t1:=wRe[wnum]*werte[absch+j+xmin+dist+i*params.xsteps] - wIm[wnum]*ims[absch+j+dist];
- t2:=wRe[wnum]*ims[absch+j+dist] + wIm[wnum]*werte[absch+j+xmin+dist+i*params.xsteps];
-
- werte[absch+j+xmin+dist+i*params.xsteps]:=werte[absch+j+xmin+i*params.xsteps]-t1;
- ims[absch+j+dist]:=ims[absch+j]-t2;
- werte[absch+j+xmin+i*params.xsteps]:=werte[absch+j+xmin+i*params.xsteps]+t1;
- ims[absch+j]:=ims[absch+j]+t2;
- wnum:=wnum+wstep;
- end;
- absch:=absch+2*dist;
- end;
- wstep:=wstep div 2;
- dist:=dist*2;
+ sizeof(double):
+ for i:=smin to smax do begin
+ algo.laden(invers,pDouble(@(werte[i*sstep])),pstep);
+ algo.summen(vorher,in0);
+ algo.ausfuehren;
+ algo.summen(nachher,out0);
+ algo.speichern(invers,pDouble(@(werte[i*sstep])),pstep);
end;
- end;
- if senkrecht then begin
- case nach of // werte[ i + {tmin..tmax}*xsteps] & ims[{0..tmax-tmin}] -> werte[ i + {0..tmax-tmin}*xsteps ]
- doResIms: begin // Re_0 Re_1 Re_2 ... Re_n/2 Im_1 Im_2 ... Im_n/2-1
- for j:=0 to haL do // Realteile zusammensuchen
- werte[i+j*params.xsteps]:=werte[i+(j+tmin)*params.xsteps]/sqrt(tmax+1-tmin);
- for j:=1 to haL - 1 do // Imaginärteile zusammensuchen
- werte[i+(j+haL)*params.xsteps]:=ims[j]/sqrt(tmax+1-tmin);
- end;
- doResSmi: begin // Re_0 Re_1 Re_2 ... Re_n/2 Im_n/2-1 Im_n/2-2 ... Im_1
- for j:=0 to haL do // Realteile zusammensuchen
- werte[i+j*params.xsteps]:=werte[i+(j+tmin)*params.xsteps]/sqrt(tmax+1-tmin);
- for j:=1 to haL - 1 do // Imaginärteile zusammensuchen
- werte[i+(2*haL-j)*params.xsteps]:=ims[j]/sqrt(tmax+1-tmin);
- end;
- doBetr: // Abs(Re + i*Im) = Sqrt(Re^2 + Im^2)
- for j:=0 to tmax-tmin do
- werte[i+j*params.xsteps]:=sqrt((werte[i+(j+tmin)*params.xsteps]*werte[i+(j+tmin)*params.xsteps]+sqr(ims[j]))/(tmax+1-tmin));
- doBetrQdr: // Re^2 + Im^2
- for j:=0 to tmax-tmin do
- werte[i+j*params.xsteps]:=(werte[i+(j+tmin)*params.xsteps]*werte[i+(j+tmin)*params.xsteps]+sqr(ims[j]))/(tmax+1-tmin);
- doRes: // Re
- for j:=0 to tmax-tmin do
- werte[i+j*params.xsteps]:=werte[i+(j+tmin)*params.xsteps]/sqrt(tmax+1-tmin);
- else begin
- gibAus('Flieht, ihr Narren! (diese Meldung sollte eigentlich nicht auftreten können)',3);
- exit;
- end;
- end;
- for j:=0 to tmax-tmin do
- out0:=out0 and (werte[i+j*params.xsteps]=0);
- end
- else begin
- case nach of // werte[ {xtmin..xmax}+i*xsteps] & ims[{0..xmax-xmin}] -> werte[ {0..xmax-xmin}+i*xsteps ]
- doResIms: begin // Re_0 Re_1 Re_2 ... Re_n/2 Im_1 Im_2 ... Im_n/2-1
- for j:=0 to haL do // Realteile zusammensuchen
- werte[j+i*params.xsteps]:=werte[j+xmin+i*params.xsteps]/sqrt(xmax+1-xmin);
- for j:=1 to haL - 1 do // Imaginärteile zusammensuchen
- werte[j+haL+i*params.xsteps]:=ims[j]/sqrt(xmax+1-xmin);
- end;
- doResSmi: begin // Re_0 Re_1 Re_2 ... Re_n/2 Im_n/2-1 Im_n/2-2 ... Im_1
- for j:=0 to haL do // Realteile zusammensuchen
- werte[j+i*params.xsteps]:=werte[j+xmin+i*params.xsteps]/sqrt(xmax+1-xmin);
- for j:=1 to haL - 1 do // Imaginärteile zusammensuchen
- werte[2*haL-j+i*params.xsteps]:=ims[j]/sqrt(xmax+1-xmin);
- end;
- doBetr: // Abs(Re + i*Im) = Sqrt(Re^2 + Im^2)
- for j:=0 to xmax-xmin do
- werte[j+i*params.xsteps]:=sqrt((werte[j+xmin+i*params.xsteps]*werte[j+xmin+i*params.xsteps]+sqr(ims[j]))/(xmax+1-xmin));
- doBetrQdr: // Re^2 + Im^2
- for j:=0 to xmax-xmin do
- werte[j+i*params.xsteps]:=(werte[j+xmin+i*params.xsteps]*werte[j+xmin+i*params.xsteps]+sqr(ims[j]))/(xmax+1-xmin);
- doRes: // Re
- for j:=0 to xmax-xmin do
- werte[j+i*params.xsteps]:=werte[j+xmin+i*params.xsteps]/sqrt(xmax+1-xmin);
- else begin
- gibAus('Flieht, ihr Narren! (diese Meldung sollte eigentlich nicht auftreten können)',3);
- exit;
- end;
+ sizeof(extended):
+ for i:=smin to smax do begin
+ algo.laden(invers,pExtended(@(werte[i*sstep])),pstep);
+ algo.summen(vorher,in0);
+ algo.ausfuehren;
+ algo.summen(nachher,out0);
+ algo.speichern(invers,pExtended(@(werte[i*sstep])),pstep);
end;
- for j:=0 to xmax-xmin do
- out0:=out0 and (werte[j+i*params.xsteps]=0);
- end;
- end;
+ end{of case};
- gibAus(inttostr(xmin)+'-'+inttostr(xmax)+'x'+inttostr(tmin)+'-'+inttostr(tmax),1);
- if senkrecht then begin
- tmax:=tmax-tmin;
- tmin:=0;
- end
- else begin
- xmax:=xmax-xmin;
- xmin:=0;
- end;
- nachher:=0;
- if nach=doBetrQdr then begin
- for i:=xmin to xmax do
- for j:=tmin to tmax do
- nachher:=nachher + werte[i+j*params.xsteps];
- end
- else begin
- for i:=xmin to xmax do
- for j:=tmin to tmax do
- nachher:=nachher + werte[i+j*params.xsteps]*werte[i+j*params.xsteps]
- *faktoren[nach,
- Byte(((not senkrecht) and (i=xmin)) or (senkrecht and (j=tmin)))
- +2*Byte(((not senkrecht) and (2*i=xmin+xmax+1)) or (senkrecht and (2*j=tmin+tmax+1)))];
- end;
if (nachher=0) and (vorher=0) then pvFehler:=0
else pvFehler:=abs(nachher-vorher)/(nachher+vorher);
- gibAus(inttostr(byte(senkrecht))+' '+inttostr(byte(vor))+' '+inttostr(byte(nach))+' '+inttostr(byte(invers))+' (Parseval-Fehler = '+floattostr(pvFehler)+') ... nämlich von '+floattostr(vorher)+' zu '+floattostr(nachher),1);
+ gibAus(inttostr(byte(senkrecht))+' '+inttostr(byte(invers))+' (Parseval-Fehler = '+floattostr(pvFehler)+') ... nämlich von '+floattostr(vorher)+' zu '+floattostr(nachher),1);
if in0 then gibAus('Nur Nullen im Input der FFT!',1);
if out0 then gibAus('Nur Nullen im Output der FFT!',1);
result:=not (in0 or out0);
@@ -1479,10 +1234,11 @@ end;
function tWavelet.berechneWerte: boolean;
var
- i: longint;
- tmp: extended;
- fenster: tFenster;
- nur0: boolean;
+ i: longint;
+ tmp: extended;
+ fenster: tFenster;
+ nur0: boolean;
+ tmpFFTAlgo: tFFTAlgorithmus;
begin
result:=false;
werte.params.xsteps:=2;
@@ -1517,7 +1273,9 @@ begin
// fenster.Breite:=0;
// fenster.berechneWerte(werte.tsiz);
gibAus(inttostr(werte.params.xsteps)+' '+inttostr(werte.params.tsiz)+' '+inttostr(length(werte.werte)),1);
- werte.fft(0,1,0,werte.params.tsiz-1,true,false,doRes,doResIms,fenster,pvFehler);
+ tmpFFTAlgo:=createFFTAlgorithmus(2*hlen,doRes,doResIms);
+ werte.fft(true,false,tmpFFTAlgo,fenster,pvFehler);
+ tmpFFTAlgo.free;
end;
wtFrequenzfenster: begin
hlen:=werte.params.tsiz div 2;
generated by cgit v1.2.3-70-g09d2 (git 2.47.0) at 2024-11-02 10:40:23 +0000