function [naglowki,dane]=wczytaj_sp107(nazwa_pliku)
// wczytuje plik TAB zapisany w programie SP107
// naglowki -- wektor łańcuchów z nagłówkami kolumn
// dane -- macierz próbek
	f=mopen(nazwa_pliku,'r')
	sep=ascii(9)
	naglowki=(tokens(mgetl(f,1),sep))'
	lb_kolumn = size(naglowki,'c')
	format_linii = ''
	for i=1:lb_kolumn do
        format_linii=format_linii+'%e'+sep
	end
	l=strsubst(mgetl(f,1),',','.')
	dane=msscanf(l,format_linii)
	wiersz=1
	while meof(f)==0 do
		wiersz=wiersz+1
    l=strsubst(mgetl(f,1),',','.');
		if isempty(l)==%f then
            dane(wiersz,:)=msscanf(l,format_linii)
    end
	end
	mclose(f)
endfunction

function [Xk,xn1,nf1]=fft_el_f1(xn,kol_t,kol_x,f1)
	Ts=xn(2,kol_t)-xn(1,kol_t)
	N=size(xn,'r')
	T=1/f1
	nT=T/Ts
	lb_T=floor((N-1)/nT)
	nsup=1+round(lb_T*nT)
	xn1=xn(1:nsup,kol_t)
	xn1(:,2)=xn(1:nsup,kol_x)
	Xk=fft_el(xn1,1,2)
	nf1=lb_T
endfunction

function Xk=fft_el(xn,kol_t,kol_x)
	Ts=xn(2,kol_t)-xn(1,kol_t)
	Xk1=fft(xn(:,kol_x),-1)
	N=size(xn,'r')
	N2=floor(N/2)
	Xk=[0:(N2-1)]'
	delta_f=1/(N*Ts)
	Xk(:,2)=[0:delta_f:((N/2-1)*delta_f)]'
	Xk(:,3)=Xk1(1:N2)/(N/2)
	Xk(1,3)=Xk1(1)/N
endfunction

function rms=rmsIS_fft(Xk,harm_inf,harm_sup)
// wyznacza wartość skuteczną biorąc pod uwagę harmoniczne od harm_inf do harm_sup
  suma = 0
  for k=(harm_inf+1):(harm_sup+1) do
    suma=suma+abs(Xk(k,3))^2/2
  end
  rms=sqrt(suma)
endfunction

function rms=rms_fft(Xk)
// wyznacza wartość skuteczną sygnału
  rms=rmsIS_fft(Xk,0,size(Xk,'r')-1)
endfunction

function rms=rmsAC_fft(Xk)
// wyznacza wartość skuteczną składowej przemiennej sygnału
  rms=rmsIS_fft(Xk,1,size(Xk,'r')-1)
endfunction

function rms=xrms_f1(xn,kol_t,kol_x,kx,f1)
  [Xk,xn1,nf1]=fft_el_f1(xn,kol_t,kol_x,f1)
  rms=rms_fft(Xk)*kx
endfunction

function rms=xrmsAC_f1(xn,kol_t,kol_x,kx,f1)
  [Xk,xn1,nf1]=fft_el_f1(xn,kol_t,kol_x,f1)
  rms=rmsAC_fft(Xk)*kx
endfunction

function pk=xpk(xn,kol_x,kx)
  pk=max(xn(:,kol_x))*kx
endfunction

function pk=xpkAV(xn,kol_x,kx)
  pk=( max(xn(:,kol_x)) - min(xn(:,kol_x)) )/2 * kx
endfunction

function av=xav_f1(xn,kol_t,kol_x,kx,f1)
  [Xk,xn1,nf1]=fft_el_f1(xn,kol_t,kol_x,f1)
  av=Xk(1,3)*kx
endfunction

function av=xavabs_f1(xn,kol_t,kol_x,kx,f1)
  xn(:,kol_x)=abs(xn(:,kol_x))
  [Xk,xn1,nf1]=fft_el_f1(xn,kol_t,kol_x,f1)
  av=Xk(1,3)*kx
endfunction

//function [Ikg,wykres1,wykres2,osie]=widmo_ii_en(in,kol_t,kol_x,ki,Pwe,st)
function dhi40=widmo_ii_en(in,kol_t,kol_x,ki,Pwe)
  [Ik,in1,nf1]=fft_el_f1(in,kol_t,kol_x,50)
  Ik(:,3)=Ik(:,3)*ki
  Ikg1=widmo_grupuj(Ik,50,nf1)
  Ikg40=[]
  en=[0,0,0,3.4,0,1.9,0,1,0,0.5,0,0.35,0]
  for k=[13:2:39] do
    en=[en,3.85/k,0]
  end
  for k=[3:2:39] do
    Ikg40=[Ikg40;k,Ikg1(1+k,3)/1e-3/Pwe,en(1+k)]
  end
  Ikg=Ikg1(1:41,:)
  clf
  plot(Ikg40(:,1),Ikg40(:,2))
	e=gce()
	wykres1=e.children(1)
	wykres1.polyline_style=3
	wykres1.thickness=3
  plot(Ikg40(:,1),Ikg40(:,3))
	e=gce()
	wykres2=e.children(1)
	wykres2.foreground=5
	osie=get('current_axes')
	osie.sub_ticks=[4,1]
  xtitle('','k','|Ii(k)|/Pi [mA/W]')
  dhi40=dh1(Ikg,40)
endfunction

function Xgk=widmo_grupuj(Xk,f1,nf1)
  N=size(Xk,'r')
  maxm=abs(floor((N-1)/nf1))
  Xgk=zeros(maxm,3)
  Xgk(1,:)=Xk(1,:)
  deltan=floor((nf1-1)/2)
  if modulo(nf1-1,2)<>0 then
    sumuj_polowki=%T
  else
    sumuj_polowki=%F
  end
  for m=[1:maxm] do
    ncen=1+m*nf1
    if m==1 then
      ninf=1+1
      nsup=ncen+deltan
    elseif m==maxm then
      ninf=ncen-deltan
      nsup=1+N-1
    else
      ninf=ncen-deltan
      nsup=ncen+deltan
    end
    suma=0
    for n=[ninf:nsup] do
      suma=suma+abs(Xk(n,3))^2
    end
    if sumuj_polowki==%T then
      if m>1 then
        suma=suma+0.5*abs(Xk(abs(ncen-deltan-1),3))^2
      end
      if m<maxm then
        suma=suma+0.5*abs(Xk(abs(ncen+deltan+1),3))^2
      end
    end
    Xgk(m+1,:)=[m,m*f1,sqrt(suma)]
  end
endfunction

function dh=dh1(Xk,lb_harm)
// wyznacza THD z odniesieniem do pierwszej harmonicznej biorąc pod uwagę
// pierwszych lb_harm (licząc od 1), lub wszystkie harmoniczne jeżeli
// lb_harm jest 0
  if lb_harm==0 then
    lb_harm=size(Xk,'r')-1
  end
  rms_bez_1 = rms_fft_infsup(Xk,2,lb_harm)
  rms_1 = Xkrms(Xk,1)
//  disp(rms_1)
//  suma_harm = 0
//  for k=(2+1):(lb_harm+1) do
//    suma_harm=suma_harm+abs(Xk(k,3))^2
//  end
//  dh=sqrt(suma_harm)/abs(Xk(1+1,3))
  dh=rms_bez_1/rms_1
endfunction

function rms=rms_fft_infsup(Xk,harm_inf,harm_sup)
// wyznacza wartość skuteczną biorąc pod uwagę harmoniczne od harm_inf do harm_sup
  suma = 0
  for k=(harm_inf+1):(harm_sup+1) do
    if k==1 then
      suma=suma+abs(Xk(k,3))^2
    else
      suma=suma+abs(Xk(k,3))^2/2
    end
  end
  rms=sqrt(suma)
endfunction

function Xkrms=Xkrms(Xk,k)
  Xkrms=abs(Xk(k+1,3))/sqrt(2)
endfunction

function [lambda,phi1]=lambda_phi(in,kol_t,kol_u,kol_i,ki)
  [Ik,in1,nf1]=fft_el_f1(in,kol_t,kol_i,50)
  Ik(:,3)=Ik(:,3)*ki
  [Uk,in1,nf1]=fft_el_f1(in,kol_t,kol_u,50)
  Urms=rms_fft(Uk)
  Irms=rms_fft(Ik)
  P=p_fft(Uk,Ik)
  S=Irms*Urms
  lambda=P/S
  phi1=(atan(imag(Uk(nf1+1,3)),real(Uk(nf1+1,3)))-atan(imag(Ik(nf1+1,3)),real(Ik(nf1+1,3))))*180/3.1416
endfunction

function p=p_fft(Uk,Ik)
	N=size(Uk,'r')
  Uk_rms=zeros(N)
	Uk_rms(1)=Uk(1,3)
	Uk_rms(2:N)=abs(Uk(2:N,3))/sqrt(2)
  Ik_rms=zeros(N)
	Ik_rms(1)=Ik(1,3)
	Ik_rms(2:N)=abs(Ik(2:N,3))/sqrt(2)
  phik=zeros(N)
//	phik(1)=0
	phik(2:N)=atan(imag(Uk(2:N,3)),real(Uk(2:N,3)))-atan(imag(Ik(2:N,3)),real(Ik(2:N,3)))
  pk=zeros(N)
	pk=Uk_rms.*Ik_rms.*cos(phik)
	p=sum(pk)
endfunction

