% make_hilb.m
% David Rowe May 2015
%
% creates Hilber Transformer FIR coeffs

graphics_toolkit ("gnuplot");

% from https://www.dsprelated.com/freebooks/sasp/Hilbert_Transform_Design_Example.html

M = 257;    % window length = FIR filter length (Window Method)
fs = 8000;  % sampling rate assumed (Hz)
f1 = 100;   % lower pass-band limit = transition bandwidth (Hz)
beta = 8;   % beta for Kaiser window for decent side-lobe rejection
fn = fs/2;             % Nyquist limit (Hz)
f2 = fn - f1;          % upper pass-band limit
N = 2^(nextpow2(8*M)); % large FFT for interpolated display
k1 = round(N*f1/fs);   % lower band edge in bins
if k1<2, k1=2; end;    % cannot have dc or fn response
kn = N/2 + 1;          % bin index at Nyquist limit (1-based)
k2 = kn-k1+1;          % high-frequency band edge
f1 = k1*fs/N           % quantized band-edge frequencies
f2 = k2*fs/N
w = kaiser(M,beta)'; % Kaiser window in "linear phase form"
H = [ ([0:k1-2]/(k1-1)).^8,ones(1,k2-k1+1),...
    ([k1-2:-1:0]/(k1-1)).^8, zeros(1,N/2-1)];
h = ifft(H);               % desired impulse response
hodd = imag(h(1:2:N));     % This should be zero

% put window in zero-phase form:
wzp = [w((M+1)/2:M), zeros(1,N-M), w(1:(M-1)/2)];
hw = wzp .* h;             % single-sideband FIR filter, zero-centered
Hw = fft(hw);
hh = [hw(N-(M-1)/2+1:N),hw(1:(M+1)/2)]; % causal FIR
hh *= 2;

figure(1);
HH = fft([hh,zeros(1,N-M)]);
plot(20*log10(abs(HH)));
figure(2);
subplot(211); plot(real(hh)); title('real imp resp');
subplot(212); plot(imag(hh)); title('imag imp resp');

% save coeffs to a C header file

f=fopen("../src/ht_coeff.h","wt");
fprintf(f,"/* Hilbert Transform FIR filter coeffs */\n");
fprintf(f,"/* Generated by make_hilb Octave script */\n");

fprintf(f,"\n#define HT_N %d\n\n", M);

fprintf(f,"COMP ht_coeff[]={\n");
for r=1:M
  if r < M
    fprintf(f, "  {%f,%f},\n",  real(hh(r)), imag(hh(r)));
  else
    fprintf(f, "  {%f,%f}\n};", real(hh(r)), imag(hh(r)));
  end
end

fclose(f);
