How can I design Nyquist interpolation filters with the Parks-McClellan algorithm?

Question

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 Harris¹ 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}

Answer 1

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

Answer 2

One method to get your desired zero crossings is to do a hybrid design.

Start out with a Parks-McLellan/Remez half-band filter given equal weight to passband and stopband. Since it is a halfband filter, it will have zeros at alternate samples. You can then interpolate the time domain by sin(x)/x by zero-stuffing in the frequency domain.

Example: creating a fs/12 lowpass filter with zero crossings every 6 samples.

% prototype Remez filter 
taps=18; 
b = remez(taps,[0 .4 .6 1],[1 1 0 0])';  
% force halfband condition of zeros at every other sample
b(2:2:end)=0;  b(taps/2+1)=.5; 

% zero pad the time domain to give the Gibbs ripple some deadspace
B=fft(b,4*(taps+1) ); 
% split the frequency domain into two halves, split the Nyquist bin
Blo = [ B(1:length(B)/2) 0.5*B(length(B)/2+1) ]; 
Bhi = [ 0.5*B(length(B)/2+1) B(length(B)/2+2:length(B))  ]; 

% insert padding at pi to increase size 3x
Bpad = [ Blo zeros(1,3*length(B)-length(Blo)-length(Bhi) ) Bhi];  
bint = real( ifft(Bpad) ); % this has zeros every 6 samples

The resulting filter is close to, but not quite as good as, the prototype in terms of stopband/passband ripple. The sin(x)/x interpolation does introduce some low level ringing. You might need to slightly overdesign the prototype filter to get the required level of attenuation in the interpolated filter.