在R程序中使用lanczos低通滤波器

Question

我想知道有没有什么包可以让我们使用Lanczos过滤器。我找到了其他滤波器，如巴特沃斯，但我正在寻找Lanczos低通滤波器。

Lanczos过滤器和巴特沃斯过滤器有什么不同？如有任何建议或提示，欢迎光临。

谢谢。

Answer 1

使用web，我找到了这个MATLAB implementation.

如果您跳过了第一部分(参数检查)，那么编写它的R等效项看起来很简单。

#      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
}

要测试它，请执行以下操作：

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

我不太了解MATLAB(很多年前就用过了)，所以我用this link来做R/MATLAB的类比。我希望有更多R/MATLAB/Scilab知识的人可以测试我的代码。

Answer 2

我使用了这个链接https://www.atmos.umd.edu/~ekalnay/syllabi/AOSC630/METO630ClassNotes13.pdf中提供的方法，并编写了这个函数：`

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

我将我的权重与使用NCL filwgts_lanczos获得的权重进行了比较，它们是完全相同的。