Использование фильтра низких частот Ланцоша в программе R

Question

Мне интересно, есть ли какой-нибудь пакет, который позволяет нам использовать фильтр Lanczos. Я нашел другие фильтры, такие как Butterworth, но я ищу фильтр низких частот Lanczos.

Насколько отличается фильтр Ланцоша от фильтра Баттерворта? Любые предложения или советы приветствуются.

Благодарю.

Answer 1

Используя Интернет, я нахожу эту реализацию MATLAB.

Если вы пропустили первую часть (проверка аргументов), выглядит просто написать ее 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

PS Я не очень хорошо знаю MATLAB(использовал его много лет назад), поэтому я использовал эту ссылку для аналогии с R/MATLAB. Я надеюсь, что кто-то с большим знанием R/MATLAB/Scilab сможет проверить мой код.

Answer 2

Я использовал метод, представленный в этой ссылке https://www.atmos.umd.edu/~ekalnay/syllabi/AOSC630/METO630ClassNotes13.pdf и написал эту функцию (больше на основе периодов, а не частот) lanczos_weights=function(N=50,dt=1,highper=90,lowper=10,filter="lowpass"){ n=-N:N W_low=ifelse(n!=0,(sin(2*pi*dt*n/lowper)*sin(pi*n/N))/(pi*n*pi*n/N),2*pi*dt/lowper) W_high=ifelse(n!=0,(sin(2*pi*dt*n/highper)*sin(pi*n/N))/(pi*n*pi*n/N),2*pi*dt/highper) W_high[ceiling(length(W_high)/2)]=1-W_high[ceiling(length(W_high)/2)];W_high[-ceiling(length(W_high)/2)]=-W_high[-ceiling(length(W_high)/2)] W_band=W_low-W_high if(filter=="lowpass"){return(W_low)}else if(filter=="highpass"){return(W_high)} else if (filter=="bandpass"){return(W_band)}else {print("Error: please specify the filter type")} } Входные данные: N(для количества весов), dt(интервал выборки), верхний и нижний периоды (максимальный и минимальный периоды) и фильтр (тип фильтра). Это работает очень хорошо, за исключением того, что, когда я сравниваю его с NCL filwgts_lanczos Я не получаю одинаковые значения.