Arbitrary Phase FIR design using FFT

Question

I want to design a digital FIR filter to cancel the phase of an analog Butterworth filter (i.e. the phase is smooth) i.e. the filter has unity gain with just an inverted phase component.

The following standard method for designing linear phase filters via FFT guarantees that the phase is linear (without deviation) even after applying a window (the window is obviously useful to prevent ripples occurring if the required impulse response is longer than the desired length): FIR filter design with frequency sampling method - setting proper phase response for IDFT

Is there a similar method where by careful maths you can ensure that the magnitude is guaranteed to be unity (even after windowing to smooth any "ripples" in the resulting phase if the impulse is longer than the desired length)?

Answer 1

You are trying to design a FIR allpass filter. That's impossible (other than a simple delay) since true allpass filters consist of zero/pole pairs that are inverse of each others and the FIR poles are always at the origin.

However, with enough resources, you can always get "close enough" for your specific application. Of course, you would have to define what "close enough" means.

The actual design is simple enough

1. Sample your filter target on a very dense FFT grid
2. Do an inverse FFT
3. Window/truncate the resulting impulse response until you meet your "close enough" requirements. Typically the impulse response will decay exponentially, so truncation often gives better results that any type windowing (for the same filter length)
4. In most cases the impulse response will be non-causal, so you will need to add bulk delay.

Below is a piece of Matlab code that does that

%% start with a Butter filter
fs = 44100; % sample rate
bwCut = 1000; % 1000 Hz cut off frequency
bwOrder = 3;  % thrid order
%% design the filter
[b,a] = butter(bwOrder,bwCut/(fs/2));

%% sample on a large FFT grid
nFFT = 16384;
delta = zeros(nFFT,1); delta(1) = 1;
bwIR = filter(b,a,delta); % impulse response
bwTF = fft(bwIR);         % transer fundction

%% design the target
targetTF = exp(-1i.*angle(bwTF)); % that works every where except for Nyquist, where
                                  % the phase is undefined, we'll manually smootht this out
targetTF(nFFT/2+1) = 1;           % Nyquist phase repair                                  
h0 = real(ifft(targetTF));        % target impulse response

%% now we have to cut it to a practical size
% let's see if we can capture all samples that are above a certain
% thershold. That's in essence a rectangular window
thresh = max(abs(h0)).*0.01;      % -40dB re max
tmp = circshift(h0,nFFT/2);       % shift to the middle of FFT window
i1 = find(abs(tmp) > thresh);     % find the area above threshold
h1 = tmp(min(i1):max(i1));        % grab that area

%% analyse the result. This really depends on the application.
% In this example we will look at the difference of the resulting
% transfer function to the "ideal" answer. The difference looks at both the
% real and imaginary parts, so it's a decent metric that combines magnitude
% and phase. It coould be interpreted as the amount of "noise" induced by
% the truncation error.
idealTF = abs(bwTF);

% model the results, rectangular window
h2 = conv(h1,bwIR);                   % convolution
h2 = h2(1:nFFT);                      % cut off the excess length
h2 = circshift(h2,min(i1)-nFFT/2-1);  % time align
diffTFRect = fft(h2)-idealTF;         % calculate difference
% same length hanning window
h1w = h1.*hanning(length(h1));
h2 = conv(h1w,bwIR);                  % convolution
h2 = h2(1:nFFT);                      % cut off the excess length
h2 = circshift(h2,min(i1)-nFFT/2-1);  % time align
diffTFHanning = fft(h2)-idealTF;


%% plot it
clf;
i = (1:nFFT/2);
fx = [diffTFRect diffTFHanning];
plot(i*fs/nFFT, 20*log10(abs(fx(1+i,:))));
xlabel('Frequency in Hz ->');
set(gca,'xlim',fs*[1/nFFT 1/2]);
set(gca,'xscale','log');
set(gca,'ylim',[-120, 10]);
ylabel('Difference in dB');
grid('on');
legend('Rectangle','Hanning');
title('Error');