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

Question

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 ?

Thanks in advance, regards, Bul.

Answer 1

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

Answer 2

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).

Answer 3

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)}]

Answer 4

Your code logic may be right. I have the same problem with you, at last I found my C code was wrong, for example ,in the C code I use int type to divide 2, however it should use 2.0 or 2.0f not 2.