Design a comb notch filter when sampling rate is not an integer multiplicative of the cutoff frequency
You can design that filter manually without problems. Matlab just uses a very simplistic approach to comb filtering with a delay line.

In order to keeps things as simple as possible I would recommend you use a series of notch filters to remove each partial of your harmonic noise separately. That also gives you more control over how much of each harmonic really has to be removed.

A simple notch filter is one with unit gain at DC, i.e. $H(1) = 1$, and the zeros on the unit circle at the desired notch frequency, positive and negative. The poles have to be inside the unit circle, close to the zeros to cancel their effect farther away from the notch frequency. We can design this directly in the z-domain.

$$H(z) = \frac{(z-\exp(i\omega))(z-\exp(-i\omega)}{(z-r\exp(i\omega))(z-r\exp(-i\omega))}\times\frac{(1-r\exp(i\omega))(1-r\exp(-i\omega)}{(1-\exp(i\omega))(1-\exp(-i\omega))}$$

Here $\omega$ is the normalized frequency, meaning $\omega=2 \pi \frac{f_0}{f_s}$ where $f_0$ is your notch frequency and $f_s$ the sampling frequency. The parameter $r\in]0,1[$ controls the width of the notch. The second fraction makes sure we get a unit response at $\omega=0$.

To get the coefficients for a recursive filter implementation, we have to simplify the transfer function and cancel $z^2$ to turn it into a rational function of $z^{-1}$. The result is 
$$H(z) = \frac{\left(1+r^2-2r\cos(\omega)\right)\left(z^{-2}+1-2z^{-1}\cos(\omega)\right)\csc(\omega/2)^2}{4\left(r^2 z^{-2}+1-2rz^{-1}\cos(\omega) \right)}$$

So you can read off your filter coefficients to be

$$ A[n]=\left(1,-2r\cos\omega,r^2 \right)$$
$$ B[n]=\left(1,-2\cos\omega,1 \right)\cdot\frac{1}{4}\left(1+r^2-2r\cos\omega\right)\cdot\csc(\frac{\omega}{2})^2$$

The last open question is how exactly the parameter $r$ controls the notch behavior. Near the zero the transfer function on the unit circle locally behaves like

$$f(x) = \frac{x}{x-i(1-r)}$$

around $x=0$. For $r=1$ we get a singularity at $x=0$ and an ill-defined transfer function. For values $r&lt;1$ but close to $1$ we can see that for both $x\to-\infty$ and $x\to\infty$ the transfer function becomes $1$. The interesting question is, how far does the effective influence of the pole reach on the x-axis?

The square magnitude of $f(x)$ is $$|f(x)|^2=f(x)\cdot f(x)^*=\frac{x^2}{x^2+(1-r)^2}$$
which is always real and positive. We can easily equate it to a threshold and solve for $x$.The most common threshold for filter cutoffs in signal processing is $1/\sqrt{2}$, and if we solve $|f(\Delta x/2)|^2==1/\sqrt{2}$ we get the effective bandwidth of the notch to be
$$\Delta x = 2(1-r)\sqrt{1+\sqrt{2}}$$

and we can use the equivalence of our problem to the original pole/zero placement to conclude

$$ r = 1-\frac{\Delta \omega}{2\sqrt{1+\sqrt{2}}} $$ 

with the natural frequency bandwidth $\Delta \omega$. If you want this as a true frequency then again $\Delta \omega=2\pi\frac{\Delta f}{f_s}$

Note that this bandwidth calculation used a local approximation, so it is strictly only true for $(1-r)&lt;&lt;1$, or small bandwidths compared to the sampling frequency.

So this gives you a notch filter for one single harmonic of your harmonic noise. Just apply one filter at a time, one after the other, with the parameters matched to the harmonic you would like to cancel, up to fs/2. The matlab filter() function readily takes the coefficients A[n] and B[n] given above.



An attempt to implement this in scilab seems to work nicely.  Plot of before and after spectrum:



Code:

// 12740
fs = 1000;
f = 101.89798;
omega = 2*%pi*f/fs;
delta_omega = omega/2;
N= 1000;
t = [0:N];
f = [0:N]/N*fs-fs/2

r = 1 - delta_omega/(2*sqrt(1+sqrt(2)));

num = [1 -2*cos(omega) 1]*1/4*(1+r.^2-2*r*cos(omega))*csc(omega/2).^2;
den = [ 1 -2*r*cos(omega) r.^2];

X = sin(omega/2*t) + sin(omega*t) + sin(2*omega*t);
Y = filter(num,den,X);

figure(1)
clf;
plot(f,log(fftshift(abs(fft(X)))+0.01))
plot(f,log(fftshift(abs(fft(Y)))+0.01),'k:')