0

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.

rose_t
  • 93
  • 5
  • 1
    Are you trying to reverse engineer this code? – Thomas Sablik Mar 09 '21 at 22:15
  • I'm not sure what you mean. The code gives me everything I need. I'm just not sure which variables are related to the frequency response as I don't know what all of the outputs mean. – rose_t Mar 10 '21 at 00:28
  • 1
    Reverse engineering means you have an undocumented system and you want to find out how it works. How would you calculate the frequency window of Lanczos filter with pen and paper? – Thomas Sablik Mar 10 '21 at 07:48

1 Answers1

0

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.

carla.beg
  • 16
  • Hmm. I guess I'll have to ask the professor if there's some data transformation we need to do because that's exactly what I had done and it looks nothing like the graph she showed in class. Thank you though. – rose_t Mar 10 '21 at 20:33