Strange results with all pass compensation

Question

I am new to DSP and, as a learning exercise, I am trying to design an all-pass filter in order to compensate the phase of an input impulse response (heqm.mat, which you can download from here). The MATLAB code is the following; for the phase compensation I am using the time-reversed technique, also called matched-filter:

load heqm.mat; % Load the heqm vector
HEQM      = fft(heqm);
% Design the all-pass filter
HEQMP     = exp(1j.*angle(HEQM));
tmp       = ifft(HEQMP,'symmetric');
% Time-reversed technique
eqpr      = tmp(end:-1:1);
% Figure on the left
figure; plot(unwrap(angle(fft(heqm))));
hold on; plot(unwrap(angle(fft(eqpr))));
% Convolution of the original impulse response with the all-pass filter
o = conv(heqm,eqpr);
% Figure on the right
figure; plot(unwrap(angle(fft(o))));

As you can see from the left picture, the phase of the all-pass filter is the opposite of the original one. But if I convolve the two vectors, why the resulting phase (right figure) isn't flat as it should be? What am I missing here?

Answer 1

You need to distinguish between circular and linear convolution. The DFT considers all signals to be periodic with the FFT length.

Your all pass filter is non-causal. Due to the periodic nature of the FFT the non-causal part for $t < 0$ shows up at the end of the buffer. Applying linear convolution with $conv$ makes a mess: you get two copies of your signal. One with the causal part and one with the non-causal part which is delayed by one FFT length.

The way to fix non-causalities is to apply bulk delay: Circulate the impulse response until the peak is in the middle, than convolve. Note that this will also increase the length. This method will always create some amount of error and your system will have overall linear phase behavior.

%% rotate so that the all pass becomes causal
n = length(heqm);
eqprRot = circshift(eqpr,n/2);
% linear cponvolution
result = conv(eqprRot,heqm);
% append a zero to get back to 2^n length
result = [result; 0];
% plot the phase minuse the bulk delay
fResult = fft(circshift(result,-n/2+1));
figure; clf;
plot(angle(fResult(1:n+1)));