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]=maxi(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]=maxi(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)
  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 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