# FIR filter design with frequency sampling method - setting proper phase response for IDFT

I'd like to implement freuqency sampling method for linear phase FIR filter design using IDFT transform.

My procedure goes like this :

1. determine desired magnitude response value in frequency points
2. I add linear phase response function with group delay (N-1)/2 to get complex frequency response values and equidistant points in frequency
3. use IDFT to calculate system's impulse response.

This approach works for odd-length FIR filter but I have problems with even-length FIR filters - I just don't get proper result.

I guess I'm doing something wrong maybe with phase function or something else - but what ?

Your procedure is of course correct. It is in the details where things usually go wrong. One important thing is how to extend the desired frequency response beyond Nyquist taking the required symmetry into account. In order for the filter coefficients to be real-valued, the desired frequency response must satisfy

$$D[k]=D^*[N-k],\quad k=0,1,\ldots,N-1\tag{1}$$

where $*$ denotes complex conjugation, and $N$ is the FFT length (= filter length). From (1) it follows that for even $N$, the elements of the complex desired frequency response are

$$D_0,D_1,\ldots,D_{N/2},D^*_{N/2-1},\ldots,D^*_1$$

For odd $N$ you get

$$D_0,D_1,\ldots,D_{(N-1)/2},D^*_{(N-1)/2},\ldots,D^*_1$$ Note that in order for (1) to be satisfied, $D_0$ and, for even $N$, $D_{N/2}$ must be real-valued.

This little Matlab/Octave code works for even and odd filter lengths.

% frequency sampling design of linear phase FIR filter

N = 64;                         % FFT length = filter length
np = floor(N/2) + 1;            % number of independent frequency points
n = 0:np-1;
w = n*2*pi/N;                   % frequency vector
M = sin(n*pi/(np-1));           % some desired magnitude response
D = M.*exp(-1i*(N-1)/2*w);      % desired complex frequency response (linear phase)
D = [D,conj(D(N-np+1:-1:2))];   % append redundant points for IFFT
h = ifft(D);                    % compute impulse response
max(abs(imag(h)))               % should be very close to 0
h = real(h);                    % remove numerical inaccuracies

% check result
[H,w2] = freqz(h,1,4*N);
plot(w/2/pi,abs(D(1:np)),'.',w2/2/pi,abs(H))


• I don't understand why you are multiplying by the term $exp(-1i*(N-1)/2*w)$ when generating the desired frequency response. What's the purpose of that? – keith Oct 3 '17 at 11:24
• @keith: That's the linear phase term. If you just have a magnitude, i.e. a real-valued frequency response, the impulse response would be non-causal. – Matt L. Oct 3 '17 at 15:51
• That's odd. I'm using your code snippet to design a filter with arbitrary frequency and phase, and the $exp(-1i*(N-1)/2*w)$ makes it work. I don't understand why it's a "linear phase" term given it seems necessary for designing any arbitrary magnitude/phase FIR. Maybe I just got lucky using the code and a more complicate bit of code is needed for arbitrary magnitude/phase FIR? – keith Oct 3 '17 at 16:16
• @keith: That term basically adds some delay to the desired frequency response to make it easier for a causal filter to approximate the desired response. It all depends on the definition of your desired (non-linear) phase if that's necessary or not. If the ideal filter (with the specified desired frequency response) is non-causal, such a delay term will greatly improve the approximation by the actual causal filter. – Matt L. Oct 4 '17 at 8:06
• If things are still unclear it might make sense to formulate a new question. – Matt L. Oct 4 '17 at 8:06

I would like to suggest a variation on the procedure given by Matt L.

I think a better approach is to use a large number of frequency samples (> 1024) no matter the number of taps needed. This prevents the loss of so much frequency domain information when only a few taps are required.

If the frequency samples describe a linear phase filter, then the desired taps will be located in the center bins of the Inverse FFT output. Simply truncate the output to the desired length and apply an appropriate time domain window, such as a Kaiser.

If one works through a simple example of a windowed filter, this use of over sampling gives the same results, to several decimal places, as the commonly used sinc pulse formula.

This link describes this method of frequency over sampling in some detail. In particular, it shows how to set the frequency domain phase for the various filter types (even or odd tap count, low pass, band pass, etc).

For N being odd, h(n)=1/N[H(0)+ 2∑_(k=1)^((N-1)/2)▒〖Re{H(k) e^(j2πnk/N )}〗]

For N being even and H(N/2)=0 , h(n)=1/N [H(0)+2∑_1^(N/2-1)▒Re{H(k) e^(j2πnk/N)}]