unit fft;

INTERFACE
USES
Math,sysUtils,common;
type Tmodul = procedure(var xr,yi:array of double; out virt_vect,angle:array of double;count:word;npdb,rms:boolean);// procedure(var xr,yi,virt_vect:rAryType;count:word;npdb,rms:boolean);
var Yzero,windowWeights,xr0,xr,yi,angle:TdbArr;
           jump,rectWinFill : realtype;
          checkBreak : boolean; //crt, reagowac natychmiast na ctrl-break
                  cp : word;    //z disco2
               modul : Tmodul;
          gaussStDev : double=maxDouble;
            rectFill : word=127;
         gaussScaler : word=high(word);
 gaussStDevUnscalled : longWord;
procedure Y0(var Yzero:array of double; count:word);
procedure modul_NpdBrms(var xr,yi:array of double; out virt_vect,angle:array of double;count:word;npdb,rms:boolean);
procedure filter_modul_NpdBrms(var xr,yi:array of double; out virt_vect,angle:array of double;count:word;npdb,rms:boolean);//filter_modul_NpdBrms(var xr,yi,virt_vect:rAryType;count:word;npdb,rms:boolean);
procedure WindowsWeights(const wn_nr:byte; const gaussStdev:single; out weights:array of double;out mean:extended;
                    const count:longint);
type Tindexa = array of word;
var index_mix,indexSort : Tindexa;
     maxpower : word;
procedure ind_mix_Arr(var index_mix:Tindexa;numdat:word);
procedure power2(var maxpower:word; numdat:longint);
var cosar,sinar : TdbArr;
procedure SinCosArrays(var cosary,sinary:array of double;numdat:word);
PROCEDURE rfft_calc(const index_mix:Tindexa;maxpower:word;
                var sinar,cosar:array of double; var xreal,yimag:array of double;numdat:word);
PROCEDURE fft_calc(var xreal:rAryType; var yimag:rAryType;numdat:integer);
procedure windowsWeightsSpectrDraw(const weights:array of double);

IMPLEMENTATION

uses newSpectrum,visualization,drawings,grafiti4,unit1,graphics;
procedure Y0(var Yzero:array of double; count:word);
 var i:word;
 begin
  for i:=0 to count-1 do Yzero[i]:=0;
 end;{Y0}

procedure modul_NpdBrms(var xr,yi:array of double; out virt_vect,angle:array of double;count:word;npdb,rms:boolean);
{
 procedura oblicza modul mocy w dB albo w neperach albo rms
 xr,yi - część rzeczywista i urojona, nie modyfikowane!
 virt_vect - wynik!
}
 var i:integer;
     t,t1,t2:extended;
 begin
  temp:=high(angle);   //debug
  for i:=0 to cp  do
   begin
    t1:=xr[i]; t2:=yi[i];
    try
     angle[i]:=arcTan(t2/t1);
    except
     angle[i]:=pi*sign(t1)/2
    end;
    t:=(sqr(t1)+sqr(t2))/sqrCount;    //kwadrat modułu
    if rms then virt_vect[i]:=sqrt(t)
    else
     begin
      if t<minDouble then
       t:=minDouble;         {ln(2)=0.693147181  obliczamy modul}
      if npdb then  virt_vect[i]:=0.5*ln(t)          { in nephers; 0.5, because ln(pierwiastek(t))!}
      else virt_vect[i]:=4.34294482*ln(t);           { in decibels, because 1dB=20*log(rms)=0.5*20*log(sqr(rms))= 10/ln(10)*ln(sqr(rms))=4.3429*ln(sqr(rms)); było  4.3429*ln(2)=3.010299385, bo mnożymy widmo przez 2 za wyjątkiem składowej zerowej, ale pomnożono 4*t}
     end;{else}
   end;{for}                                         { 10/ln(10)=4.3429...}
   for i:=cp+1 to count-1 do     //second half regain 06042013
    begin
     virt_vect[i]:=virt_vect[count-i];
     angle[i]:=-angle[count-i];
    end;
 end;{modul_NpdBrms}

procedure filter_modul_NpdBrms(var xr,yi:array of double; out virt_vect,angle:array of double;count:word;npdb,rms:boolean);
{
 procedura oblicza modul mocy w dB albo w neperach albo rms z uwzględnieniem funkcji filtrującej
 virt_vect - wynik!
}
 var i:integer;
     t,t1,t2:extended;
 begin
  for i:=0 to cp  do
   begin
    t1:=xr[i]; t2:=yi[i];
    try
     angle[i]:=arcTan(t2/t1);
    except
     angle[i]:=pi*sign(t1)/2
    end;
    t:=(sqr(t1)+sqr(t2))/sqrCount;       //kwadrat modułu
    if i<=cp then t:=t*filter[i]         //wprowadzić filtr
    else t:=t*filter[count-i];
    if rms then virt_vect[i]:=sqrt(t)
    else
     begin
      if t<1e-9 then t:=1e-9;    // if t<minSingle then t:=minSingle;       {ln(2)=0.693147181  obliczamy modul}// {ln(2)=0.693147181  obliczamy modul}
      if npdb then  virt_vect[i]:=0.5*ln(t)            { w neperach}
      else virt_vect[i]:=4.34294482*ln(t);             { w decybelach}
     end;{else}
   end;{for}                                           { 10/ln(10)=4.3429...}
   for i:=cp+1 to count-1 do   //second half regain
    begin
     virt_vect[i]:=virt_vect[count-i];
     angle[i]:=-angle[count-i];
    end;
 end;{filter_modul_NpdBrms}


 procedure windowsWeightsSpectrDraw(const weights:array of double);
  var i:word;
  begin
   for i:=0 to count-1 do begin xr[i]:=weights[i]; yi[i]:=0 end;
   rfft_calc(index_mix,maxpower,sinar,cosar,xr,yi,count);     //dividing by count is carried out in the "modul" procedure
   modul_NpdBrms(xr,yi,virt_vect,angle,count,form2.CheckBox3.Checked,true);    //true - obliczamy tu rms, aby móc włącza check dB; xr, yi - data, not modified here! virt_vect, angle - results form2.CheckBox3.Checked - true nepher, false - decibel
   form2.label18.Caption:='Windows weights spectrum';
   label18Caption:=form2.label18.Caption;
   filterDraw(virt_vect,form2,form2.panel6,form2.label18,clRed,wspx6,wspy6,centr6);
  end;{}

procedure WindowsWeights(const wn_nr:byte; const gaussStdev:single; out weights:array of double; out mean:extended;
                    const count:longint);
{
 Oblicza tabele wag okiennych;
}
 var zmp,zmp1,zmp3,min,max,stDev:extended;
     i:longint;
 begin //-------------------------WindowsWeights--------------------------------
   for i:=0 to count-1 do weights[i]:=yZero[i]; //200052013 jawnie z powodu niezgodności typów
  case wn_nr of

 0 :begin //Hamming
        zmp1:=2*pi/(count-1);
        i:=count div 2;
        max:=(0.53836-0.46164*cos(zmp1*i));
        for i:=0 to count-1 do
         begin
          zmp:=(0.53836-0.46164*cos(zmp1*i))/max;
          weights[i]:=zmp;
         end;{for}
       end;{Hamming}

 1 :begin  //gauss normalised, 1 standard dev.
     i:=count div 2;
     max:=exp(-0.5*sqr((2*i-count+1)/(gaussStdev*(count-1))));
     for i:=0 to count-1 do
      begin
       weights[i]:=exp(-0.5*sqr((2*i-count+1)/(gaussStdev*(count-1))))/max;
      end;{for}
    end;{gauss}

 2 :begin //Hann
        zmp1:=2*pi/(count-1);
        i:=count div 2;
        max:=0.5*(1-cos(zmp1*i));
        for i:=0 to count-1 do
         begin
          zmp:=0.5*(1-cos(zmp1*i))/max;
          weights[i]:=zmp;
         end;{for}
       end;{Hann}

 3 :begin //Flat Top
     zmp1:=2*pi/(count-1);
     i:=count div 2;
     max:=(1-1.933*cos(zmp1*i)+1.286*cos(2*zmp1*i)-0.388*cos(3*zmp1*i)+
              0.0322*cos(4*zmp1*i));
     for i:=0 to count-1 do
      begin
       zmp:= (1-1.933*cos(zmp1*i)+1.286*cos(2*zmp1*i)-0.388*cos(3*zmp1*i)+
              0.0322*cos(4*zmp1*i))/max;
       weights[i]:=zmp;
      end;{for}
    end;{Flat Top}

 4 :begin //Keiser-Bessel
     zmp1:=2*pi/(count-1);
     i:=count div 2;
     max:=1-1.24*cos(zmp1*i)+0.244*cos(2*zmp1*i)-0.00305*cos(3*zmp1*i);
     for i:=0 to count-1 do
      begin
       zmp:=((1-1.24*cos(zmp1*i)+0.244*cos(2*zmp1*i)-0.00305*cos(3*zmp1*i)))/max;
       weights[i]:=zmp;
      end;{for}
    end;{Keiser-Bessel}

 5 :begin  //Blackman-Harris
     zmp1:=2*pi/(count-1);
      i:=(count div 2);
     max:=0.35875-0.48829*cos(zmp1*i)+0.14128*cos(2*zmp1*i)-0.01168*cos(3*zmp1*i);
     for i:=0 to count-1 do
      begin
       zmp:=(0.35875-0.48829*cos(zmp1*i)+0.14128*cos(2*zmp1*i)-0.01168*cos(3*zmp1*i))/max;
       //zmp:= (1-1.36*cos(zmp1*i)+0.39*cos(2*zmp1*i)-0.032*cos(3*zmp1*i)); poprzednia wersja 01042013 (wyrazy podzielone przez 0.35875)
       weights[i]:=zmp;
      end;{for}
    end;{Blackman-Harris}

 6 :begin  //Blackman-Nutall
       zmp1:=2*pi/(count-1);
        i:=(count div 2);
       max:=(0.3635819-0.4891755*cos(zmp1*i)+0.1365995*cos(2*zmp1*i)-0.0106411*cos(3*zmp1*i));
       for i:=0 to count-1 do
        begin
         zmp:=(0.3635819-0.4891755*cos(zmp1*i)+0.1365995*cos(2*zmp1*i)-0.0106411*cos(3*zmp1*i))/max;
         weights[i]:=zmp;
        end;{for}
      end;{Blackman-Nutall}

 7 :begin   //rectangle
     for i:=0 to count div 2 -1-rectFill do weights[i]:=0;
     for i:=count div 2-rectFill to count div 2+rectFill-1 do weights[i]:=1;
     for i:=count div 2+rectFill to count-1 do weights[i]:=0;
    end;{rectangle}

  end;{case}
  MeanAndStdDev(weights,mean,stDev);
  WeightsDraw(weights,form1,form1.panel6,form1.label5,form1.radiogroup7,f1wspx6,f1wspy6,f1centr6);
  WeightsDraw(weights,form2,form2.panel8,form2.label19,form2.radiogroup7,wspx8,wspy8,centr8);
  form2.Label32.Caption:=floatTostrF(mean,ffFixed,8,3);
  form2.Label34.Caption:=floatTostrF(stdev,ffFixed,8,3);

  //-------------------weights spectrum drawings----------------------------------------
  if not form2.checkbox78.Checked then  windowsWeightsSpectrDraw(weights);
 end;{WindowsWeights}

procedure ind_mix_Arr(var index_mix:Tindexa;numdat:word);
 var i,j,k,l:word;
 begin
  j:=0;   l:=high(index_mix);//l - debug
  for i:=0 to numdat-2 do
  begin
   index_mix[i]:=j;
   k:=numdat div 2;
   while k<=j do
   begin
    j:=j-k;
    k:=k div 2;
   end;{while}
   j:=j+k;
  end;{for}
 end;{ind_mix_Arr}

procedure power2(var maxpower:word; numdat:longint);
 var i:longint;
 begin
  maxpower :=0;
  i:=numdat;
  while i<>1 do
  begin
   maxpower:=maxpower+1;
   i:=i div 2;
  end;
 end;{power2}

procedure SinCosArrays(var cosary,sinary:array of double;numdat:word);
 var i:word; harm:real; ce,se:Extended;
 begin
  harm:=2*pi/numdat;
  for i:=0 to numdat-1 do
   begin
    SinCos(harm*i,se,ce);
    cosary[i]:=ce;// cos(harm*i);
    sinary[i]:=se;// sin(harm*i);
   end;
 end;{SinCosArrays}

PROCEDURE swap(VAR s1,s2: realtype);
VAR
  temp: realType;

BEGIN
  temp := s1;
  s1 := s2;
  s2 := temp;
END;


{---------------------SZYBKA TRANSFORMATA DLA DANYCH ZESPOLONYCH --}

PROCEDURE fft_calc(var xreal:rAryType;var yimag:rAryType;numdat:integer);

VAR
  maxpower,arg,cntr,pnt0,pnt1,i,j,a,b,k: integer;
                    prodreal,prodimag,harm: real;
                                 cosary,sinary: rAryType;

PROCEDURE swap(VAR s1,s2: realtype);
VAR
  temp: realType;

BEGIN
  temp := s1;
  s1 := s2;
  s2 := temp;
END;

BEGIN
//  new(cosary);
//  new(sinary);
  j := 0;
  for  i := 0 to numdat - 2 do
  begin
    if i < j then
    BEGIN
      swap(xreal[i],xreal[j]); //usunięto ^ ^
      swap(yimag[i],yimag[j]); //j.w.
    END;
    k := numdat div 2;
    while k <= j do
    begin
       j := j - k;
       k := k div 2;
    end;
    j := j + k;
  end;
  maxpower := 0;
  i := numdat;
  while i <> 1 do
  begin
    maxpower := maxpower + 1;
    i := i div 2;
  end;
  harm := 2 * pi / numdat;
  FOR i:= 0 TO numdat - 1 DO
  BEGIN
      sinary[i] := -sin(harm * i );   //usunięto ^ ^
      cosary[i] := cos(harm * i );    //j.w.
  END;
  a:=2 ; b := 1;
  FOR cntr:= 1 TO maxpower DO
  BEGIN
    pnt0 := numdat DIV a;
    pnt1 := 0;
    FOR k := 0 TO b-1 DO
    BEGIN
      i:= k;
      WHILE i<numdat DO
      BEGIN
        arg := i + b;
        IF k = 0 THEN
        BEGIN
          prodreal := xreal[arg];  //usunięto ^ ^
          prodimag := yimag[arg];
        END
        ELSE
        BEGIN
          prodreal := xreal[arg] * cosary[pnt1] - yimag[arg] * sinary[pnt1];     //usunięto ^ ^ ^ ^
          prodimag := xreal[arg] * sinary[pnt1] + yimag[arg] * cosary[pnt1];     // j.w.
        END;
        xreal[arg] := xreal[i] - prodreal;     //usunięto ^ ^
        yimag[arg] := yimag[i] - prodimag;      //usunięto ^ ^
        xreal[i]   := xreal[i] + prodreal;      //usunięto ^ ^
        yimag[i]   := yimag[i] + prodimag;       //usunięto ^ ^
        i  := i + a;
      END;
      pnt1 := pnt1 + pnt0;
    END;
    a := 2 * a;
    b := b * 2;

  END;
 // setLength(cosary,0);
 // setLength(sinary,0);
END;


{---------------------SZYBKA TRANSFORMATA DLA sygnalow RZECZYWISTYCH---------}
{tablice sin  cos i indeksow mix potega 2 sa przekazywane z zewnatrz}
{co pozwala urzymywac je jako parametry stale w sesji }


PROCEDURE rfft_calc(const index_mix:Tindexa;maxpower:word;
                var sinar,cosar:array of double; var xreal,yimag:array of double;numdat:word);

VAR
  arg,
  cntr,pnt0,pnt1,i,
  j,a,b,k:                       integer;
  prodreal,prodimag:      real;
BEGIN
  for  i := 0 to numdat - 2 do
  begin
    j:=index_mix[i];
    if i < j then
     begin
      swap(xreal[i],xreal[j]);
      swap(yimag[i],yimag[j]);// blokada 14.12.12, bo yimag nie dostarcza danych; usunięta 210113
     end
  end;{for i}
  a:=2 ; b := 1;
  FOR cntr:= 1 TO maxpower DO
  BEGIN
    pnt0 := numdat DIV a;
    pnt1 := 0;
    FOR k := 0 TO b-1 DO
    BEGIN
      i:= k;
      WHILE i<numdat DO
      BEGIN
        arg := i + b;
        IF k = 0 THEN
        BEGIN
          prodreal := xreal[arg];
          prodimag := yimag[arg];
        END
        ELSE
        BEGIN
          prodreal :=  xreal[arg] * cosar[pnt1] + yimag[arg] * sinar[pnt1];
          prodimag := -xreal[arg] * sinar[pnt1] + yimag[arg] * cosar[pnt1];
        END;
        xreal[arg] := xreal[i] - prodreal;
        yimag[arg] := yimag[i] - prodimag;
        xreal[i]   := xreal[i] + prodreal;
        yimag[i]   := yimag[i] + prodimag;
        inc(i,a);
      END;
      inc(pnt1,pnt0);
    END;
    a := 2 * a;
    b := b * 2;
  END;
END;

Initialization
 checkBreak:=false;
END.



