I am trying to apply a high-pass filter to a signal (column or row vector) consisting of 1-pixel-wide lines taken from a black-and-white image. I know the resolution of the image (res
in the code below, given in mm/pixel). How can I filter these line data in MATLAB to discard certain low frequencies (waviness) or large wavelengths, say >10 mm, using a Butterworth filter or any other?
Line data are not centered at zero.
Fs = 1; % I do not know if this assumption is correct for the image.
Fn = Fs/2; % Nyquist frequency.
lambda = 10; % Cut-off wavelength in mm, given.
samples_in_lambda = lambda/res; % divide by resolution to get samples.
fc = 1/samples_in_lambda; % Cut-off frequency from lambda.
I tried : [z, p, k] = butter(9, fc/fn, 'high'); % I see the filter is high pass on plotting.
Can I filter the line data using the above given and assumed values? If not, is there a way that I can filter the data using a cut-off wavelength?
The highest linear spatial frequency you can represent without aliasing is 1 wave cycle per 2 pixels. This means a spatial Nyquist frequency of 1 wave cycle per
2*(res*1e-3)
meters, or1000/(res*2)
reciprocal meters. (Confront this with temporal frequencies, which are measured in reciprocal seconds a.k.a. hertz).In terms of wavelengths: the shortest wave you can represent without aliasing is 2 pixels long per wave cycle. This means a spatial "Nyquist wavelength" of
res*2e-3
meters. (Confront this with temporal "wavelengths" a.k.a. periods, which are measured in seconds.)If you want to set a cutoff wavelength of 10 mm, that corresponds to a spatial frequency of 100 reciprocal meters. Since the
butter()
function takes as its second input argument (Wn
, the cutoff frequency) an arbitrary fraction of the (spatial) Nyquist frequency (the MATLAB documentation calls it "half the sampling rate"), you merely need to setWn=100/(1000/(res*2))
, i.e.Wn=res/5
.Even though your definition of the spatial sampling frequency is not quite correct (unless you are intentionally measuring it in reciprocal pixels), your final result ended up being equivalent to
Wn=res/5
, so you should be fine using the call tobutter()
that you indicated.