Do small but non-zero coefficients in symmetric Hilbert FIR found from firpm indicate a bug?
3 views (last 30 days)
Show older comments
I noticed that depending on how you generate the coefficients of a Hilbert Transform FIR, one method (firpm using 'Hilbert') gives small non-zero (but significant at 16 bit bitwidths) coefficients and the other method does not (since they are 0 by construction). Why does the first method do this? Is it a bug? The following code indicates the behavior. The magnitude response plotted is nearly the same, but the impulse response of the first (at 16 bits) can be seen to not have desired 0 coefficients every other tap for the imaginary side. This makes it not possible to use the first method to design a filter for hardware that doesn't have these multipliers.
%%is_this_a_hilbert_bug.m
% This code shows that for shorter filters and wider transition bandwidths, the Hilbert filter design method for firpm
% can give non-zero coefficients in places they are undesired (here at wordlenght = 16 bits, all these non-zero taps are
% quantized to a number > 0, implying that multipliers are needed in hardware. However, the halfband design performs
% just as well and has no such issue. Shouldn't the firpm 'Hilbert' be deprecated (or fixed)?
clear, clc, close all
format long
%%Assign config
fs = 1; % sample frquency
Nfft = 2^12; % FFT length
f = fs*(0:Nfft - 1)/Nfft; % frequency plotting vector
transBW = 0.2; % transition bandwidth (passband to stopband)
L = 13;
%%Generate FIR from firpm using 'Hilbert' option
coeffIm = firpm(L-1, [transBW 1-transBW], [1 1], 'Hilbert').';
coeffIm = 1j*coeffIm;
coeffRe = [zeros(floor(L/2), 1); 1; zeros(floor(L/2), 1)];
h1 = coeffRe + coeffIm;
H1 = mag2db(abs(fft(h1, Nfft)));
%%Generate FIR using halfband design method and then freq shift the impulse response to get a Hilbert FIR
halfband_spec = fdesign.halfband('n,tw' , L - 1, 2*transBW);
halfband_design = design(halfband_spec,'equiripple');
h_halfband = halfband_design.Numerator.';
h_halfband = h_halfband/max(abs(h_halfband));
n = (-(L-1)/2:1:(L-1)/2).';
h2 = h_halfband.*exp(1j*pi*n/2);
H2 = mag2db(abs(fft(h2, Nfft)));
%%Display impulse response and plot magnitude response
coeff_float = [real(h1) imag(h1) real(h2) imag(h2)];
coeff_fix = round(2^15*coeff_float);
coeff_method_1_re_im = coeff_fix(:, 1:2) %#ok<NOPTS>
coeff_method_2_re_im = coeff_fix(:, 3:4) %#ok<NOPTS>
figure, plot(f, H1, 'b', f, H2, 'g')
grid on
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
title('Hilbert Transform FIR Magnitude Response')
legend('firpm Hilbert', 'freq shift halfband')
0 Comments
Answers (0)
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!