FIR filter design with frequency sampling method - setting proper phase response for IDFT
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))