//{$F+,O+,R-,S+,I-,D+-,L-,N+}

unit spectanalisier;
INTERFACE
uses //dos,crt,graph,printer,
sysUtils,common,friend,fft,disco,disco2,grafiti0,grafiti1,
     grafiti4,porto,grafiti3;

var
                             signal,fazR1,fazR2,fazI1,fazI2 : TdbArr;
   unwrap_iter,prtscr,cpstr,phase_corr,
      cmplxCpstr,spectr,Tribolett,pas_odd,pas_even,cut_init,
          interp,ToDiskAvail,firstIn,fall,sort,buzz_to_disk,
                      hiss_buzz_header,TSD,TSP,phTest,quiet : boolean;
       power_cmplx_cpstr,cr_m,buzz_rms_skok,config,init_cut,
                                          scr_nr_as,inv,poz : byte;
                         pix_as,pix_as1,pixY,pix1Y,minLevel,
                 freq_interval,rms1,rms,buzz_1,f_low,f_high : real;
                                            dOm,dOmKw,delta : realType;
                                       zm_1,zm_2,weightMean : extended;
              corr_spectr_nb,ph_corr_nb,time_change,nr_konc : longint;
                                              Cpik,thld,dev : single;
                                               F0idx1,F0idx : integer;
                                                        lst : text;
                                               printFileDir : string='C:\apl\spectrum\results\lst.txt';

procedure spl(var xr:array of double; var rms1,rms:real);
procedure druk(cynk:string);
procedure spectr_moment;
procedure buzz_hiss_files_open(storedDataNb:longint);
procedure buzz_hiss_files_close;
procedure standard_Ins(k:word);
procedure getCpstr;
procedure PutCpstr;
procedure cepstrum;
procedure filtr(var xr,yi:array of double; count:word; var fcpstr:realType;
                 var j:integer;var l,m:word);
procedure tran_rev(var stat : Text);
procedure invers(var xr,yi:array of double);

IMPLEMENTATION
 uses unit1, visualization;
var                  t1 : extended;
     p_rms,p_buzz,d_var : real;


type crctype =  array[1..50] of realType;
var circleArr : ^crcType;

procedure spl(var xr:array of double; var rms1,rms:real);
 var i:word;
 begin
 p_rms:=rms; rms:=0;
  for i:=0 to count-1 do rms:=rms+sqr(xr[i]);
  rms:=sqrt(rms)/(WeightMean*count);
  if rms=0 then rms:=rms1;
 end;{spl}

procedure quicksort(var a:array of double; var indexSort:Tindexa;Low,High: integer);
 var i:integer;
procedure sort(l,r: integer);
var
  i,j,k : integer; x : realType;
begin
  i:=l; j:=r; x:=a[indexSort[(l+r) DIV 2]];
  repeat
    while a[indexSort[i]]<x do i:=i+1;
    while x<a[indexSort[j]] do j:=j-1;
    if i<=j then
    begin
      k:=indexSort[i]; indexSort[i]:=indexSort[j]; indexSort[j]:=k;
      i:=i+1; j:=j-1;
    end;
  until i>j;
  if l<j then sort(l,j);
  if i<r then sort(i,r);
end; {sort}

begin
 for i:=0 to count-1 do indexSort[i]:=i;
  sort(Low,High);
end {quicksort};

procedure Tribolett_Phase_Unwrap(var virt_vect:array of double);
{
rozwija faze metoda Triboletta - Bir Bham - McClellan

IEEE ASSP-25,2, 1977, 170-177 Tribolett
IEEE ASSP-26,1, 1978, 104-105 Bonzanigo
IEEE ASSP-28,5, 1980, 583-585 Bir Bham & McClellan
}

 var
 x0,y0,x1,y1,x2,y2,sqrx,sqry,modul,
                  A0,ph,arg,dOmAcc : extended;
                               i : word;
                             logik : boolean;
                 infoStr1,infostr2 : DirStr;
                           adaptNb : word;
 const
       stlim=31;
       accuracy:longint=1073741824; {2 do potegi stlim-1}

procedure infos(zew:byte);
 begin
  inc(adaptNb);
  if not quiet then
   begin
    str(zew:1,InfoStr1);
    str(i:1,InfoStr2);
    InfoStr2:=infoStr1+' i='+InfoStr2;
    str(adaptNb:1,InfoStr1);
    infoStr1:='; lbaad='+InfoStr1;
    informator('Iteracje, zew '+InfoStr2+InfoStr1,120*(zew div 10)+80,12,0,0);
//    if KeyPressed then if ReadKey='q' then
     begin
      quiet:=true;
//      nosound
     end;
   end;{zew}
   end;{infos}

procedure DFT(i,k:longint);
{
oblicza spectrum i I-wsza pochodna w pkcie i-1+(k-1)/accuracy
i-1, bo obliczenia dotycza poprzedn. przedzialu czestotl.
}
var j:word;
 begin
  A0:=pi2*(i-1+(k-1)/accuracy)/count;
    x1:=signal[0];  y1:=0;
    x2:=0;          y2:=0;
    for j:=0 to count-1 do
     begin
      arg:=A0*j;
      x0:=cos(arg)*signal[j];
      x1:=x1+x0; x2:=x2+x0*j;
      y0:=-sin(arg)*signal[j];
      y1:=y1+y0; y2:=y2+y0*j;
     end;{for}
     x1:=x1/(WeightMean*count);   y1:=y1/(WeightMean*count);  {spectrum}
     x2:=x2/(WeightMean*count);   y2:=y2/(WeightMean*count);  {I-wsza poch. po czestotl.}
 end;{DFT}

procedure adapt(i:word;prevDeriv,actDeriv,oldPhase:extended;
           var signal:array of double; var ph:extended;var adaptNb:word;zew:byte);

{
oblicza posrednie wartosci fazy
}
 label 100,110,120,130,140;


 var                            isk : array[1..stlim] of longint;
                                sk2 : array[1..stlim] of realType;
                        phinc,delta : extended;
                               j,sp : word;
                               ib,k : longint;
                              b1,b2 : realType;

 begin

  {---------zainicjowac stos}
  sp:=1;   j:=0;
  isk[1]:=accuracy+1; {dokladnosc}
  sk2[1]:=actDeriv;   {pochodna}

  {--------zainicjowac rejestry}
  ib:=1;
  b1:=oldPhase;      {poprzednia faza}
  b2:=prevDeriv;     {poprzednia pochodna}
100:
  if isk[sp]-ib>1 then goto 110
    else
     begin
      str(round(ln(count)/ln(2)+1):1,InfoStr1);
      informator('Accuracy too small for correct phase unwrap -'+'press '+
       infoStr1,80,12,200,200);
      goto 140;
     end;{else}

  {--------czestotliwosci posrednie w=2pi/N*(i-1+(k-1)/accuracy)}
110:
  infos(10*zew+j);
  k:=(isk[sp]+ib) div 2;
//  if not quiet then sound(round((i-1+(k-1)/accuracy)*rate/count));

  {--------obliczyc DFT dla czestotliwosci posrednich}
    DFT(i,k);
  {--------obliczyc pochodna}
  if sp<stlim then inc(sp)
  else
   begin
    informator('Tablice stosu za male powieksz stlim',80,12,200,50);
    goto 140;
   end;{else}
  isk[sp]:=k;
  sk2[sp]:=-(x1*x2+y1*y2)/(sqr(x1)+sqr(y1));
  logik:=(isk[sp]-ib>1) and (sp<stlim);

130:{--------obliczyc faze}
  delta:=dOmacc*(isk[sp]-ib);
  phinc:=delta*(b2+sk2[sp]);
  if abs(phinc)>thld then begin j:=1; goto 100 end;
  ph:=oldPhase+phinc;
140:
  ib:=isk[sp];
  b1:=ph;
  b2:=sk2[sp];
  dec(sp);
  {
  -------------jesli stos jest pusty, to faza rozwinieta dla pktu
               w=2pi*i/N jest trzymana w rejestrze b1
  }
  if sp>0 then goto 130;
 end;{adapt}

procedure fft(var fazR,fazI:array of double);
var  i:word;

 begin
  for i:=0 to count-1 do fazI[i]:=0;
  rfft_calc(index_mix,maxpower,sinar,cosar,fazR,fazI,count);
  for i:=0 to count-1 do
   begin
    fazR[i]:=fazR[i]/(WeightMean*count);
    fazI[i]:=fazI[i]/(WeightMean*count);  {zmiana znaku}
   end;{for}
  end;{fft}

var
oldPhase,prev_deriv,act_deriv,phinc : extended;

label 1000;

 begin
  if grafika and moovie and cpstr_to_scr then
//   informator('phase unwrap processing',8,getMaxY-24,200,200);
  dOmAcc:=0.5*dOm/accuracy;
  for i:=0 to count-1 do                                   {Pochodne, dane do}
   begin
    fazI1[i]:=i*signal[i];
    (*fazR2^[i]:=sqr(i)*signal^[i];*)
   end;{for}
                      {pierwsza pochodna}
  fft(fazI1,fazR1);   {zamienic tablice miejscami, bo I poch.=-j*fft(i*x[i])}
(*fft(fazR2,fazI2);   {druga pochodna}*)
                      {calka}
  x0:=0;
  for i:=0 to cp do                   {I-wsza pochodna; ladowac do vitr_vect}
   begin
    x0:=xr[i];        y0:=yi[i];   sqrx:=sqr(x0); sqry:=sqr(y0);
    x1:=-fazR1[i];    y1:=fazI1[i];
    modul:=sqrx+sqry;
    if modul>1e-39 then
     begin
      ph:=-(x0*y1-y0*x1)/modul;                  {I-wsza pochodna}
      x0:=x0+ph;
      virt_vect[i]:=ph;
(*
      fazR2^[i]:=                                {II-ga  pochodna}
       (modul*(-x0*fazI2^[i]+y0*fazR2^[i])+
       +2*x0*y0*(sqr(x1)-sqr(y1))+
       +2*x1*y1*(sqry-sqrx))/sqr(modul)
*)
     end
    else begin virt_vect[i]:=0; fazR2[i]:=0 end;
   end;
  ph:=0;
  oldPhase:=ph;   adaptNb:=0;   prev_deriv:=virt_vect[0];
  oldPhase:=atan(xr[1],yi[1]);   adaptNb:=0;   prev_deriv:=virt_vect[1];
  virt_vect[1]:=oldPhase;
  if not quiet then informator('q=cicho!',0,12,100,200);
  for i:=2 to cp do
   begin
    act_deriv:=virt_vect[i];
    phinc:=0.5*dOm*(act_deriv+prev_deriv);
(*          -(dOmKw/12)*(fazR2^[i]-fazR2^[i-1]);*)
    if unwrap_iter then
     if abs(phinc)>thld then
      begin
       adapt(i,prev_Deriv,act_Deriv,oldPhase,signal,ph,adaptNb,1);
       goto 1000
      end;
    ph:=oldPhase+phinc;
1000:
    prev_deriv:=act_deriv;
    virt_vect[i]:=ph;
    oldPhase:=ph;
  end;
  virt_vect[0]:=0;
  x0:=virt_vect[cp]/cp;
  for i:=1 to cp do                         {usunac skladnik liniowy}
   virt_vect[i]:=virt_vect[i]-i*x0;
  for i:=cp to count-1 do
   virt_vect[i]:=-virt_vect[count-i];        {uzupelnic druga polowe}
//  nosound;
  str(adaptNb:1,InfoStr1);
  informator('liczba wyw. DFT= '+iNfoStr1,80,12,0,0);
  if grafika and moovie and cpstr_to_scr then
//   informator('                       ',8,getMaxY-24,0,0);
 end;{Tribolett_Phase_Unwrap}

procedure APL_Phase_Unwrap(var virt_vector:array of double);
var i:word;
 oldPhase,phi,phi1,phi2 : extended;
 begin
  virt_vect[0]:=0;
  oldPhase:=0;
  for i:=1 to cp do
   begin
    phi:=int(oldPhase/pi);
    phi:=phi*pi+atan(xr[i],yi[i]);
    phi:=phi+int((oldPhase-phi)/pi)*pi;
    if abs(phi-oldPhase)>pi/2 then
     begin
      phi1:=phi+pi;
      phi2:=phi-pi;
      if abs(phi1-oldPhase)<abs(phi2-oldPhase) then phi:=phi1
      else phi:=phi2;
     end;
    virt_vect[i]:=phi;
    oldPhase:=phi;
   end;{unwraping}
   phi:=virt_vect[cp]/cp;
   if phi<>0 then
    for i:=1 to cp do                        {delete linear component}
     virt_vect[i]:=virt_vect[i]-i*phi;
   for i:=cp to count-1 do
    virt_vect[i]:=-virt_vect[count-i];     {create secound half}
 end;{Apl_Phase_unwrap}



procedure phase_correction(var xr,yi,virt_vect:array of double;
 var zm_1,zm_2:extended;ip,ik:word;var firstIn:boolean; sort,fall:boolean;
 var stat:Text);

 var                              // indexSort : indexa; blokada 01052014
                                     j,l,m : word;
                                       break : boolean;
                                       coeff : byte;
             pi2mc,pi2mcj,gm,gm11,gm12,delta,
                 modj,alfa,phij,gmj,maxj : extended;
                                 xtemp,ytemp : realType;

 function REV_DFT(m:word):extended;
  var s1,scp,pi2kmc:extended; k:word;
  begin
   s1:=0;
   pi2mc:=pi2*m/count;
   if odd(m) then scp:=-xr[cp] else scp:=xr[cp];
   for k:=1 to cp-1 do
    begin
     pi2kmc:=k*pi2mc;
     s1:=s1+xr[k]*cos(pi2kmc)-yi[k]*sin(pi2kmc);
    end;
    REV_DFT:=xr[0]+2*s1+scp;
  end;{REV_DFT}

 procedure correct_components;
  begin                         {---------------skorygowac xr i yi}
     xr[j]:=cos(phij)*modj; xr[count-j]:=xr[j];
     yi[j]:=sin(phij)*modj; yi[count-j]:=-yi[j];
  end;{correct_components}

 BEGIN
  {
  1. obliczyc wartosc sygnalu gm w punkcie m wprost ze wzoru DFT i roznice
  pomiedzy gm a wartoscia w tym punkcie wynikajaca z obliczenia poprzedniego
  odcinka
  }
  if ik>=count-2 then exit;  {korekta fazy potrzebuje 2 pktow spoza odcinka}
  inc(corr_spectr_nb);
  m:=ip;
  pi2mc:=pi2*m/count;
  gm:=REV_DFT(ip);
  if firstIn then begin zm_1:=gm; zm_2:=REV_DFT(ip+1);firstIn:=false end;
  delta:=gm-zm_1;   {zapewnic pierwsza wartosc=0 przez wpisanie 0 na poczatku}
                    {analizowanego sygnalu}

  {-----------------zbior dla obliczen statystycznych skutecznosci metody}
  with restat do
   if retNd then
   begin
    delta_0:=zm_0-gm;
    s_d_0:=s_d_0+delta_0; k_s_d_0:=k_s_d_0+sqr(delta_0);
    s_d:=s_d+delta;       k_s_d:=k_s_d+sqr(delta);
    abss_d_0:=abss_d_0+abs(delta_0);
    abss_d:=abss_d+abs(delta);
    zm_0:=rev_dft(ik+1);
   end;{reStat}

  {
  2. jesli delta<> 0 to posortowac indeksy wg wartosci modulu odwracanej
   transformaty, i skorygowac fazy skladowych tak, aby gm=zm_1 z zachowaniem
   pierwotnej wartosci modulu
  }
  if delta<>0 then
   begin
    if sort then quickSort(virt_vect,indexSort,1,cp);
    if fall then l:=cp else l:=1;
    repeat
     break:=false;
     if sort then j:=indexSort[l] else j:=l;
     modj:=sqrt(sqr(xr[j])+sqr(yi[j]));
     if modj>0 then
      begin
       pi2mcj:=pi2mc*j;

     {---------------obliczyc udzial j-ej skladowej w sygnale wynikowym}
       phij:=atan(xr[j],yi[j])+pi2mcj;
       gmj:=modj*cos(phij);

     {---------------obliczyc max zdolnosc przesuniecia przez j-ta skladowa}
       if delta>0 then maxj:=modj+gmj
       else            maxj:=modj-gmj;
       if maxj<>0 then
        begin
//         if not quiet then sound(round(j*rate/count));
         inc(ph_corr_nb);
         if j=cp then coeff:=1 else coeff:=2;
         maxj:=coeff*maxj;
         gmj:=coeff*gmj;
     {---------------obliczyc nowy kat fi}
         if abs(delta)<abs(maxj) then
          begin
           alfa:=arcCos((zm_1-gm+gmj)/(coeff*modj));
           phij:=atan(xr[j],yi[j]);
           if phij<>pi then
            begin
             if abs(alfa-pi2mcj-phij)<abs(pi2-alfa-pi2mcj-phij) then
              phij:=alfa-pi2mcj
             else phij:=pi2-alfa-pi2mcj;
             Correct_Components;
            end
           else
            begin
             phij:=alfa-pi2mcj;
             xtemp:=xr[j]; ytemp:=yi[j];
             Correct_Components;
             gm11:=Rev_DFT(ip+1);
             xr[j]:=xtemp; yi[j]:=ytemp;
             phij:=pi2-alfa-pi2mcj;
             Correct_Components;
             gm12:=Rev_DFT(ip+1);
             if abs(gm12-zm_2)>abs(gm11-zm_2) then
              begin
               phij:=alfa-pi2mcj;
               Correct_Components;
              end;{if}
            end;{else phij=pi}
           gm:=gm-gmj+coeff*modj*cos(alfa);
           break:=true;
          end
         else
          begin
           if delta<0 then
            begin
             phij:=-pi2mcj;
             gm:=gm+maxj;
            end
           else
            begin
             phij:=pi-pi2mcj;
             gm:=gm-maxj;
            end;
           correct_components;
          end;{else}
         delta:=gm-zm_1;
        end;{maxj<>0}
      end;{modj>0}
     if fall then dec(l) else inc(l);
     if fall then break:=(l=0) or break else break:= (l>=cp) or break;
    until (gm=zm_1) or break or (abs(delta)<2e-11);
   end;{delta<>0}
  zm_1:=REV_DFT(ik+1); zm_2:=REV_DFT(ik+2);
//  nosound;
 END;{phase_correction}


procedure druk(cynk:string);
 var i:word;
 begin
  if TtextRec(lst).mode=fmclosed then
   begin
    assignFile(lst,printFileDir);
    rewrite(lst);
   end;
  writeLn(lst,'==========================================================');
  writeln(lst,cynk);
  str(inv,lancuch); lancuch:=', invers='+lancuch;
  lancuch:='filtr'+tryb(filtrOn)+', korekcja fazy'+tryb(phase_corr)+lancuch;
  writeln(lst,lancuch);
  write(lst,'_NR_|_____ X ____|_____ Y ______|__ Virt_Vect _|'#10#13);
  for i:=0 to count-1 do
    write(lst,i:4,'|',xr[i]:11:6,' |  ',yi[i]:11:6,' |  ',virt_vect[i]:11:6,' | ',#10#13);
  write(lst,'=============================================='#10#10#13);
  flush(lst);
 end;{druk}

procedure spectr_moment;
 var i:longint;
 begin
  //if sing or only then konvert_rel(Gain,xr,count,t5,t5)    w WAV nie konwertujemy
 // else           konvert_rel(Gain,xr,count,ymin,ymax);
  if winNbr<>0 then
   for i:=0 to count-1 do xr[i]:=xr[i]*windowWeights[i];
  signal:=xr;
  spl(xr,rms1,rms);                               {natezenie dzwieku - RMS}
  if not Change_insert then
  if krok and grafika and moovie and not (sing or only) and osc_to_scr then
   sig_we_graph;                                  {wykres sygnalu wejsciowego}
  yi:=yZero;
  if count_cp_bool then jot:=count else jot:=cp;
  rfft_calc(index_mix,maxpower,sinar,cosar,xr,yi,count);
  for i:=0 to count-1 do
   begin
    xr[i]:=xr[i]/(WeightMean*count); yi[i]:=yi[i]/(WeightMean*count);
   end;{for}
  if not cpstr then jot:=count;
  modul_npdBrms(xr,yi,virt_vect,angle,count,false,false);
  if prtScr and spectr_to_scr then druk('Wyniki analizy spectrum'#10#13+
   'X - real, Y - imag, V_V - modul, dB');
 end;{spectr_moment}

 procedure Invers(var xr,yi:array of double);
  var i:word;
  begin
   for i:=1 to count-1 do yi[i]:=-yi[i];
   rfft_calc(index_mix,maxpower,sinar,cosar,xr,yi,count);
  end;

 procedure phaseTest(var xr,yi,v_v:array of double);
 {
 Generuje faze w postaci trojkata rownoramiennego 0,-range, 0,range,0
 a nastepnie czesc urojona
 }
  var       i : word;
      range,c : real;
  begin
   range:=100;
   c:=range/(cp div 2);
   for i:=0 to cp div 2-1 do
    v_v[i]:=-i*c;
   for i:=cp div 2 to cp do
   v_v[i]:= c*i-2*range;
   for i:=cp to count-1 do
    v_v[i]:=-v_v[count-i];     {create secound half}
   cpstr_gr1(true,fIdx,'Sztuczna faza');
   xr[0]:=1;
   yi[0]:=0;
   for i:=1 to cp do             {create imaginary part}
    begin
     xr[i]:=-1;
     yi[i]:=sin(v_v[i])/cos(v_v[i]);
    end;
   for i:=cp to count-1 do        {create secound half}
    begin
     yi[i]:=-yi[count-i];
     xr[i]:=xr[count-i];
    end;
   for i:=0 to count-1 do  //20052013 jawne przeładowanie z powodu niezgodności typów
    begin
     signal[i]:=xr[i];
     v_v[i]:=yi[i];       {create origin signal}
    end;
   invers(signal,v_v);
  end;{phaseTest}



 procedure xr_retriev(var stat:Text);
 var i:integer;
     tim : realType;
 begin
  for i:=0 to count-1 do          {odzyskanie czesci rzeczywistej i urojonej}
   begin
    t1:=exp(xr[i]);              {RMS}
    tim:=yi[i];
    yi[i]:=t1*sin(tim);          {imag}
    xr[i]:=t1*cos(tim);          {real}
    virt_vect[i]:=t1;            {RMS}
   end;{for}
  if phase_corr then
   begin
    phase_correction(xr,yi,virt_vect,zm_1,zm_2,ip,ik,firstIn,sort,fall,stat);
    lancuch:= 'bez korekcji fazy skladowych spectr.'#10#13;
   end
  else lancuch:='po korekcji fazy skladowych spectr.'#10#13;
  lancuch:='RETRIEV; Wektory przed transformacja '+
              'fft odwracajaca '#10#13'spectrum filtrowane cepstralnie; '+
               lancuch+
              'X - real, Y - imag, V_V - RMS';
  if prtScr and ret_to_scr then druk(lancuch);
  if krok and not sing and grafika and moovie and ret_to_scr then
   if not Change_insert then ret_gr1(freq_interval);
  invers(xr,yi);
  if not phase_corr and (ik<count-2) then
   begin zm_1:=xr[ik+1]; zm_2:=xr[ik+2] end;
  if prtScr and ret_to_scr then
   druk('RETRIEV; Wektory po transformacji odwracajacej;'+
        'sygnal filtrowany cepstralnie'+
   #10#13'X, Y - czesc rzeczyw. i uroj. wynikowej postaci czasowej'#10#13+
   'V_V - modul spectrum filtrowanego cepstralnie - RMS');
  if grafika and moovie and ret_to_scr then
   if not Change_insert then ret_gr2;
 end;{xr_retriev}

 procedure filtr(var xr,yi:array of double; count:word; var fcpstr:realType;var j:integer;var l,m:word);
  var i:longint;

  procedure LowPas;
 {
  wycina srodek zaczynajac od pozycji na ktorej stoi kursor
  symetria wokol punktu (count/2+1) ! punkt 0 nie moze byc wyciety !
 }
   var i:integer;
   begin
    for i:=j to count-j do
     begin
      xr[i]:=0; yi[i]:=0;
     end
   end; {LowPas}

  procedure HiPas;
{
pozostawia srodek, obcina boki
nie wycina punktu 0 ! wynika to z zasady symetrii
}
   var i:integer;
   begin
    for i:=1 to j-1 do
     begin
      xr[i]:=0; yi[i]:=0;
     end;{for}
     begin
      for i:=count-j+1 to count-1 do
       begin
        xr[i]:=0; yi[i]:=0;
       end{for}
     end;
   end;{HiPas}

  procedure EvenPas;
  {wycina skladowe nieparzyste}
   var i:integer;
   begin
    i:=1;
    while i<=count-1 do
     begin
      xr[i]:=0; yi[i]:=0;
      inc(i,2);
     end;
   end; {EvenPas}

  procedure OddPas;
  {wycina skladowe parzyste}
   var i:integer;
   begin
    i:=0;
    while i<=count-1 do
     begin
      xr[i]:=0; yi[i]:=0;
      inc(i,2);
     end;
   end; {OddPas}

 procedure StartPas;
  var i:integer;
   begin
    for i:=1 to j do
     begin
      xr[i]:=0;
      yi[i]:=0;
     end;{for}
   end;{StartPas}

 procedure endPas;
  var i:integer;
   begin
    for i:=j to count-1 do
     begin
      xr[i]:=0;
      yi[i]:=0;
     end;{for}
   end;{endPas}

  begin
    t1:=0;                             {utrata mocy}
    if pas_odd then OddPas;
    if pas_even then EvenPas;
    if cut_init then
     begin
      xr[init_cut]:=0;
      yi[init_cut]:=0;
      xr[count-init_cut]:=0;
      yi[count-init_cut]:=0;
     end;
   if fcpstr>0 then
    if cpstr then
     begin
      if double_side_cpstr_filtr then
       if j>cp then begin j:=cp; fcpstr:=rate/cp end
       else
      else if j>count-1 then begin j:=count-1; fcpstr:=rate/(count-1) end;
      if j<=0 then begin fcpstr:=rate/count-1; j:=count-1 end;
      if count_cp_bool then begin l:=count-j;m:=count-1 end
      else begin l:=cp-1;m:=cp-1 end;
      if double_side_cpstr_filtr then
       if HiFiltr then LowPas else HiPas
      else if HiFiltr then StartPas else EndPas;
     end
    else
     begin
      if count_cp_bool then begin l:=count-j; m:=count-1 end
      else begin l:=cp-1; m:=cp-1 end;
      t1:=0;
      if HiFiltr then HiPas else LowPas;
     end
   else
    begin
     j:=cp;
    for i:=0 to count-1 do  //20052013 jawne przeładowanie z powodu niezgodności typów
     begin
      xr[i]:=yZero[i];
      yi[i]:=yZero[i];
     end;
    end;
   if j<0 then j:=0 else if j>count then j:=count;
   if prtScr then druk('FILTR; Wektory j.w. po filtracji');
  end;{filtr}

procedure I_O_Ring(i:longint);
 {
 pierscien opozniajacy
 pobiera i laduje dane do bufora pierscieniowego
 }
  var wsk:byte; tim : realType;
  begin
   wsk:=(i-1) mod cr_m+1;
   tim:=circleArr^[wsk];            {odczyt}
   circleArr^[wsk]:=rms;     {zapis}
  end;{I_O_Ring}


procedure buzz_hiss;

 var i,j,k:word;
 indexSort:Tindexa;
     t1,t2:extended;

function averCircle:realType;
{
srednia antyfluktuacyjna
}
 var i:byte;
 begin
  t1:=0;
  for i:=1 to cr_m do t1:=t1+circleArr^[i];
  averCircle:=t1/cr_m
 end;{averCircle}

function deviation(min:word):extended;
{
Oblicza odchylenie od rozkladu rownomiernego na podstawie dystrybuanty;
dystrybuanta[i]=max_index-min+1
}
var      i,j : integer;
         tim : extended;

function max_index(i:word):word;
{$B-} {bo dla i=count-1 nastapiloby przesterowanie indeksu!}
 begin
  while (i<cp) and (virt_vect[indexSort[i]]=virt_vect[indexSort[i+1]]) do
  inc(i);
  max_index:=i;
 end;{max_index}

 begin
  i:=min;
  quickSort(virt_vect,indexSort,i,cp);
  tim:=0;
  i:=max_index(min);
  repeat                               {obliczenie calki po dystrybuancie}
   j:=max_index(i+1);
   tim:=tim+(i+j+2)*(virt_vect[indexSort[j]]-virt_vect[indexSort[i]]);
   i:=j
  until i>=cp;
  i:=-(min-1)+cp;                                 {calkowita liczba obserwacji}
  t1:=virt_vect[indexSort[min]];
  t2:=virt_vect[indexSort[cp]];
  if (t1<>t2) and (tim<>0) then tim:=(i+1)*(t2-t1)/tim
  else tim:=0;
  dev:=tim;
  deviation:=tim;
 end;{deviation}

var mx : realType;
   tim : realType;
 procedure rc_buzz;
 {
 oblicza wartosc buzz sprawdzajac zakresy zmiennych -
 performs range checking
 }
 
  begin
   if mx<=sqrt(1e38) then
      if t2<=1e38/4 then
       if 4*t2<=sqr(mx)*1e38 then
        if abs(t2)>0 then
         if 4*t2>=sqr(mx)*1e-38 then  buzz:=(4*t2)/sqr(mx)              {***}
         else buzz:=1e-38
        else buzz:=1e-38
       else buzz:=1e38
      else begin buzz:=sqr(mx)/4; buzz:=t2/buzz end
     else
      begin
       buzz:=mx/2; buzz:=sqrt(t2)/buzz;
       if buzz<=sqrt(1e38) then buzz:=sqr(buzz)
       else buzz:=1e38
      end;
  end;{rc_buzz}

  begin
   p_buzz:=buzz;
   mx:=virt_vect[cp-1];
   i:=round(rate/f_high); k:=cIdx;         {pozycja w cepstrum}
   if k>i then k:=i;
   j:=cp-1-k;
   if j>0 then                      {stosunek wariancji do kwadratu rozstepu}
    begin
     t1:=0; t2:=0; powerLost:=0;
     for i:=cp-1 downto k do
      begin
       tim:=virt_vect[i];
       powerLost:=powerLost+tim;
       if mx<tim then mx:=tim;
       t1:=t1+tim; t2:=t2+tim*tim;
      end;
     t1:=t1/j;
     t2:=t2/j-t1*t1;
     rc_buzz;
     t1:=powerLost;
     for i:=0 to k-1 do t1:=t1+virt_vect[i];
     powerLost:=powerLost/t1;
     d_var:=buzz;
     buzz:=buzz*deviation(cIdx);
     t2:=1.15+0.04*(WeightMean*count)+(0.007+0.04*(WeightMean*count))*buzz_rms_skok;  {to jest parametr staly!!}// 18102014 fill_win zastąpiono przez suma_wag to może być źle!
     if not hiss then rms1:=averCircle;                           {antyflukt}
     if rms1/rms>t2 then  hiss:=true
     else
      begin
       if hiss then                                {nastapila zmiana stanu}
        begin
         for i:=1 to cr_m do circleArr^[i]:=rms;   {zapoczatkowac bufor}
         hiss:=false;
        end;
       buzz_1:=buzz;
       I_O_Ring(nr_konc);    {odpowiednik rms1:=rms; nr_konc - to nr prazka!}
      end;{else}
    end{if j>0}
   else buzz:=0;
  end;{buzz_hiss}

 procedure buzz_hiss_tun_print(var i:integer; tim:realType);
 var
             logo : boolean;
 begin
  {$I-}
  repeat
   if buzz_to_disk then
    begin
//     close(lst);                       {wydruki na dysk}
     logo:=fileExists(patternPath+'.buz');
//      assignFile(lst,patternPath+'.buz');
//     if logo then append(lst) else rewrite(lst);
    end;{buzz_to_disk}
//   writeln(lst);
   i:=IOResult;
   if i<>0 then
    begin
//     bar(56,getMaxY-8,422,getMaxY);
//     if buzz_to_disk then outTextXY(56,getMaxY-8,'Disk not ready!')
//     else outTextXY(56,getMaxY-8,'Printer is off!;');
//     outTextXY(184,getMaxY-8,'  esc, ^F7 lub ^C - uciekamy!');
     write(#7);
//     if KeyPressed then  {wylaczyc wydruki?}
      begin
//       p_ex:=readkey;
//       if (p_ex=#0)  and (readkey=#100) or (p_ex=#3) or (p_ex=#27) then
        begin tsd:=false; config:=4; config_show:=true; exit end;
      end;{exit}
    end;
  until i=0;
  if hiss_buzz_header then
   begin
    hiss_buzz_header:=false;
//    getDate(ye,mo,da,daw);
//    writeln(lst);
//    writeln(lst,'----------------------------------');
//    write(lst,'dn. ',da,'.',mo,'.',ye);
//    getTime(ye,mo,da,daw);
//    writeln(lst,' o godz. ',ye,'.',mo);
//    write(lst,'lba pkt. fft ',count,' okno ');
(*    case winNbr of
     0 : write(lst, 'prostokat');
     1 : write(lst,'Gaussa');
     2 : write(lst,'Keisser-Bessel');
     3 : write(lst,'Blackman-Harris');
     4 : write(lst,'top-flat')
    end;{case}*)
//    writeln(lst,' skok ',skok:1:2,' ms',' filtr:',f_low:1:0,', ',f_high:1:0,' Hz'#13#10,opis);
//    writeln(lst,'=================================='); writeln(lst);
  //  writeln(lst,'nr |','  RMS  |RMS2/RMS1|','RMS1/RMS2|','odch/rozst|',' dev.|',
    //' ton_szum|','t_s2/t_s1|','t_s1/t_s2|',' F0  | stan | powod');
//    repeat inc(i); write(lst,'=') until i>=105; writeln(lst);
   end;{hiss_buzz_header}
//  write(lst,nr_konc:3,'|',rms:7:4,'|',p_rms/rms:9:3,' |',rms/p_rms:9:3,' |',
//   d_var:9:4,' |',tim:5:3,'|',buzz:8:3,' |',buzz/p_buzz:8:3,' |',
//   p_buzz/buzz:8:3,' |',F0:5:0,'|');
//  if hiss then  write(lst,' szum |')
//  else write(lst,' ton  |');
 end;{buzz_hiss_tun_print}

procedure buzz_hiss_files_open(storedDataNb:longint);
 var zmp:array[1..5] of single; i:byte;
 begin
  assignFile(rms_file,patternPath+'.rm'+resFileNb);
  assignFile(Cpik_file,patternPath+'.pk'+resFileNb);
  assignFile(F0_file,patternPath+'.F0'+resFileNb);
  rewrite(rms_file);
  rewrite(Cpik_file);
  rewrite(F0_file);

  {---------------------------------zapisac parametry analiz}
  zmp[1]:=jump;
  zmp[2]:=rate;
  zmp[3]:=storedDataNb;
  zmp[4]:=count;
  zmp[5]:=winNbr;
  for i:=1 to 5 do
   begin
    write(rms_file,zmp[i]);
    write(Cpik_file,zmp[i]);
    write(F0_file,zmp[i]);
   end;
 end;{buzz_hiss_files}

procedure buzz_hiss_files_close;
 begin
  close(rms_file);
  close(Cpik_file);
  close(F0_file);
 end;{buzz_hiss_files_close}

procedure buzz_hiss_files_write;
 begin
  dev:=rms; write(rms_file,dev);
  dev:=F0;  write(F0_file,dev);
  dev:=Cpik; write(Cpik_file,dev);
 end;{buzz_hiss_files_write}

 procedure forDown(var i:integer; j:integer);
  begin
   while (abs(xr[i+1])<=abs(xr[i])) and (i<j) do inc(i); {najbl. minimum}
  end;{ForDown}

 procedure BackDown(var i:integer; j:integer);
  begin                                 {Cofajac sie dalej znalezc minimum}
   while (abs(xr[i-1])<abs(xr[i])) and (i>=j) do dec(i);
  end;{backDown}

 procedure forMin(var i:integer; j:word);
  begin
   while (abs(xr[i+1])>=abs(xr[i])) and (i<j) do inc(i); {przeskocz maks.}
   forDown(i,j);                                           {najbl. minimum}
  end;{forMin}

 procedure backMin(var i:integer; j:word);
  begin                                 {cofajac sie przeskoczyc maximum}
   while (abs(xr[i-1])>=abs(xr[i])) and (i>=j) do dec(i);
   backDown(i,j);                       {Cofajac sie dalej znalezc minimum}
  end;{backMin}

 procedure Auto_Tun(f_low,f_high:real; var k,idx1,idx:integer);
 {
 Automatyczne przestrajanie filtru cepstralnego.
 Procedura Auto_Tun wyznacza punkt graniczny rozdzielajacy pobudzenie
 od odpowiedzi impulsowej.

 Wszystkie operacje sa wykonywane na obszarach okreslonych indeksami;
 korzystamy tu z zaleznosci:

 cIdx = round(rate/freq)       = round(count/fIdx)  - pozycja w cepstrum
 fIdx = round(freq*count/rate) = round(count/cIdx)  - pozycja w spectrum
 freq = fIdx*rate/count        = rate/cIdx          - czestotliwosc
 quef = count/(fIdx*rate)      = cIdx/rate          - quefrency

 }

  const c1=200.0; {[Hz]; c1 real !!}
        c2=6.0;   {[ms]}
        c3=0.6;   {ampli Slope}
        c4=0.6;   {ampli Subh}
        c5=0.4;   {wzgledny przyrost}


  type proc= procedure(var i:integer; j:word);
  var
    i,j,hcp,lcp,hsp : integer;
                 t1 : RealType;
       logic,logico : boolean;

  procedure maximum(l,m:longint; var n:integer);
   begin
    t1:=abs(xr[l]); n:=l; inc(l);
    while l<m do                   {poszukac wartosci najwiekszej w cepstrum}
     begin                         {sposrod rownych bierzemy maksimum}
      if abs(xr[l])>t1 then       {najnizej lezace}
      begin t1:=abs(xr[l]); n:=l end;
      inc(l);
     end;
   end;{maximum}

   procedure maxiS(l,m:longint; var n:integer);
   begin
    t1:=virt_vect[l]; n:=l; inc(l);
    while l<m do                   {poszukac wartosci najwiekszej w spectrum}
     begin                         {sposrod rownych bierzemy maksimum}
      if virt_vect[l]>t1 then     {najnizej lezace}
      begin t1:=virt_vect[l]; n:=l end;
      inc(l);
     end;
   end;{maxiS}

  procedure walker(i,j,k:integer; var m:integer; subr:proc);
   begin
    subr(i,j);
    if (@subr=@forMin) then maximum(i,k,m)
    else maximum(k,i,m);
   end;{walker}

  function forWalk(k:integer; var j:integer):boolean;
   begin
    walker(k,cp,cp,j,forMin);
    forWalk:=(j>=(-2+k)*2) and (j<=(k+2)*2)
   end;{forWalk}

   function backWalk(k,l,m:integer; var n:integer):boolean;
    begin
     walker(k,l,m,n,backMin);
     backWalk:=(k>=(-2+n)*2) and (k<=(n+2)*2)
    end;{backWalk}

  procedure subh(k,l:integer; var m:integer);
   var j:integer;
   begin
    m:=k;
    if backWalk(k,l,l,j) then
     begin
      walker(k,cp,cp,l,formin);
      if (abs(xr[j])>c4*abs(xr[k])) and (abs(idx1-idx)>abs(j-idx1)) or
       ((l>=(-2+j)*3) and (l<=(j+2)*3) ) then
      begin m:=j;kuk:=#16 end
     end;
   end;{subh}

  procedure climb(var j:integer;k:integer);
   begin
    t1:=virt_vect[j];
    while (j<k) and (t1<=virt_vect[j+1]) do  {wspiac sie na rosnace}
     begin inc(j); t1:=virt_vect[j] end;     {zbocze od j do <= k}
   end;{climb}

  procedure slope(j:integer; k:integer);
    var i:integer;
    begin  kuk:='^';                 {znajdujemy sie na tylnym zboczu spktr}
     i:=j;                           {w pkcie j}
     climb(j,k);                     {wspiac sie na najbl.szczyt}
     if j=k then i:=round(count/j) else
      begin
       i:=round(count/k);
       j:=round(count/j);
       walker(j,i,i,i,backMin);
      end;
     if (abs(idx1-idx)>abs(i-idx1)) and (abs(xr[i])>=c3*abs(xr[idx]))
     then begin idx:=i; logic:=true end                         {nie wpadac}
    end;{slope}                                                  {w wysokie}

  procedure BigRise;
   begin
    kuk:=#31;
    if not forWalk(idx,k) then          {sprawdzic, czy sa subharmoniki}
     begin           {dla idx nie ma subharm. wiec wprowadzic poprawke idx}
      i:=idx;
      forMin(i,cp);
      subh(k,i,j);
      if j<>k then idx:=j
      else
       if (forWalk(k,j) or (abs(idx1-idx)>abs(k-idx1))) and
        (k<hcp) then begin idx:=k; kuk:=#174 end;
     end;
   end;{bigRise}

  BEGIN
 {
 1. POSZUKIWANIE MAKSIMUM W CEPSTRUM
 }
   logico:=false; logic:=false;
   t1:=0;
   hcp:=round(c2*rate/1000+idx);
   backDown(hcp,2);
   if hcp>round(rate/f_low) then hcp:=round(rate/f_low);
   if hcp>=cp then begin hcp:=cp-1; f_low:=rate/hcp end;
   i:=round(rate/f_high);                        {pozycja f_high w cepstrum}
   if i<2 then begin  i:=2; f_high:=rate/2 end;
   j:=round(f_high*count/rate);                  {pozycja f_high w spectrum}
   t1:=2*rate*idx1*(idx/(rate*(idx1+idx)+(2*c1*idx1)*idx));
   lcp:=round(t1);
   if lcp<i then lcp:=i;
   forDown(lcp,hcp);
   hsp:=round(count/t1);
   if hsp>j then hsp:=j;
   idx1:=idx;
   maximum(lcp,hcp,idx);           {poszukac w cepstrum wartosci najwiekszej}
 {
   2.  POROWNAC LOKALIZACJE AMPLITUDY MODULU CEPSTRUM W OBSZARZE POBUDZENIA
       Z LOKALIZACJA AMPLITUDY MODULU specTRUM
 }
   maxiS(2,hsp,k);
   j:=idx; backMin(j,round(count/k));
   j:=round(count/j);
   if j<k then slope(j,k);               {jestesmy na tylnym zboczu }
   j:=round(count/k); BackMin(j,2);            {spectrum}
   if (j>idx) and not logic then
    begin                            {slaba zawartosc harmonicznych w pobudz}
     idx:=j; j:=round(1.5*j); if j>=cp then j:=cp-1;
     walker(idx,j,j,i,forMin);
     if abs(idx1-idx)>abs(i-idx) then idx:=i
     else
      begin
       forDown(idx,hcp);
       while (abs(xr[idx])<abs(xr[idx+1])) and (idx<hcp) do inc(idx);
      end;
     if idx>count/k then slope(round(count/idx),k);
     logico:=true
    end;
   subh(idx,lcp,idx);
   if logico then kuk:=#174 else if logic then kuk:='>' else kuk:=#15;
 {
  3.  ZBYT DUZY PRZYROST F0
 }
   if (rate/idx>rate/idx1+c1) or ((idx1-idx)/idx>=c5) then BigRise;
   if (idx-idx1)/idx>=c5 then subh(idx,lcp,idx);
   if idx>=hcp then begin walker(hcp,lcp,lcp,idx,backMin); kuk:=#175 end;
 {
  4. OKRESLIC GORNA CZESTOTLIWOSC PASMA POBUDZENIA:
     Poszukac najmniejszej wartosci w module spectrum w przedziale F0, 2F0
 }
   j:=round(count/idx);                    {j, to pozycja akt. F0 w spectrum}
   i:=round(2*count/idx);                  {j:=fcpstr*count/rate}
   if i>hsp then i:=hsp;
   k:=j;
   climb(j,i);                                {wspiac sie na rosnace zbocze}
   while j<=i do                              {szukac najmn. wartosci w spectr}
    begin                                     {tylko do przodu, do 2*F0}
     if t1>=virt_vect[j] then begin t1:=virt_vect[j]; k:=j end;
     inc(j)
    end;

 {
  5.  KOSMETYKA: poszukac najblizszego minimum w cepstrum idac w dol
               nie przekraczajac granic <F0,2F0>
 }

   k:=round(count/k);                                    {pozycja w cepstrum}
   if k>=1 then
    begin
     i:=idx div 2;
     if i<lcp then i:=lcp;
     backMin(k,i);
    end;{if}                             {wycofac sie z kosmetyki, wkroczono}
   if k<=lcp then forMin(k,idx);         {na obsz. odp. imp.}
   fcpstr:=rate/k;
   F0:=rate/idx;
   fIdx:=round(count/k);
   Cpik:=abs(xr[idx]);
   for i:=0 to count-1 do virt_vect[i]:=abs(xr[i]);     {cepstrum na ekran}
  (* buzz_hiss;
   if jump and hiss then F0:=F0_1;
   if tsd then buzz_hiss_tun_print(i,t1);*)
   if tsp then buzz_hiss_files_write;  {shift-F7, wyniki analizy ton/szum na}
  END;{Auto_Tun}                       { dysk}

procedure standard_Ins(k:word);
{
Nadaje standardowa wartosc tablicy ReIns - pik na poczatku i na koncu -
a tablicy ImIns - zera.
Tablice te sa uzywane do przechowywania wstawek do cepstrum.
}
 var cynk:nameStr; l:word;
 begin
    if F0Fix<rate/count then F0Fix:=rate/count;
    F0Fix:=round(rate*round(count/rate*F0Fix)/count);
    str(F0Fix:6:1,cynk);
    informator('Laduje pobudzenie sztuczne F0:='+cynk,80,12,k,100);
    k:=round(rate/F0Fix);
    if InsPunct>k then
     begin
      InsPunct:=k;
//      sound(50);
      informator('Zmiana punktu wstawiania pobudzenia!',80,12,k,100);
//      nosound;
     end;
    l:=1;
    if k>=cp then k:=cp-1;
    ReIns:=Yzero; ImIns:=Yzero;
    repeat
     ReIns[k]:=400/l;                    {F0, pik dla 300 hz}
     ReIns[count-k]:=400/l;
     inc(l);
     k:=round(l*rate/F0Fix);
    until k>=cp;
 end;{standard_ReIns}

procedure getCpstr;
 begin
  informator('Laduje pobudzenie naturalne',80,12,255,200);
  reIns:=xr; ImIns:=yi;
 end;{getCpstr}

procedure PutCpstr;
 var k,l,m:word;
 begin
  if hiFiltr then begin l:=insPunct; m:=count-l end
  else
   begin
    for k:=count-l+1 to count-1 do
     begin
      xr[k]:=reIns[k];
      yi[k]:=ImIns[k];
     end;
    l:=0;
    if insPunct>0 then m:=InsPunct-1
    else m:=0;
   end;{else}
   for k:=l to m do
    begin
     xr[k]:=reIns[k];
     yi[k]:=ImIns[k];
    end;
 end;{PutCpstr}

procedure GetPutCpstr;
{
Zaladowac aktualne tablice xr^ i  yi^ do ReIns i do ImIns odpowiednio
}
 begin
  getCpstr;
 end;{GetPutCpstr}


procedure ChangeInsert;

  begin
   if art_nat then  GetPutCpstr
   else standard_Ins(255);
   change_insert:=false
  end;


procedure mix_phase_cepstrum;

 begin
  if Tribolett then Tribolett_Phase_Unwrap(virt_vect)
  else Apl_Phase_Unwrap(virt_vect);
  if not Change_insert then
   cpstr_gr1(cpstr,fIdx,'Faza rozwinieta');
  modul_npdbrms(xr,yi,xr,angle,count,true,false);          {modul transf. fft w neperach}
  yi:=virt_vect;                           {podstawic faze pod czesc urojona}
  if auto_freq and filtrOn then virt_vect:=xr;{zachowac modul dla automatyki}
  if prtScr and cpstr_to_scr then
   druk('mix_phase_cepstrum - vektory przed transformacja odwrotna'#10#13' '+
        'X - modul w neperach, Y, V_V - faza w radianach');
  invers(xr,yi);
 end;{mix_phase_cepstrum}


procedure min_phase_cepstrum;
{oblicza cepstrum minimalno-fazowe}
 var i:word;
 begin
  for i:=0 to count-1 do                          {obliczyc faze w radianach}
   if yi[i]=0 then virt_vect[i]:=0 else virt_vect[i]:=aTan(xr[i],yi[i]);
  modul_npdbrms(xr,yi,xr,angle,count,true,false);          {modul transf. fft w neperach}
  yi:=virt_vect;                           {podstawic faze pod czesc urojona}
  if auto_freq and filtrOn then virt_vect:=xr;{zachowac modul dla automatyki}
  if prtScr and cpstr_to_scr then
   druk('min_phase_cepstrum - vektory przed transformacja odwrotna'#10#13' '+
        'X - modul w neperach, Y, V_V - faza w radianach');
  invers(xr,yi);
 end;{min_phase_cepstrum}


procedure power_cepstrum;
{oblicza cepstrum mocy}
 begin
  modul_npdbrms(xr,yi,xr,angle,count,true,false);          {modul transf. fft w neperach}
  yi:=yZero;                                   {--------------wylaczyc faze!}
  if auto_freq and filtrOn then virt_vect:=xr;{zachowac modul dla automatyki}
  if prtScr and cpstr_to_scr then
   druk('Cepstrum min_fazowe - vektory przed transformacja odwrotna'#10#13' '+
        'X, V_V - modul w neperach, Y- faza w radianach');
  invers(xr,yi);
 end;{power_cepstrum}


procedure cepstrum;
 var logik1,logik2,logik3 : boolean;
          i : word;
 begin  //---------------------cepstrum-------------------
  logik1:=grafika and moovie and cpstr_to_scr;
  if logik1 then
  if fcpstr=0 then
   begin fcpstr:=rate/jot; cIdx:=jot;fIdx:=round(count/jot) end;
  jot:=count;
  repeat
   if PhTest then PhaseTest(xr,yi,virt_vect);
   case power_cmplx_cpstr of
    1 :  power_cepstrum;
    2 :  min_phase_cepstrum;
    3 :  mix_phase_cepstrum
    else
     informator('Blad programu, zle okreslona wersja cepstrum',0,12,255,200);;
   end;{case}

  if auto_freq and filtrOn and (time_change <>nr_konc) then
   auto_tun(f_low,f_high,cIdx,F0idx1,F0idx)
  else
   begin
    for i:=0 to count-1 do virt_vect[i]:=abs(xr[i]);      {cepstrum na ekran}
    buzz_hiss;                                              {i na dysk}
   end;
  if prtScr and cpstr_to_scr then
   case power_cmplx_cpstr of
     1 : druk('Cepstrum - vektory po transformacji odwracajacej'#10#13+
         'X - cepstrum-real, Y - cepstrum-imag, V_V - cepstrum-modul');
     2 : druk('Cepstrum - vektory po transformacji odwracajacej'#10#13+
         'X - cepstrum-real, Y - cepstrum-imag, V_V - cepstrum-modul');
     3 : druk('Cepstrum mocy - vektory po transformacji jw.'#10#13+
         'X - cepstrum-real, Y - cepstrum-imag, V_V - cepstrum-modul');
   end;{case}
  if logik1 then cpstr_gr2(cIdx,15{white},12{lightRed});
  if change_Insert then ChangeInsert
  else
   begin
    if filtrOn then filtr(xr,yi,count,fcpstr,cIdx,barL,barM);
    barJ:=cIdx
   end;
  if ins then                                        {wstawka do cepstrum}
   begin
    PutCpstr;
    if logik1 then                                   {wykres po wstawce}
     begin
      logik1:=zms; logik2:=overlap; logik3:=newCoord;
      zms:=true;   newCoord:=true;
      if krok then overlap:=false else overlap:=true;
      for i:=0 to count-1 do virt_vect[i]:=abs(xr[i]);  {cepstrum na ekran}
      cpstr_gr2(InsPunct,14{yellow},11{lightCyan});
      zms:=logik1; overlap:=logik2; newCoord:=logik3;
      if change_Insert then ChangeInsert;
     end;
   end;{ins}
  until not change_insert;
 end;{cepstrum}

procedure tran_rev(var stat : Text);
 var i:longint;
 begin
  if count_cp_bool then jot:=count else jot:=cp;
  if cpstr then
   begin                                                 {powrot do spectrum}
    if power_cmplx_cpstr=1 then   {wytworzyc sygnal minimalnofazowy}
    begin
     for i:= 1 to cp-1 do xr[i]:=2*xr[i];
     for i:= cp+1 to count-1 do xr[i]:=0;
    end;
    rfft_calc(index_mix,maxpower,sinar,cosar,xr,yi,count);
    for i:=0 to count-1 do
     begin
      xr[i]:=xr[i]/(WeightMean*count);                   {NEPERY}
      virt_vect[i]:=2*4.34294482*xr[i]+3.013;  {z neperow przechodzimy do db}
      yi[i]:=yi[i]/(WeightMean*count)                    {radiany}
     end;{for}
     virt_vect[0]:=virt_vect[0]-3.013;

    if prtScr and rev_to_scr then
     druk('TRAN_REV; Wektory po transformacji fft odwracajacej cepstrum '+
     #10#13' X - nepery, Y - radiany, V_V - decybele');
    if not Change_insert then
     if grafika and moovie and rev_to_scr then rev_gr1(freq_interval);
    if grafika and filtrOn and spectr then
     begin
      minlev:=-minlevel; maxlev:=0;
      if not Change_insert then
       kolor_spectr(scr_nr_as,freq_interval,pixY,pix1Y,pix_as,pix_as1,delta,
       interp,poz,cp,minlev,maxlev);
     end;
    if revNd and ToDiskAvail then
     //RealFileDisk(revPlik,cp,nr_konc,VminRF,VmaxRF,VsRF,VsqRF,Vdata_nb,virt_vect); blokada 19042013 z powodu braku źródła
    if retriev then xr_retriev(stat);
   end{if cpstr}
   else
    begin
     if phase_corr then
      begin
       phase_correction(xr,yi,virt_vect,zm_1,zm_2,ip,ik,firstIn,sort,fall,stat);
       lancuch:= 'bez korekcji fazy skladowych spectr.'#10#13;
       lancuch:='RETRIEV; Wektory przed transformacja '+
        'fft odwracajaca '#10#13'spectrum filtrowane pzez obciecie; '+
        lancuch+'X - real, Y - imag, V_V - RMS';
       if prtScr and ret_to_scr then druk(lancuch);
       if krok and not sing and grafika and moovie and ret_to_scr then
        if not Change_insert then ret_gr1(freq_interval);
      end
    else lancuch:='po korekcji fazy skladowych spectr.'#10#13;
     invers(xr,yi);
     if not phase_corr and (ik<count-2) then
      begin zm_1:=xr[ik+1]; zm_2:=xr[ik+2] end;
     if prtScr and rev_to_scr then
      druk('Wektory po transformacji  odwracajacej spectrum');
     if not Change_insert then
      if grafika and moovie and rev_to_scr then rev_gr2;
    end;{else}
 end;{tran_rev}

var i:byte;

BEGIN
{
 --kryteria poczatkowe dla automatycznego przestrajania filtru cepstralnego:
}
        new(circleArr);           {- bufor pierscieniowy}
          count:=256;              //określono tu, bo count=0 i zalamywała się inicjalizacja; zależy od tego parametr dOm oraz, później, phinc
       sqrCount:=count*count;
            ins:=false;           {tryb insert - wstawianie F0 do cepstrum}
          quiet:=false;           {- wlaczyc ton dla F0 i kor. fazy}
        firstIn:=true;            {- rozpoczeto korekte fazy}
      auto_freq:=false;           {-autom. przestrajanie filtru cepstralnego}
          f_low:=20;              {dolna czestotl. graniczna}
         f_high:=800;             {gorna   "          "     }
 power_cmplx_cpstr:=1;            {cepstrum mocy|rozwinieto- minimalno- fazowe}
       cut_init:=false;           {wyciac poczatkowe linie cepstrum}
       init_cut:=0;               {lba pocz. linii cpstr. do wyciecia}
        pas_odd:=false;           {wyciac linie parzyste}
       pas_even:=false;           {wyciac linie nieparzyste}
             F0:=0;               {czestotliwosc podstawowa}
            buzz:=1;              {poprzednia wartosc wsp. buzz}
         buzz_1:=1;               {aktualna     "      "    "  }
         p_buzz:=1;               {poprzednia wartosc buzz}
          p_rms:=1;               {   "                rms}
           rms1:=0;               {1-wsza wartosc do porown. (sredn. z buf.)}
            rms:=1;               {2-ga     "         "      }
           hiss:=true;            {segment bezdzwieczny}
 hiss_buzz_header:=true;          {wylaczyc druk naglowka tabeli wydrukow}
                                  {kontrolnych analiz ton/szum}
               tsd:=false;        {wydruki ton/szum}
               tsp:=false;        {zanotowac wydruki j.w. na dysku}
      buzz_to_disk:=false;        {wyniki analiz ton/szum na dysk}
               dev:=0;            {odchylka od rozkladu rownomiernego}
              jump:=5;            {- odstep probkowania spectrum w [ms]}
     buzz_rms_skok:=20;           {- odstep porownywania probek RMS'u w [ms]}
cr_m:=round(buzz_rms_skok/jump);  {- pojemnosc pierscienia opozniajacego
                                     (modyfikowana razem z modyf. skoku!!!)}
 for i:=1 to cr_m do circleArr^[i]:=rms1;
 setLength(signal,count);                      {wycinek sygnalu}
 setLength(fazR1,count);                       {pochodna spectrum}
 setLength(fazI1,count);
 setLength(fazR2,count);                       {druga pochodna spectrum}
 setLength(fazI2,count);
 dOm:=2*pi/count;
 dOmKw:=sqr(dOm);
 ph_corr_nb:=0;                   {liczba korekt fazy}
 corr_spectr_nb:=0;               {liczba korygowanych wektorow}
 sort:=false;                     {kor.fazy - sortowac indeksy wg. sily}
 fall:=false;                     {kor. fazy od slabych do silnych skladow.}
 unwrap_iter:=true;               {Tribolet, iteracje przy rozwijaniu fazy}
 thld:=pi;                        {Tribolett, prog wejscia do iteracji}
 Tribolett:=false;                {rozwijac faze metoda Triboletta|APL}
 phTest:=false;                   {generowac dane do testu rozwijania fazy}
 double_side_cpstr_filtr:=true;   {filtrowac cepstrum centralnie}
 with restat do                   {wartosci poczatkowe zbioru statystycznego}
  begin                           {dla oceny skutecznosci korekty fazy}
    zm_0:=0;
 delta_0:=0;
   s_d_0:=0;
abss_d_0:=0;
 k_s_d_0:=0;
     s_d:=0;
  abss_d:=0;
   k_s_d:=0;
  end;
    Cpik:=0;                      {amplituda cepstrum w pkcie F0}
   F0idx:=1;                      {indeks dla czestotliwosci F0}
  F0idx1:=1;                      {indeks dla poprzedniej F0}
END.


