diff options
| -rw-r--r-- | fftunit.inc | 199 | ||||
| -rw-r--r-- | fftunit.pas | 326 | ||||
| -rw-r--r-- | matheunit.pas | 73 |
3 files changed, 525 insertions, 73 deletions
diff --git a/fftunit.inc b/fftunit.inc new file mode 100644 index 0000000..7f5cc34 --- /dev/null +++ b/fftunit.inc @@ -0,0 +1,199 @@ +{$IFDEF fftunitLaden} +//procedure tFFTAlgorithmus.laden(invers: boolean; const q: single; schritt: longint); +var + i,inv: longint; + {$IFDEF fftNormierung} + lenFak: extended; + {$ENDIF} +begin + inv:=1-2*byte(invers); + {$IFDEF fftNormierung} + lenFak:=1/sqrt(len); + {$ENDIF} + case ino of + doRes: + for i:=0 to len-1 do begin + {$IFDEF tFFTCooleyTukey} + res[perm[i]]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF}; + ims[perm[i]]:=0; + {$ELSE} + res[i]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF}; + ims[i]:=0; + {$ENDIF} + end; + doAlleResIms: + for i:=0 to len-1 do begin + {$IFDEF tFFTCooleyTukey} + res[perm[i]]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF}; + ims[perm[i]]:=(q+(i+len)*schritt)^ * inv {$IFDEF fftNormierung} * lenFak{$ENDIF}; + {$ELSE} + res[i]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF}; + ims[i]:=(q+(i+len)*schritt)^ * inv {$IFDEF fftNormierung} * lenFak{$ENDIF}; + {$ENDIF} + end; + doAlleResSmi: + for i:=0 to hlen + (len and $1) - 1 do begin + {$IFDEF tFFTCooleyTukey} + res[perm[i]]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF}; + ims[perm[i]]:=(q+(2*len-1-i)*schritt)^ * inv {$IFDEF fftNormierung} * lenFak{$ENDIF}; + {$ELSE} + res[i]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF}; + ims[i]:=(q+(2*len-1-i)*schritt)^ * inv {$IFDEF fftNormierung} * lenFak{$ENDIF}; + {$ENDIF} + end; + doResIms: begin + res[0]:=q^ {$IFDEF fftNormierung} * lenFak{$ENDIF}; // perm[0]=0 + ims[0]:=0; + for i:=1 to hlen + (len and $1) - 1 do begin + {$IFDEF tFFTCooleyTukey} + res[perm[i]]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF}; + ims[perm[i]]:=(q+(i+hlen)*schritt)^ * inv {$IFDEF fftNormierung} * lenFak{$ENDIF}; + res[perm[len-i]]:=res[perm[i]]; + ims[perm[len-i]]:=-ims[perm[i]]; + {$ELSE} + res[i]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF}; + ims[i]:=(q+(i+hlen)*schritt)^ * inv {$IFDEF fftNormierung} * lenFak{$ENDIF}; + res[len-i]:=res[i]; + ims[len-i]:=-ims[i]; + {$ENDIF} + end; + if not odd(len) then begin + {$IFDEF tFFTCooleyTukey} + res[1]:=(q+hlen*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF}; // perm[hlen]=1 + ims[1]:=0; + {$ELSE} + res[hlen]:=(q+hlen*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF}; + ims[hlen]:=0; + {$ENDIF} + end; + end; + doResSmi: begin + res[0]:=q^ {$IFDEF fftNormierung} * lenFak{$ENDIF}; // perm[0]=0 + ims[0]:=0; + for i:=1 to hlen + (len and $1) - 1 do begin + {$IFDEF tFFTCooleyTukey} + res[perm[i]]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF}; + ims[perm[i]]:=(q+(len-i)*schritt)^ * inv {$IFDEF fftNormierung} * lenFak{$ENDIF}; + res[perm[len-i]]:=res[perm[i]]; + ims[perm[len-i]]:=-ims[perm[i]]; + {$ELSE} + res[i]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF}; + ims[i]:=(q+(len-i)*schritt)^ * inv {$IFDEF fftNormierung} * lenFak{$ENDIF}; + res[len-i]:=res[i]; + ims[len-i]:=-ims[i]; + {$ENDIF} + end; + if not odd(len) then begin + {$IFDEF tFFTCooleyTukey} + res[1]:=(q+hlen*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF}; // perm[hlen]=1 + ims[1]:=0; + {$ELSE} + res[hlen]:=(q+hlen*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF}; + ims[hlen]:=0; + {$ENDIF} + end; + end; + doBetr,doBetrQdr: + fehler('Ich brauche mehr als Beträge oder Betragsquadrate um eine FFT auszurechnen!'); + else + fehler('Diese Inputordnung ist in dieser Methode nicht implementiert!');; + end{of case}; +end; +{$ENDIF} + +{$IFDEF fftunitGetrenntLaden} +//procedure tFFTAlgorithmus.laden(invers: boolean; qRe,qIm: tExtendedArray); +var + i,inv: longint; + {$IFDEF fftNormierung} + lenFak: extended; + {$ENDIF} +begin + inv:=1-2*byte(invers); + {$IFDEF fftNormierung} + lenFak:=1/sqrt(len); + {$ENDIF} + case ino of + doGetrennt: + for i:=0 to len-1 do begin + {$IFDEF tFFTCooleyTukey} + res[perm[i]]:=qRe[i] {$IFDEF fftNormierung} * lenFak{$ENDIF}; + ims[perm[i]]:=qIm[i] * inv {$IFDEF fftNormierung} * lenFak{$ENDIF}; + {$ELSE} + res[i]:=qRe[i] {$IFDEF fftNormierung} * lenFak{$ENDIF}; + ims[i]:=qIm[i] * inv {$IFDEF fftNormierung} * lenFak{$ENDIF}; + {$ENDIF} + end; + else + fehler('Diese Inputordnung ist in dieser Methode nicht implementiert!');; + end{of case}; +end; +{$ENDIF} + +{$IFDEF fftunitSpeichern} +//procedure tFFTAlgorithmus.speichern(invers: boolean; var z: single; schritt: longint); +var + i,inv: longint; +begin + inv:=1-2*byte(invers); + case outo of + doRes: + for i:=0 to len-1 do + (z+i*schritt)^:=res[i]; + doAlleResIms: + for i:=0 to len-1 do begin + (z+i*schritt)^:=res[i]; + (z+(i+len)*schritt)^:=ims[i] * inv; + end; + doAlleResSmi: + for i:=0 to len-1 do begin + (z+i*schritt)^:=res[i]; + (z+(2*len-1-i)*schritt)^:=ims[i] * inv; + end; + doResIms: begin + z^:=res[0]; + for i:=1 to hlen + (len and $1) - 1 do begin + (z+i*schritt)^:=res[i]; + (z+(i+hlen)*schritt)^:=ims[i] * inv; + end; + if not odd(len) then + (z+hlen*schritt)^:=res[hlen]; + end; + doResSmi: begin + z^:=res[0]; + for i:=1 to hlen + (len and $1) - 1 do begin + (z+i*schritt)^:=res[i]; + (z+(len-i)*schritt)^:=ims[i] * inv; + end; + if not odd(len) then + (z+hlen*schritt)^:=res[hlen]; + end; + doBetr: + for i:=0 to len-1 do + (z+i*schritt)^:=sqrt(sqr(res[i])+sqr(ims[i])); + doBetrQdr: + for i:=0 to len-1 do + (z+i*schritt)^:=sqr(res[i])+sqr(ims[i]); + else + fehler('Diese Outputordnung ist in dieser Methode nicht implementiert!');; + end{of case}; +end; +{$ENDIF} + +{$IFDEF fftunitGetrenntSpeichern} +//procedure tFFTAlgorithmus.speichern(invers: boolean; zRe,zIm: tExtendedArray); +var + i,inv: longint; +begin + inv:=1-2*byte(invers); + case outo of + doGetrennt: + for i:=0 to len-1 do begin + zRe[i]:=res[i]; + zIm[i]:=ims[i] * inv; + end; + else + fehler('Diese Outputordnung ist in dieser Methode nicht implementiert!');; + end{of case}; +end; +{$ENDIF} diff --git a/fftunit.pas b/fftunit.pas new file mode 100644 index 0000000..f78935f --- /dev/null +++ b/fftunit.pas @@ -0,0 +1,326 @@ +unit fftunit; + +{$mode objfpc}{$H+} + +interface + +uses + Classes, SysUtils, lowlevelunit; + +{$DEFINE fftNormierung} + +type + tFFTDatenordnung = (doResIms,doResSmi,doAlleResIms,doAlleResSmi,doRes,doBetr,doBetrQdr,doGetrennt); + + tFFTAlgorithmus = class + private + len,hlen: longint; + res,ims: tExtendedArray; + ino,outo: tFFTDatenordnung; + public + procedure dumpResIms; + constructor create(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung); + destructor destroy; override; + procedure laden(invers: boolean; q: pSingle; schritt: longint); overload; virtual; + procedure laden(invers: boolean; q: pDouble; schritt: longint); overload; virtual; + procedure laden(invers: boolean; q: pExtended; schritt: longint); overload; virtual; + procedure laden(invers: boolean; qRe,qIm: tExtendedArray); overload; virtual; + procedure ausfuehren; virtual; abstract; + procedure speichern(invers: boolean; z: pSingle; schritt: longint); overload; virtual; + procedure speichern(invers: boolean; z: pDouble; schritt: longint); overload; virtual; + procedure speichern(invers: boolean; z: pExtended; schritt: longint); overload; virtual; + procedure speichern(invers: boolean; zRe,zIm: tExtendedArray); overload; virtual; + end; + + tFFTCooleyTukey = class(tFFTAlgorithmus) + private + wre,wim: tExtendedArray; + perm: tLongintArray; + public + constructor create(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung); + destructor destroy; override; + procedure laden(invers: boolean; q: pSingle; schritt: longint); overload; override; + procedure laden(invers: boolean; q: pDouble; schritt: longint); overload; override; + procedure laden(invers: boolean; q: pExtended; schritt: longint); overload; override; + procedure laden(invers: boolean; qRe,qIm: tExtendedArray); overload; override; + procedure ausfuehren; override; + end; + + tFFTBluestein = class(tFFTAlgorithmus) + private + rRe,rIm,sRe,sIm: array of extended; + subFFT: tFFTCooleyTukey; + public + constructor create(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung); + destructor destroy; override; + procedure ausfuehren; override; + end; + + function createFFTAlgorithmus(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung): tFFTAlgorithmus; + +implementation + +uses + math; + +// allgemeine Funktionen ******************************************************* + +function createFFTAlgorithmus(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung): tFFTAlgorithmus; +begin + if inputOrdnung in [doBetr,doBetrQdr] then + fehler('Ich brauche mehr als Beträge oder Betragsquadrate um eine FFT auszurechnen!'); + + if round(power(2,round(ln(laenge)/ln(2))))=laenge then + result:=tFFTCooleyTukey.create(laenge,inputOrdnung,outputOrdnung) + else + result:=tFFTBluestein.create(laenge,inputOrdnung,outputOrdnung); +end; + +// tFFTAlgorithmus ************************************************************* +{$DEFINE tFFTAlgorithmus} + +constructor tFFTAlgorithmus.create(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung); +begin + inherited create; + ino:=inputOrdnung; + outo:=outputOrdnung; + len:=laenge; + hlen:=laenge div 2; + fillchar(res,sizeof(res),#0); + setlength(res,0); + fillchar(ims,sizeof(ims),#0); + setlength(ims,0); +end; + +destructor tFFTAlgorithmus.destroy; +begin + setlength(res,0); + setlength(ims,0); + inherited destroy; +end; + +procedure tFFTAlgorithmus.dumpResIms; +var + i: longint; +begin + writeln('***'); + for i:=0 to length(res)-1 do + writeln(res[i],' ',ims[i]); + writeln('***'); +end; + +{$DEFINE fftunitLaden} +procedure tFFTAlgorithmus.laden(invers: boolean; q: pSingle; schritt: longint); +{$INCLUDE fftunit.inc} +procedure tFFTAlgorithmus.laden(invers: boolean; q: pDouble; schritt: longint); +{$INCLUDE fftunit.inc} +procedure tFFTAlgorithmus.laden(invers: boolean; q: pExtended; schritt: longint); +{$INCLUDE fftunit.inc} +{$UNDEF fftunitLaden} + +{$DEFINE fftunitGetrenntLaden} +procedure tFFTAlgorithmus.laden(invers: boolean; qRe,qIm: tExtendedArray); +{$INCLUDE fftunit.inc} +{$UNDEF fftunitGetrenntLaden} + +{$DEFINE fftunitSpeichern} +procedure tFFTAlgorithmus.speichern(invers: boolean; z: pSingle; schritt: longint); +{$INCLUDE fftunit.inc} +procedure tFFTAlgorithmus.speichern(invers: boolean; z: pDouble; schritt: longint); +{$INCLUDE fftunit.inc} +procedure tFFTAlgorithmus.speichern(invers: boolean; z: pExtended; schritt: longint); +{$INCLUDE fftunit.inc} +{$UNDEF fftunitSpeichern} + +{$DEFINE fftunitGetrenntSpeichern} +procedure tFFTAlgorithmus.speichern(invers: boolean; zRe,zIm: tExtendedArray); +{$INCLUDE fftunit.inc} +{$UNDEF fftunitGetrenntSpeichern} + +{$UNDEF tFFTAlgorithmus} + +// tFFTCooleyTukey ************************************************************* +{$DEFINE tFFTCooleyTukey} + +constructor tFFTCooleyTukey.create(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung); +var + i,j,loglen: longint; +begin + inherited create(laenge,inputOrdnung,outputOrdnung); + + setlength(res,len); + setlength(ims,len); + + loglen:=round(ln(len)/ln(2)); + setlength(wre,hlen); + setlength(wim,hlen); + setlength(perm,len); + for i:=0 to len-1 do begin + res[i]:=nan; + ims[i]:=nan; + perm[i]:=0; + for j:=0 to loglen-1 do + perm[i]:=perm[i] or (byte(odd(i shr j)) shl (loglen-1-j)); + end; + for i:=0 to hlen-1 do begin + wre[i]:=cos(2*pi*i/len); + wim[i]:=sin(2*pi*i/len); + end; +end; + +destructor tFFTCooleyTukey.destroy; +begin + setlength(wre,0); + setlength(wim,0); + setlength(perm,0); + inherited destroy; +end; + +{$DEFINE fftunitLaden} +procedure tFFTCooleyTukey.laden(invers: boolean; q: pSingle; schritt: longint); +{$INCLUDE fftunit.inc} +procedure tFFTCooleyTukey.laden(invers: boolean; q: pDouble; schritt: longint); +{$INCLUDE fftunit.inc} +procedure tFFTCooleyTukey.laden(invers: boolean; q: pExtended; schritt: longint); +{$INCLUDE fftunit.inc} +{$UNDEF fftunitLaden} + +{$DEFINE fftunitGetrenntLaden} +procedure tFFTCooleyTukey.laden(invers: boolean; qRe,qIm: tExtendedArray); +{$INCLUDE fftunit.inc} +{$UNDEF fftunitGetrenntLaden} + +procedure tFFTCooleyTukey.ausfuehren; +var + i,dist,absch,wnum,wstep: longint; + t1,t2: extended; +begin + dist:=1; + wstep:=hlen; + while dist<len do begin + absch:=0; + while absch<len do begin + wnum:=0; + for i:=0 to dist-1 do begin + // x_links: [absch+j] + // x_rechts: [absch+j+dist] + t1:=wRe[wnum]*res[absch+i+dist] - wIm[wnum]*ims[absch+i+dist]; + t2:=wRe[wnum]*ims[absch+i+dist] + wIm[wnum]*res[absch+i+dist]; + + res[absch+i+dist]:=res[absch+i]-t1; + ims[absch+i+dist]:=ims[absch+i]-t2; + res[absch+i]:=res[absch+i]+t1; + ims[absch+i]:=ims[absch+i]+t2; + wnum:=wnum+wstep; + end; + absch:=absch+2*dist; + end; + wstep:=wstep div 2; + dist:=dist*2; + end; +end; + +{$UNDEF tFFTCooleyTukey} + +// tFFTBluestein *************************************************************** + +constructor tFFTBluestein.create(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung); +var + preFFT: tFFTCooleyTukey; + e: extended; + i: longint; +begin + inherited create(laenge,inputOrdnung,outputOrdnung); + setlength(res,round(power(2,ceil(ln(2*len-1)/ln(2))))); + setlength(ims,length(res)); + setlength(rRe,length(res)); + setlength(rIm,length(res)); + setlength(sRe,length(res)); + setlength(sIm,length(res)); + subFFT:=tFFTCooleyTukey.create(length(res),doGetrennt,doGetrennt); + + for i:=0 to subFFT.len-1 do begin + e:=sqr(min(i,subFFT.len-i))*pi/len; + rRe[i]:=Cos(e); + rIm[i]:=Sin(e); + end; + + preFFT:=tFFTCooleyTukey.create(length(rRe),doGetrennt,doGetrennt); + preFFT.laden(false,rRe,rIm); + preFFT.ausfuehren; + preFFT.speichern(false,sRe,sIm); + preFFT.free; + + {$IFDEF fftNormierung} + for i:=0 to subFFT.len-1 do begin + sRe[i]:=sRe[i]*sqrt(subFFT.len); + sIm[i]:=sIm[i]*sqrt(subFFT.len); + end; + {$ENDIF} +end; + +destructor tFFTBluestein.destroy; +begin + setlength(rRe,0); + setlength(rIm,0); + setlength(sRe,0); + setlength(sIm,0); + subFFT.free; + res:=nil; + ims:=nil; + inherited destroy; +end; + +procedure tFFTBluestein.ausfuehren; +var + i: longint; + tmp: extended; +begin + // a[j] := werte[j] * CC(raumFaktoren[j]) + for i:=0 to len-1 do begin + tmp:=res[i]; + res[i]:= tmp*rRe[i] + ims[i]*rIm[i]; + ims[i]:=- tmp*rIm[i] + ims[i]*rRe[i]; + end; + for i:=len to length(res)-1 do begin + res[i]:=0; + ims[i]:=0; + end; + + // a[j] := iFFT(a[j]) + subFFT.laden(true,res,ims); + subFFT.ausfuehren; + subFFT.speichern(true,res,ims); + + for i:=0 to length(res)-1 do begin + res[i]:=res[i]; + ims[i]:=ims[i]; + end; + + // a[j] := a[j] * spektralFaktoren[j] + for i:=0 to length(res)-1 do begin + tmp:=res[i]; + res[i]:=tmp*sRe[i] - ims[i]*sIm[i]; + ims[i]:=tmp*sIm[i] + ims[i]*sRe[i]; + end; + + // a[j] := iFFT(a[j]) + subFFT.laden(false,res,ims); + subFFT.ausfuehren; + subFFT.speichern(false,res,ims); + + for i:=0 to length(res)-1 do begin + res[i]:=res[i]; + ims[i]:=ims[i]; + end; + + // werte[j] := a[j] * CC ( raumFaktoren[j] ) + for i:=0 to len-1 do begin + tmp:=res[i]; + res[i]:= tmp*rRe[i] + ims[i]*rIm[i]; + ims[i]:=- tmp*rIm[i] + ims[i]*rRe[i]; + end; +end; + +end. + diff --git a/matheunit.pas b/matheunit.pas index 431345c..61c1eea 100644 --- a/matheunit.pas +++ b/matheunit.pas @@ -49,7 +49,6 @@ function nullfunktion(x: extended): extended; function exprtofloat(st: boolean; s: string; kvs: tKnownValues; cbgv: tCallBackGetValue): extended; function formelnAuswerten(s: string; kvs: tKnownValues; cbgv: tCallBackGetValue): string; function knownValue(nam: string; val: extended): tKnownValue; -procedure fft(var res,ims: tExtendedArray; inv: boolean); inline; implementation @@ -593,77 +592,5 @@ begin result.value:=val; end; -procedure fft(var res,ims: tExtendedArray; inv: boolean); -var - haL,i,j,k,n,dist,wstep,absch,wnum: longint; - umsortierung: tLongintArray; - wRe,wIm: tExtendedArray; - t1,t2: extended; -begin - haL:=length(res) div 2; - n:=round(ln(length(res))/ln(2)); - - if round(power(2,n))<>length(res) then - raise exception.create(inttostr(length(res))+' ist keine Zweierpotenz!'); - - setlength(umsortierung,2*haL); - for j:=0 to 2*haL-1 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(wRe,haL); - setlength(wIm,haL); - for i:=0 to haL-1 do begin - wRe[i]:=cos(pi*i/haL); - wIm[i]:=sin(pi*i/haL); - end; - - for j:=0 to 2*haL-1 do begin - if umsortierung[j]>j then begin - t1:=res[j]; - res[j]:=res[umsortierung[j]]; - res[umsortierung[j]]:=t1; - t1:=ims[j]; - ims[j]:=ims[umsortierung[j]] * (1-2*byte(inv)); - ims[umsortierung[j]]:=t1 * (1-2*byte(inv)); - end; - if (umsortierung[j]=j) and inv then - ims[j]:=-ims[j]; - end; - - dist:=1; - wstep:=haL; - - while dist<2*haL do begin - absch:=0; - while absch<2*haL do begin - wnum:=0; - for j:=0 to dist-1 do begin - // x_links: [absch+j] - // x_rechts: [absch+j+dist] - t1:=wRe[wnum]*res[absch+j+dist] - wIm[wnum]*ims[absch+j+dist]; - t2:=wRe[wnum]*ims[absch+j+dist] + wIm[wnum]*res[absch+j+dist]; - - res[absch+j+dist]:=res[absch+j]-t1; - ims[absch+j+dist]:=ims[absch+j]-t2; - res[absch+j]:=res[absch+j]+t1; - ims[absch+j]:=ims[absch+j]+t2; - wnum:=wnum+wstep; - end; - absch:=absch+2*dist; - end; - wstep:=wstep div 2; - dist:=dist*2; - end; - - t1:=1/sqrt(length(res)); - for i:=0 to 2*haL-1 do begin - res[i]:=res[i]*t1; - ims[i]:=ims[i]*t1; - end; -end; - end. |