%% Hilbert-Algo-Test



pkg load signal


% Hilbert Trafo für Controller-Implementation nach: https://www.mikrocontroller.net/topic/480404#7103779
% sample_in : das datensignal links am Eingang
% delay_out : Der aktuelle Ausgang des Verzögerers (das kann nicht hier in der Funktion gemacht werden)
% coeff     : Der Koeffizient
function [sample_out,delay_in]=block2_calc(sample_in,delay_out,coeff)

  signal_top_right = coeff * ( delay_out - sample_in );
  delay_in         = signal_top_right + sample_in;
  sample_out       = signal_top_right + delay_out;
  
endfunction


% Hilbert Trafo für Controller-Implementation nach: https://www.mikrocontroller.net/topic/480404#7103779
% mit reellen Koeffizienten (wie im Bild) 

function sample_z=hilbert_algo_1(sample_in)

  % delay storage
  persistent i1=0;           %  inphase Z-1 
  persistent i2a=0,i2b=0;    %  inphase 0.3900 block
  persistent i3a=0,i3b=0;    %  inphase 0.8906 block

  persistent q1a=0,q1b=0;    %  quadratur 0.1206 block
  persistent q2a=0,q2b=0;    %  quadratur 0.6628 block

  % inphase
  
  % der inphase 0.3900 Teil
  [i2_out,delay_in]=block2_calc(i1,i2b,0.3900);
  i2b=i2a;
  i2a=delay_in; 
  % der inphase 0.8906 Teil
  [sample_out_i,delay_in]=block2_calc(i2_out,i3b,0.8906);
  i3b=i3a;
  i3a=delay_in; 
  % der inphase Z-1 Teil
  i1=sample_in;
  
  % quadratur  

  % der quadratur 0.1206 Teil
  [q1_out,delay_in]=block2_calc(sample_in,q1b,0.1206);
  q1b=q1a;
  q1a=delay_in; 

  % der quadratur 0.6628 Teil
  [sample_out_q,delay_in]=block2_calc(q1_out,q2b,0.6628);
  q2b=q2a;
  q2a=delay_in; 
  
  
  
  sample_z=sample_out_i+i*sample_out_q;
  
endfunction









%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                                                                                                         %%%
%%%                                               Main                                                      %%%
%%%                                                                                                         %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




filename='wav/zahlen1_8kHz.wav';  
[data,samplerate]=audioread(filename);
data=data';             % Spaltenvektor draus machen
%data=data(1:65000);     % Ausschnitt zum testen ---------------
max_amplitude=max(abs(data));   % Amplitude auf 1 normieren
data=data/max_amplitude;
%samplerate=8000;    
    
datalen=length(data);
n=1:datalen;
freq_axis=0.001*samplerate*n/datalen;
freq_axis=freq_axis-0.001*samplerate/2;   % Frequenzachse in kHz



% eingebaute Hilbert-Trafo

data_z_intern=hilbert(data);
fft_data_z_intern=fft(data_z_intern);




%% eigene Hilbert-Trafo berechnen

data_z=zeros(1,datalen);

for n=1:datalen,   %  im Streaming-Stil wie es im controller laufen würde

  data_z(n)=hilbert_algo_1(data(n));
  
  if bitand(n,65535)==0   % Fortschrittsanzeige
    n
  endif


end;


fft_data_z=fft(data_z);


%%%%%%%%%%%%%%%%%%%%%% Anzeige


figure(101)
semilogy(freq_axis,abs(fftshift(fft_data_z))+1e-6,freq_axis,abs(fftshift(fft_data_z_intern))+1e-6)  % verschiedene Hilbert-Algo vergleichen
%xlabel('Freq [kHz]');
title('FFT der Hilbert-Trafo')


figure(102)  % XY-Daten nach der Hilbert-Trafo
plot(data_z)