# wandelt ein Frequenzgang-Messfile in ein für LTspice anwendbares .subckt um.
# In REW kann ein .frd Messfile erstellt werden wie folgt:
# -> Export -> Export measurement as text -> 
# than: Use custom resolution: for example 96 PPO 
#  © juergen@micka.de 2017 

pkg load signal; % signal package for hilbert function 

# mit Minimalphasenberechnung
# gewünschte Anzahl Datenpunkte für das zu berechnende File:
Datenpunkte=1000;

FRD_Messfile="test"; %REW Messfile  ohne .frd Endung


fid = fopen([FRD_Messfile ".frd"]);
i=1;
j=1;
while(! feof(fid))
txt = fgetl (fid);
if (findstr(txt,"*"))
  zeilen_kommentar=j;
  kommentar(j++)={txt};
else
   zeilen_messdaten=i;
   [x(i), y(i), z(i)] = strread (txt, "%f %f %f");
   i++;
endif
endwhile
# C(i) = textscan(fid, '%s %*[^\n]');
fclose(fid);

startfrequenz=x(1)
endfrequenz=x(end)

# Anfang minimale Phase berechnen:
flin=linspace(startfrequenz,endfrequenz,1000000); # z.B. 20000 Frequenzpunkte (linear)
ylin=interp1(x,y,flin);
ylinlin=10.^(ylin./20);
phase=hilbert(log(abs(ylinlin)));
min_phase = unwrap(-imag(hilbert(log(abs([ylinlin, fliplr(ylinlin(2:end))])))), pi);
min_phase2 = unwrap(-imag(hilbert(log(abs([ylinlin(1:end), fliplr(ylinlin(1:end))])))), pi);
mp=180*min_phase2(1:numel(flin))/pi;
# Ende minimale Phase berechnen

figure(1);
semilogx(flin(2:end), 180 * min_phase(2:numel(flin)) / pi, flin(1:end), 180 * min_phase2(1:numel(flin)) / pi);
title("min phase response in degrees");

# Anfang Messpunkte verringern
fa=log10(startfrequenz);
fe=log10(endfrequenz);
flog=logspace(fa,fe,Datenpunkte);
ylog=interp1(flin,ylin,flog);
mp_log=interp1(flin,mp,flog);
# Ende Messpunkte verringern

zeilen= length(flog);  % Zeilenanzahl
fid = fopen ([FRD_Messfile "_frd_mp_ltspice.mod"], "w");
for i=1:zeilen_kommentar
  fprintf(fid,"%s\n",char(kommentar(i)))
endfor
fprintf(fid,[".subckt " FRD_Messfile "_frd_mp_ltspice 1 2 3 4\n"]);
fprintf(fid ,"B100 3 4 V=V(1,2) FREQ=\n");

for i = 1:zeilen-1
#  fprintf(fid ,"+(%f, %f, %f)\n",flin(i),ylin(i),mp(i));
  fprintf(fid ,"+(%f, %f, %f)\n",flog(i),ylog(i),mp_log(i));
endfor


fprintf(fid,".ends\n");
fclose (fid);



