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')
	naglowki=(tokens(mgetl(f,1),ascii(9)))'
	l=strsubst(mgetl(f,1),',','.');
	dane=msscanf(l,'%e'+ascii(9)+'%e'+ascii(9)+'%e'+ascii(9))
	wiersz=1
	while meof(f)==0 do
		wiersz=wiersz+1
    l=strsubst(mgetl(f,1),',','.');
		if isempty(l)==%f then
		  dane(wiersz,:)=msscanf(l,'%e'+ascii(9)+'%e'+ascii(9)+'%e'+ascii(9))
    end
	end
	mclose(f)
endfunction

function [Xk,xn1,nf1]=fft_el_f1(xn,kol_t,kol_x,f1)
// na podstawie elementów od poczatku przebiegu do najwiekszej wielokrotnosci
// okresu obliczonego jako 1/f1, z kolumny numer kol_x macierzy xn,
// uzupe³niajšc próbki zerami do najbli¿szej potêgi 2
// oblicza FFT
// wynik znajduje siê w kolumnie 3 macierzy Xk
// nastêpnie na podstawie kolumny kol_t wyznacza rozdzielczo?æ widma w hercach
// i na tej podstawie umieszcza warto?ci czêstotliwo?ci w kolumnie 2 macierzy Xk
// w kolumnie 1 macierzy Xk znajdujš siê numery harmonicznych 
// (0 -- sk³adowa sta³a -- wiersz 1)
// nf1 -- numer rzeczywistej sk³adowej podstawowej (liczšc sk³adowš sta³š jako 0)

  // wyznaczenie ca³kowitej liczby okresów
	Ts=xn(2,kol_t)-xn(1,kol_t)
	N=size(xn,'r')
	T=1/f1
	nT=T/Ts // liczba próbek na okres
	lb_T=floor((N-1)/nT) // ca³kowita liczba okresów w przedziale próbkowania
	nsup=1+round(lb_T*nT) // indeks na którym nale¿y obcišæ próbki; 1 bo tak sš numerowane wiersze macierzy
	xn1=xn(1:nsup,kol_t)
	xn1(:,2)=xn(1:nsup,kol_x)
	// obliczenie FFT
	Xk=fft_el(xn1,1,2)
	nf1=lb_T
endfunction

function Xk=fft_el(xn,kol_t,kol_x)
// oblicza i przeskalowuje FFT bioršc pod uwagê ca³š macierz próbek xn
// pozosta³e dzia³anie jak fft_el_f1
	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 Xk=dft_czt(xn,kol_t,kol_x,theta,phi,L)
// na podstawie kolumny kol_x macierzy xn wyznacza L wyrazów DFT jako transformatê
// chirp-z obliczonš w punktach na okrêgu jednostkowym od kšta fazowego theta 
// do theta+(L-1)phi
  N=size(xn,'r')
  Xk1=czt(xn(:,kol_x)',L,1,phi,1,theta)
  L2=floor(L/2)
  Xk=[0:(L2-1)]'
  N=size(xn,'r')
	Ts=xn(2,kol_t)-xn(1,kol_t)
  Xk(:,2)=1/(2*%pi*Ts)*(theta+Xk(:,1)*phi)
  Xk(:,3)=Xk1(1:L2)'/(N/2)
  if theta==0 then
    Xk(1,3)=Xk1(1)/N
  end
endfunction

function [Xk,xn1,f1,nf1]=fft_autof1(xn,kol_t,kol_x)
// ustala czêstotliwo?æ sk³adowej podstawowej a nastêpnie obcina próbki do
// ca³kowitej wielokrotno?ci okresu i oblicza FFT
// xn -- macierz próbek w dziedzinie czasu
// kol_t -- numer kolumny z warto?ciami czasu
// kol_x -- numer kolumny z danymi do przetworzenia
 
  //zgrubne wyznaczenie czêstotliwo?ci sk³adowej podstawowej
  Xk1=fft_el(xn,kol_t,kol_x)
  [m,nf1]=max(abs(Xk1(:,3)),'r')
  f1=Xk1(nf1,2)
  //dok³adne wyznaczenie czêstotliwo?ci sk³adowej podstawowej
	Ts=xn(2,kol_t)-xn(1,kol_t)
	fmin=0.5*f1
	fmax=1.5*f1
	deltaf=f1/100
  theta=2*%pi*fmin*Ts
  phi=2*%pi*deltaf*Ts
  L=abs(ceil(4*%pi*(fmax-fmin)*Ts/phi))
  Xk2=dft_czt(xn,kol_t,kol_x,theta,phi,L)
  [m,nf1]=max(abs(Xk2(:,3)),'r')
  f1=Xk2(nf1,2)
  //FFT z obciêciem próbek
  [Xk,xn1,nf1]=fft_el_f1(xn,kol_t,kol_x,f1)
endfunction

function Xgk=dft_grupuj_ampl(Xk,f1,nf1)
// grupuje widmo w prš¿ki odpowiadajšce sk³adowym harmonicznym wzglêdem
// rzeczywistej pierwszej harmonicznej
// Xk -- widmo wej?ciowe
// f1 -- rzeczywista czêstotliwo?æ 1. harmonicznej
// nf1 -- numer harmonicznej w widmie Xk, która odpowiada rzeczywistej
//        1. harmonicznej
  N=size(Xk,'r')
  maxm=abs(floor((N-1)/nf1))
  Xgk=Xk(1,:)
  deltan=floor((nf1-1)/2)
  a=nf1-1
  b=2
  mod_test=a-b*int(a/b)
  if mod_test<>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 fft_wykres(Xk,lb_harm)
// wy?wietla w postaci prš¿ków modu³ kolumny 3 z Xk wzglêdem kolumny 2, do
// harmonicznej nr lb_harm w³šcznie lub wszystkie prš¿ki je¿eli lb_harm jest 0
	clf()
	if (lb_harm>0) then
		plot2d3(Xk(1:(lb_harm+1),2),abs(Xk(1:(lb_harm+1),3)))
	else
		plot2d3(Xk(:,2),abs(Xk(:,3)))
	end
	e=gce()
	e.children(1).thickness=3
	e.children(1).foreground=2
endfunction

function fft_kopiuj(Xk)
// kopiuje widmo Xk do schowka w postaci tekstowej z tabulacjš, zamieniajšc
// kropki na przecinki
  s='Numer harmonicznej'+ascii(9)+'Czêstotliwo?æ'+ascii(9)+'Amplituda'+ascii(13)+'k'+ascii(9)+'f(k)'+ascii(9)+'abs(X(k))'+ascii(13)
  for n=[1:size(Xk,'r')] do
    s=s+msprintf('%e'+ascii(9)+'%e'+ascii(9)+'%e'+ascii(13),Xk(n,1),Xk(n,2),abs(Xk(n,3)))
  end
  clipboard('copy',strsubst(s,'.',','))
endfunction

function [Iok,fs]=fftfal(dane,kol_t,kol_uio,ki)
// kompletna procedura wyznaczajšca widmo pršdu wyj?ciowego falownika Iok i 
// czêstotliwo?æ sterowania fs
// dane -- macierz próbek w dziedzinie czasu
// kol_t -- numer kolumny zawierajšcej warto?ci czasu
// kol_uio -- numer kolumny zawierajšcej warto?ci napiêcia na boczniku
// ki -- mno¿nik pozwalajšcy uzyskaæ rzeczywiste warto?ci pršdu w amperach
  ion=dane(:,kol_t)
  ion(:,2)=ki*dane(:,kol_uio)
  [Xk1,xn2,f1,nf1]=fft_autof1(ion,1,2)
  Iok=dft_grupuj_ampl(Xk1,f1,nf1)
  fs=f1
endfunction
