How to Get Frequency Window from Lanczos Filter

438 views Asked by At

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.

1

There are 1 answers

1
carla.beg On BEST ANSWER

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.