1% make_hilb.m 2% David Rowe May 2015 3% 4% creates Hilber Transformer FIR coeffs 5 6graphics_toolkit ("gnuplot"); 7 8% from https://www.dsprelated.com/freebooks/sasp/Hilbert_Transform_Design_Example.html 9 10M = 257; % window length = FIR filter length (Window Method) 11fs = 8000; % sampling rate assumed (Hz) 12f1 = 100; % lower pass-band limit = transition bandwidth (Hz) 13beta = 8; % beta for Kaiser window for decent side-lobe rejection 14fn = fs/2; % Nyquist limit (Hz) 15f2 = fn - f1; % upper pass-band limit 16N = 2^(nextpow2(8*M)); % large FFT for interpolated display 17k1 = round(N*f1/fs); % lower band edge in bins 18if k1<2, k1=2; end; % cannot have dc or fn response 19kn = N/2 + 1; % bin index at Nyquist limit (1-based) 20k2 = kn-k1+1; % high-frequency band edge 21f1 = k1*fs/N % quantized band-edge frequencies 22f2 = k2*fs/N 23w = kaiser(M,beta)'; % Kaiser window in "linear phase form" 24H = [ ([0:k1-2]/(k1-1)).^8,ones(1,k2-k1+1),... 25 ([k1-2:-1:0]/(k1-1)).^8, zeros(1,N/2-1)]; 26h = ifft(H); % desired impulse response 27hodd = imag(h(1:2:N)); % This should be zero 28 29% put window in zero-phase form: 30wzp = [w((M+1)/2:M), zeros(1,N-M), w(1:(M-1)/2)]; 31hw = wzp .* h; % single-sideband FIR filter, zero-centered 32Hw = fft(hw); 33hh = [hw(N-(M-1)/2+1:N),hw(1:(M+1)/2)]; % causal FIR 34hh *= 2; 35 36figure(1); 37HH = fft([hh,zeros(1,N-M)]); 38plot(20*log10(abs(HH))); 39figure(2); 40subplot(211); plot(real(hh)); title('real imp resp'); 41subplot(212); plot(imag(hh)); title('imag imp resp'); 42 43% save coeffs to a C header file 44 45f=fopen("../src/ht_coeff.h","wt"); 46fprintf(f,"/* Hilbert Transform FIR filter coeffs */\n"); 47fprintf(f,"/* Generated by make_hilb Octave script */\n"); 48 49fprintf(f,"\n#define HT_N %d\n\n", M); 50 51fprintf(f,"COMP ht_coeff[]={\n"); 52for r=1:M 53 if r < M 54 fprintf(f, " {%f,%f},\n", real(hh(r)), imag(hh(r))); 55 else 56 fprintf(f, " {%f,%f}\n};", real(hh(r)), imag(hh(r))); 57 end 58end 59 60fclose(f); 61