Arbitrary Phase FIR design using FFT

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');