One way to perform this is to decimate the signal prior to filtering. If your application permits a reduction in sample rate, then this is a very nice way to go.
Without decimation, you are looking to achieve a lowpass filter with a normalized cutoff of:
fc = 60 * 2/10000 = 0.012 (very small!)
Suppose you were going to design an FIR filter using REMEZ. With a 10 kHz sample rate, a cutoff of 60 Hz, a stopband at 75 Hz, and some typical ripple specifications, the filter order, N, is approximately:
N = remezord([60 75],[1 0],[.01 .1],10000)
N =
898
If, for example, you had a sample rate of 150 Hz, your filter order would be approximately:
N = remezord([60 75],[1 0],[.01 .1],150)
N =
13
You can use RESAMPLE to achieve the reduction in sample rates.
Here's an example:
% ------------------------------------------
% Load a test signal
% ------------------------------------------
% This has a rate of 7418 Hz, but we'll pretend it's at 10kHz anyway:
load mtlb; x=mtlb;
% ------------------------------------------
% Filter and decimate:
% ------------------------------------------
f1 = 10000; % initial rate
f2 = 150; % final rate
% Remove common factors:
c = gcd(f1,f2);
p = f2/c; % Upsample
q = f1/c; % Downsample
% Simplest approach - use single stage resampling
% Resample signal from f1 to f2 Hz:
xr = resample(x,p,q);
% 60Hz LPF filter
b = fir1(32, 60*2/f2); % Or whatever you'd like
y = filter(b,1,xr);
% ------------------------------------------
% All the rest is plotting!
% ------------------------------------------
% Setup time vectors:
t1 = (0:length(x)-1)/10000;
t2 = (0:length(y)-1)/150;
% Plot results:
subplot(2,1,1);
plot(t1,x, t2,y);
title('Decimated and Filtered')
% Compensated for filter delay: "17" is the starting index for the "remainder"
% of the data, which means I've chopped off the first 16 samples.
% Now, the filter order is 32, meaning the filter length is N=33 taps. For a
% linear phase FIR filter, the delay in samples is (N-1)/2, or exactly 16 samples.
% Hence, I cut off the first 16 samples.
% Note that an even length (odd order) filter would have a fractional amount of
% delay, so I couldn't easily illustrate that case. Compensation can be achieved,
% it's just not as trivial due to the required interpolations.
subplot(2,1,2);
plot(t1,x, t2,[y(17:end);zeros(16,1)]);
title('Compensated for filter delay')
You could use the FFT method which is also known as "frequency sampling". What you would do is take the FFT, remove any components above 60 Hz and then perform the IFFT.
However, this is not recommened. The FFT method, when implemented correctly is equivalent to a particular FIR filter whose frequency response is exactly 1 at the DFT coefficients 60 Hz and below, and exactly 0 at the DFT coefficients greater than 60 Hz. However, if you look at the DTFT (on a continuum of frequency points from 0 to half the sampling frequency), the frequency sampling filter's response deviates significantly from 1 below 60 Hz, and from 0 above, everywhere BESIDES those N points where it is constrained exactly.
So, usually, it is better to use a technique which minimizes the error over a continuous range of frequencies rather than at discrete points. Frequency sampling is a simplistic filter design method which is not too effective. It's unfortunate that it is the ONLY method suggested in Numerical Recipes.
Note the issue here is that the passband is extremely narrow.
