If you have the DSP System Toolbox, you can design a Comb Filter. All you need do is specify your normalised frequencies for the harmonics.
The Signal Processing Toolbox has dfilt.parallel that will do essentially the same operation, but you have to build the various filters yourself.
You can automate the filter design process in a loop. This is archive code for a different problem, but you can likely adapt it easily to your harmonics problem. The plot in figure(1) plots the passbands/stopbands:
Fs = 8200; % Samping Frequency (Hz)
Fn = Fs/2; % Nyquist Frequency
pf = linspace(20,4000,17); % Passband Frequencies
cf = pf(1:end-1)+(pf(2)-pf(1))/2; % Centre Frequencies
for k1 = 1:length(cf)
[z(k1,:),p(k1,:),k(k1)] = butter(7, [pf(k1) pf(k1+1)]/Fn);
[sos{k1},g{k1}] = zp2sos(z(k1,:),p(k1,:),k(k1));
[h(k1,:),w(k1,:)] = freqz(sos{k1},512,Fs);
end
figure(1)
freqz(sos{1})
hold on
for k1 = 2:16
freqz(sos{k1})
end
hold off
