matlab - FIR to linear-phase FIR

I am trying to convert an FIR impulse into linear-phase FIR impulse, I need to remove the phase response, so when it is applied to a sound it only changes the amplitude response.

Here is my code below:

clear all; 
close all;        
[x,fs]=audioread('impulse.wav')


s=fft(x)    
ph=angle(s)

amp=abs(s)    
s2=real(ifft(pol2cart(0,amp))) ; setting phase to zero

s2= circshift(s2,(round(length(s2)/2)))
plot(s2)


s2=s2/max(max(abs(s2)))

audiowrite('LinearPhaseImpulse.wav',s2,fs,'BitsPerSample',24)

• Is this the correct way of converting an FIR to a linear-phase FIR impulse response?


• How would one negate the phase so the FIR becomes linear-phase?

Answer

Converting an arbitrary FIR filter into a linear-phase FIR with the same magnitude response is generally impossible. As mentioned in msm's answer, this conversion must be done with an allpass filter, which must be an IIR filter, and only in those rare cases where you get appropriate pole-zero cancellations will the total filter be FIR. In general you'll end up with an IIR filter.

Having said that, there are possibilities to convert a given FIR filter into a linear-phase FIR filter with approximately the same magnitude response. By increasing the length of the linear-phase FIR filter, the approximation can be made arbitrarily good.

The method is similar to the one you implemented, but the crucial difference is zero-padding. You need to use an FFT length that is much greater than the length of the given FIR filter.

Below is an octave example explaining how it works. I chose a non-linear phase filter of length $30$, and approximated it by a length $39$ linear-phase FIR filter:

% come up with some non-linear phase FIR filter
[b,a] = butter(5,.3);
h = impz(b,a);
h = h(1:30);    % length 30
h = h(:)';


% do the FFT trick
N = 1024;       % FFT length (even!)
H = fft(h,N);
hlp = real(ifft(abs(H)));
hlp = [hlp(N/2+2:N),hlp(1:N/2)];

% windowing
Nw = 39;        % choose window length (odd!) = filter length
w = ones(1,Nw); % rect. window (choose another if you wish)
w = w(:)';

Nw2 = (Nw-1)/2;
hlp = hlp(N/2-Nw2:N/2+Nw2) .* w;

% plot
Hlp = fft(hlp,N);
n = (0:N/2)*2/N;
subplot(2,1,1), stem(h), hold on, stem(hlp,'r'), hold off
legend('NLP FIR', 'LP FIR')
subplot(2,1,2), 
plot(n,20*log10(abs(H(1:N/2+1))),n,20*log10(abs(Hlp(1:N/2+1))),'r')

legend('NLP FIR', 'LP FIR') ,grid on

The plot below shows the impulse responses of the two FIR filters, and the magnitude responses. The magnitude approximation could be further improved by choosing a longer filter length for the linear-phase FIR filter.