Testing Hilbert transforms in Octave

simulation/hilbert/end2end.m

@@ -0,0 +1,82 @@
+%% Generate Hilbert coeffecients
+order = 8;
+if 1
+  order = floor(order/2)*2; %Make sure order is even
+  %See instructions for generating poles/zeros at https://dsp.stackexchange.com/questions/8692/hilbert-transform-filter-for-audio-applications-using-iir-half-band-parallel-al
+  Ik = e.^(pi./(2.^(1:1:order)));
+  Ik = Ik.*((-1).**(1:order));
+  Isos = zp2sos([Ik], 1./[Ik], 1);
+  Qsos = zp2sos([-Ik], 1./[-Ik], 1);
+  [b a] = sos2tf(Isos);
+  h1 = freqz(b,a);
+  Isos = zp2sos([Ik], 1./[Ik], 1/max(abs(h1)));
+  Qsos = zp2sos([-Ik], 1./[-Ik], 1/max(abs(h1)));
+else
+   order = floor(order/2)*2+1; %Make sure order is odd
+   if 1 % Type III FIR
+    window = hamming(order*2+1);
+    a = (-order:order);
+    b = abs(rem(a,2)); %Remove even values
+    c = 2/pi./a; %Scale and invert
+    c(order+1)=0; % Remove NaN
+    c = c.* window' .* b;
+  else % Type IV FIR
+    window = hamming(order+1);
+    a = (-order:2:order);
+    c = 2/pi./a; %Scale and invert
+    c = c.* window';
+  end
+  Qsos = tf2sos(c,[1]);
+  % delay compensation
+  Isos = zeros(floor(order/2),6);
+  Isos(:,3:4)=1;
+end
+
+%% Create test audio
+fs = 20e3;
+start_freq = 10; %Hz
+stop_freq = 10e3; %Hz
+duration = 3;    %Seconds
+w = logspace(log10(start_freq),log10(stop_freq),fs*duration);
+%t = linspace(start_freq, stop_freq, fs*duration);
+t = sosfilt([1 0 0 1 -1 0], w); %Numerical integration
+audio1 = cos(2*pi/fs*t);
+scale = 2^11;
+audio = round(audio1 * scale)/scale;
+%audio = chirp([0:0.001:5]);
+%% Hilbert transform
+if 0
+  h = hilbert(audio);
+  I = real(h);
+  Q = imag(h);
+else
+  I = sosfilt(Isos, audio);
+  Q = sosfilt(Qsos, audio);
+end
+%% Reconstruction
+if 0
+  A = real(hilbert(I)); %To compensate for group delay
+  B = imag(hilbert(Q)); 
+  audio2 = (A-B)/2;
+elseif 1
+  A = sosfilt(Isos, I);
+  B = sosfilt(Qsos, Q);
+  audio2 = (A-B)/2;
+else
+  A = sosfilt(Qsos, I);
+  B = sosfilt(Isos, Q);
+  audio2 = (A+B)/2;
+end
+%%Playback
+figure(1);
+specgram(audio);
+p1=audioplayer(audio1,fs,24);
+play(p1);
+pause(duration);
+figure(2);
+specgram(audio2);
+p2=audioplayer(audio2,fs,8);
+play(p2);
+pause(duration);
+figure(3);
+semilogx(w,audio2);

simulation/hilbert/freq_test.m

@@ -0,0 +1,31 @@
+fs = 5e3; % Sample frequency
+len = 7; % half-filter length
+%window = flattopwin(len*2+1)';
+
+%%%% Hilbert Transform
+if 1 % Type III FIR
+  window = hamming(len*2+1);
+  a = (-len:len);
+  b = abs(rem(a,2)); %Remove even values
+  c = 2/pi./a; %Scale and invert
+  c(len+1)=0; % Remove NaN
+  c = c.* window' .* b;
+else % Type IV FIR
+  len = 2*len+1;
+  window = hamming(len+1);
+  a = (-len:2:len);
+  c = 2/pi./a; %Scale and invert
+  c = c.* window';
+end
+
+B = c;
+A = [1];
+W = logspace(-4, -0.7, 512);
+h = freqz(B, A, W);
+figure(1);
+subplot(2,1,1);
+loglog(W,abs(h));
+subplot(2,1,2);
+g = unwrap(-angle(h));
+g -= len * W; % Subtract constant group delay
+semilogx(W,g/pi*180);

simulation/hilbert/iir_hilbert.m

@@ -0,0 +1,20 @@
+order = 12;
+I = e.^(pi./(2.^(1:2:order)));
+Q = e.^(pi./(2.^(2:2:order)));
+I = [I -I];
+Q = [Q -Q];
+W = logspace(-4, -.7, 512);
+[B A] = zp2tf(I, 1./I, 1);
+h1 = freqz(B,A,W);
+h1gain = 1/abs(h1(end));
+[B A] = zp2tf(I, 1./I, h1gain);
+h1 = freqz(B,A,W);
+[B A] = zp2tf(Q, 1./Q, 1);
+h2 = freqz(B,A,W);
+h2gain = 1/abs(h2(end));
+[B A] = zp2tf(Q, 1./Q, h2gain);
+h2 = freqz(B,A,W);
+figure(1);
+semilogx(W, 180/pi.*[unwrap(angle(h1))- unwrap(angle(h2))]);
+figure(2);
+semilogx(W, abs([h1;h2]));

simulation/hilbert/test.m

@@ -0,0 +1,54 @@
+fs = 5e3; % Sample frequency
+len = 7; % half-filter length
+%window = flattopwin(len*2+1)';
+
+%%%% Hilbert Transform
+if 0 % Type III FIR
+  window = hamming(len*2+1);
+  a = (-len:len);
+  b = abs(rem(a,2)); %Remove even values
+  c = 2/pi./a; %Scale and invert
+  c(len+1)=0; % Remove NaN
+  c = c.* window' .* b;
+else % Type IV FIR
+  len = 2*len+1;
+  window = hamming(len+1);
+  a = (-len:2:len);
+  c = 2/pi./a; %Scale and invert
+  c = c.* window';
+end
+
+function mag = hilmag(test_f, fs, hilbert)
+  test_f /= fs;
+  test_len = max(100, 2/test_f + 2*length(hilbert));
+  x = sin(2*pi*test_f.*(1:test_len));
+  y = conv(x, hilbert);
+  y = y(length(hilbert)+1:end-length(hilbert)); % Only take the middle section
+  mag1 = rms(y);
+  mag2 = rms(x);
+  mag = [mag1 mag2];
+end
+
+function p = phase(test_f, fs, hilbert)
+  test_f /= fs;
+  test_len = max(100, 2/test_f + 2*length(hilbert));
+  x = sin(2*pi*test_f.*(1:test_len));
+  y = x;
+  y = conv(y, hilbert);
+  test_len = floor(test_len/2);
+  b = floor(length(hilbert));
+  amp = rms(x)/rms(y);
+  p = atan2(x(test_len), y(test_len+b));
+end
+
+f = logspace(log10(10),log10(2400),200);
+%f = linspace(10, 5000, 200);
+m = zeros(length(f),2);
+p = zeros(length(f),1);
+for i=1:length(f)
+  m(i,:)=hilmag(f(i), fs, c);
+%  p(i,:)=phase(f(i), fs, c);
+end
+loglog(f,m)
+%p = unwrap(p);
+%plot(f,p)

simulation/hilbert/test1.m

@@ -0,0 +1,52 @@
+fs = 5e3; % Sample frequency
+len = 7; % half-filter length
+window = hamming(len*2+1)';
+
+%%%% Hilbert Transform
+a = (-len:len);
+b = abs(rem(a,2)); %Remove even values
+c = 2/pi./a.*b; %Scale and invert
+c(len+1)=0; %Fix the divide-by-zero
+c = c.* window;
+figure(1);
+subplot(1,2,1)
+bar(c) % This is an FIR hilbert transform
+
+%%%% Differentiation
+d = (-1).**a;
+d = d./a;
+d(len+1) = 0; %Fix the divide by zero
+d = d.* window;
+%subplot(1,3,2)
+%bar(d)
+
+test_f = 0.01;
+test_len = 200;
+x = sin(2*pi*test_f.*(1:test_len));
+%y = sosfilt([1 0 0 1 0 -1],x);
+y = conv(x, c); % Hilbert transform
+y = y(len+1:end-len); % Only take the middle section
+x = conv(x, [-1 0 2 0 -1]); % Compensation Filter
+x = x(3:end-2);
+subplot(1,2,2)
+plot(y);
+hold on;
+plot(x);
+hold off;
+
+phase = atan2(x,y);
+phase = unwrap(phase);
+%f = conv(phase,d); % differentiate
+%f = f(len*2+1:end-len*2);
+f = phase(2:end)-phase(1:end-1);
+f *= 1/2/pi;
+figure(2);
+subplot(1,2,1)
+plot(phase);
+subplot(1,2,2)
+plot(f);
+
+##h = conv(c,d);
+##h = round(h)
+##subplot(1,3,3)
+##bar(h)