% lanczos_example.m

% demonstrates the application of a Lanczos filter to a time series

% written by Steven Greybush for AOSC 630 at UMD

% =========================================================================



nsteps=100; % number of data points in the time series

x=1:nsteps;

% generate a sample time series

% (the user can replace this with their own variation)

series=sin(2*x)+sin(x/4)+sin(x/8)+rand(1,nsteps); % the original time series

series_filtered = zeros(length(series),1); % the filtered time series



% =========================================================================

% compute Lanczos Coefficients



delta_t = 40/60; % time step between data points (hours)

tau_crit = 3; % critical period (hours)

theta_crit = 2*pi*delta_t/tau_crit; % critical frequency



M=8; % number of time steps before and after available in filtering window

h = zeros(2*M+1,1); % initialize the weights



for n=(-1*M):1:M % get the lanczos coefficients

    if n == 0

        h(M+n+1) = theta_crit*delta_t/pi;

    else

        % from Kalnay (2003)

          h(M+n+1) = (sin(n*theta_crit*delta_t)/(pi*n)) * (sin(n*pi/M)/(n*pi/M));

        % from Lynch and Huang (1992)

        % h(M+n+1) = (sin(n*theta_crit*delta_t)/(pi*n)) * (sin(n*pi/(M+1))/(n*pi/(M+1)));

    end

end



% need to adjust the weights so that they add up to 1

h_unbiased = h/sum(h);

h=h_unbiased;



% =========================================================================



% step through each time for which we'd like to create a filtered value

for t = 1+M:length(series)-M

    for n=(-1*M):1:M

        % apply the cofficients to create the filtered series

        series_filtered(t) = series_filtered(t) + series(t+n)*h(M+n+1);

    end

end



% =========================================================================

% plot the series (original and filtered)

figure(1)

hold on

plot((1:length(series))*delta_t,series)

plot((M+1:length(series)-M)*delta_t,series_filtered(M+1:length(series)-M),'r')

set(gca,'xtick',[0:6:nsteps*delta_t])

xlabel('time (hours)')

ylabel('Amplitude')

legend('original time series','filtered time series');

title('Application of Lanczos Filter to Time Series');



% plot the response function

weights = h;

figure(2)

hold on

nu_deltat = [0:pi/100:pi];

for i=1:length(nu_deltat)

    R(i) = weights(M+1);

    for n=1:1:M

        R(i) = R(i) + 2 * weights(M+n+1) * cos(n*nu_deltat(i));

    end

end

plot(nu_deltat,R)

set(gca,'xtick',[0:pi/4:pi])

ylabel('Response')

xlabel('nu * deltaT')

title('Response Function of Lanczos Filter');
