How can I design Nyquist interpolation filters with the Parks-McClellan algorithm? We can easily design interpolation filters that obey certain frequency-domain constraints using the Parks-McClellan algorithm.  However, it's not immediately clear how to enforce time-domain constraints; in particular, I'm interested in generating Nyquist filters.  So if I'm oversampling by a factor of N, I want the filter to have zero-crossings at kN, for non-zero integer k (this ensures that the input samples to my interpolator will appear in the output sequence).

I've seen Harris1 talk about a technique for designing half-band filters, i.e. the special case where N=2.  Is there a general solution for this?  (I know that we can easily design filters with the window method, but that doesn't give us the same control.)


[1] Multirate Signal Processing for Communication Systems, pp. 208-209
One design method, albeit one that is limited to powers of two, would be to start with one halfband filter, insert zeros at every other (creates a spectral replica), then convolve it with a second halfband filter having a wider transition band.  Repeat the process until you get to the required power of 2.

Here's an example that creates a lowpass filter with Fc=fs/8 and zero crossings every 4 samples:

b0=remez(34,[0 .45 .55 1],[1 1 0 0])';
b1=remez(6,[0 .25 .75 1],[1 1 0 0])';
b0up = zeros(1,2*length(b0)-1);
b0up(1:2:end) = b0;
B0up=freqz(b0up);
b2 = conv(b0up,b1);  % length = 34*2+1 + 6 = 75 coefficients