1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
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);
|
