I have this code for a Lanczos filter:
dT = 1 % sampling interval
Cf = 1/40 % cutoff frequency
fl = 100 % ?
M = 100 % number of coefficients ? not sure about number
LoH = 1 % low pass
Nf=1/(2*dT); %Nyquist frequency
% Normalize the cut off frequency with the Nyquist frequency:
Cf = Cf/Nf;
% Lanczos cosine coeficients:
coef = lanczos_filter_coef(Cf,M); coef = coef(:,LoH);
% Filter in frequency space:
[window,Ff] = spectral_window(coef,length(vel)); Ff = Ff*Nf;
% Filtering:
[y,Cx] = spectral_filtering(vel,window);
function coef = lanczos_filter_coef(Cf,M)
% Positive coeficients of Lanczos [low high]-pass.
hkcs = lowpass_cosine_filter_coef(Cf,M);
sigma = [1 sin(pi*(1:M)/M)./(pi*(1:M)/M)];
hkB = hkcs.*sigma;
hkA = -hkB; hkA(1) = hkA(1)+1;
coef = [hkB(:) hkA(:)];
end
function coef = lowpass_cosine_filter_coef(Cf,M)
% Positive coeficients of cosine filter low-pass.
coef = Cf*[1 sin(pi*(1:M)*Cf)./(pi*(1:M)*Cf)];
end
function [window,Ff] = spectral_window(coef,N)
% Window of cosine filter in frequency space.
Ff = 0:2/N:1; window = zeros(length(Ff),1);
for i = 1:length(Ff)
window(i) = coef(1) + 2.*sum(coef(2:end).*cos((1:length(coef)-1)'*pi*Ff(i)));
end
end
function [y,Cx] = spectral_filtering(vel,window)
% Filtering in frequency space is multiplication, (convolution in time
% space).
Nx = length(vel);
Cx = fft(vel(:)); Cx = Cx(1:floor(Nx/2)+1);
CxH = Cx.*window(:);
CxH(length(CxH)+1:Nx) = conj(CxH(Nx-length(CxH)+1:-1:2));
y = real(ifft(CxH));
end
I need to graph both the time window and the frequency window from this. I have gotten the time window, which comes from coef, but I cannot figure out which of the outputs would give me the frequency window. I've plotted every possible combination of the output variables and tried doing Fourier transform on some of them and nothing I try is giving me the expected figure.
The frequency window is in the output "window", so would need to plot(Ff,window). You get a graph with a strong drop around Cf (so 1/50), which is the cutoff frequency you choose that separates between the frequencies you are going to use for a low-pass filter VS the one for the high-pass.