Using lanczos low pass filter in R program

Question

I am wondering if there is any package which allows us to use the Lanczos filter. I found other filters such as butterworth but I am looking for Lanczos low pass filter.

How different is Lanczos filter from butterworth filter ? Any suggestions or hints is appreciated.

Thanks.

Please help us help you by providing us with a reproducible example (i.e. code and example data), see http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example for details. — Paul Hiemstra, Jun 23 '13 at 18:57

score 5 · Accepted Answer · answered Jun 23 '13 at 19:45

Using the web I find this MATLAB implementation.

If you skipped the first part(arguments check), it looks simple to write its R equivalent.

#      Cf   - Cut-off frequency       (default: half Nyquist)
#      M    - Number of coefficients  (default: 100)
lanczos_filter_coef <- function(Cf,M=100){
  lowpass_cosine_filter_coef <- function(Cf,M)
    coef <- Cf*c(1,sin(pi*seq(M)*Cf)/(pi*seq(M)*Cf))
  hkcs <- lowpass_cosine_filter_coef(Cf,M)
  sigma <- c(1,sin(pi*seq(M)/M)/(pi*seq(M)/M))
  hkB <- hkcs*sigma
  hkA <- -hkB
  hkA[1] <- hkA[1]+1
  coef <- cbind(hkB, hkA)
  coef
}

To test it for example:

dT <- 1
Nf <- 1/(2*dT)
Cf <- Nf/2
Cf <- Cf/Nf
lanczos_filter_coef(Cf,5)

               hkB           hkA
[1,]  5.000000e-01  5.000000e-01
[2,]  2.977755e-01 -2.977755e-01
[3,]  1.475072e-17 -1.475072e-17
[4,] -5.353454e-02  5.353454e-02
[5,] -4.558222e-18  4.558222e-18
[6,]  2.481571e-18 -2.481571e-18

PS I don't know very well MATLAB(used it many years ago), so I I used this link For the R/MATLAB analogy. I hope that someone with more R/MATLAB/Scilab knowledge can test my code.

Nice find! I recommend the OP run a comparison test of your `R` code vs. native Matlab to see if the results agree. Good homework for him to do :-) — Carl Witthoft, Jun 24 '13 at 12:24

KISWENDSIDA · Answer 2 · 2018-12-24T11:14:23.533

I used the method provided in this link https://www.atmos.umd.edu/~ekalnay/syllabi/AOSC630/METO630ClassNotes13.pdf and wrote this function: `

lanczos_weights<-function(window=101,sampl_rate=1,type="lowpass",low_freq=1/100,high_freq=1/10){
low_freq=sampl_rate*low_freq
high_freq=sampl_rate*high_freq

if (type=="lowpass"){
    order = ((window - 1) %/% 2 ) + 1
    nwts = 2 * order + 1
    fc=low_freq
    w = seq(0,0,length=nwts)
    n = nwts %/% 2
    w[n+1] = 2 * fc
    k = seq(1, n-1)
    sigma = sin(pi * k / n) * n / (pi * k)
    firstfactor = sin(2 *pi * fc * k) / (pi * k)
    w[n:2] = firstfactor * sigma
    w[(n+2):(length(w)-1)] = firstfactor * sigma
    w=w[-c(1,length(w))]}
else if (type=="highpass"){
    order = ((window - 1) %/% 2 ) + 1
    nwts = 2 * order + 1
    fc=high_freq
    w = seq(0,0,length=nwts)
    n = nwts %/% 2
    w[n+1] = 2 * fc
    k = seq(1, n-1)
    sigma = sin(pi * k / n) * n / (pi * k)
    firstfactor = sin(2 *pi * fc * k) / (pi * k)
    w[n:2] = firstfactor * sigma
    w[(n+2):(length(w)-1)] = firstfactor * sigma
    w=w[-c(1,length(w))]
    w=-w
    w[order]=1-2*fc }
else if (type=="bandpass"){
    order = ((window - 1) %/% 2 ) + 1
    nwts = 2 * order + 1
    fc=low_freq
    w = seq(0,0,length=nwts)
    n = nwts %/% 2
    w[n+1] = 2 * fc
    k = seq(1, n-1)
    sigma = sin(pi * k / n) * n / (pi * k)
    firstfactor = sin(2 *pi * fc * k) / (pi * k)
    w[n:2] = firstfactor * sigma
    w[(n+2):(length(w)-1)] = firstfactor * sigma
    w1=w[-c(1,length(w))]

    order = ((window - 1) %/% 2 ) + 1
    nwts = 2 * order + 1
    fc=high_freq
    w = seq(0,0,length=nwts)
    n = nwts %/% 2
    w[n+1] = 2 * fc
    k = seq(1, n-1)
    sigma = sin(pi * k / n) * n / (pi * k)
    firstfactor = sin(2 *pi * fc * k) / (pi * k)
    w[n:2] = firstfactor * sigma
    w[(n+2):(length(w)-1)] = firstfactor * sigma
    w2=w[-c(1,length(w))]
    w=w2-w1}
else {print("Please specify a valid filter type: either 'lowpass', 'highpass' or 'bandpass'")}
return(w)}

`

#### the inputs are:
#### window: Filter length=number of weights. Corresponds to the total number of points to be lost. Should be odd: window=2N-1. The formula for N is taken from Poan et al. (2013)
#### sampl_rate: sampling rate=number of observation per time unit. ( eg: if time unit is one day, hourly data have sampl_rate=1/24)
#### type= one of "lowpass", "highpass" and "bandpass"
#### low_freq: the lowest frequency 
#### high_freq: the highest frequency

I have compared my weights to those obtained using NCL filwgts_lanczos and they are exactly the same.

Using lanczos low pass filter in R program

2 Answers2