shallow zip-file copy from codec2 e9d726bf20

1 files changed, 60 insertions, 0 deletions

diff --git a/octave/make_hilb.m b/octave/make_hilb.m
new file mode 100644
index 0000000..6c6323c
--- /dev/null
+++ b/octave/make_hilb.m
@@ -0,0 +1,60 @@
+% 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);