How to remove a range of frequencies from a time series using the Lanczos filter?

8 views (last 30 days)
Tracka Lowd
Tracka Lowd on 8 Apr 2022
Commented: William Rose on 8 Apr 2022
I would like to use this Lanczos filter code to remove the 3-9 day time scale from a time series. My data is at a 12-hourly resolution. So I computed the cut-off frequencies for the low-pass and high-pass filters as:
time_series = .5 + 1.5*rand(700,1) % dummy time series
Cf_lowpass = 3/(12/24)
Cf_highpass = 9/(12/24)
From here I can apply the filter (both low and high-pass). However, I am not sure if that is the right way to compute the cut-off frequencies. Also, after getting the low-pass and high-pass filtered time series, do I add or subtract them in order to remove the 3-9 time scale? Thanks.
  1 Comment
William Rose
William Rose on 8 Apr 2022
The frequencis I computed above are normalized the way Matlab's filters expect: by the Nyquist frequency. However, I looked at LanczosFilter.m, on the File Exchange. One of the inputs to the function is the sampling interval, dT. The author says: "Units of the cut-off frequency, [Cf], must be [dT]^{-1}". My initial interpretation of this is that if dT=0.01 seconds (Fs=100 Hz), then the units for Cf should be 1/(0.01 seconds). Therefore, if you wanted a filter with a 20 Hz cutoff, you would use Cf=20 Hz/(1/.01sec)= 20 Hz/100 Hz = 0.2. However, my reading of the input argument error-checking code at the top of LanczosFilter.m suggests that this is not the correct interpretation. When you use a funciton from the FileExchange, test it carefully to be sure you understand it fully. In this case, I would generate simulated signals and see if LanczosFilter.m removes the frequencies as expected.
My reading of LanczosFilter.m is that you should enter dT in units of hours or days, and you should enter the cutoff frequency in units of inverse hours or inverse days, respectively.
If you choose hours, dT=12 (units=hours) and Cf1=1/72 (units=1/hours), and Cf2=1/216.
If you choose days, dT=0.5 (units=days) and Cf1=1/3 (units=1/days), and Cf2=1/9.
That is my interpretation of the LanczosFilter.m code.

Sign in to comment.

Answers (1)

William Rose
William Rose on 8 Apr 2022
I recommend using the day as the unit of time. A three day time scale corresponds to a frequency of (1/3) per day, etc. Then you have:
fs=2; %sampling rate (units=1/day)
fco1=1/3; %cutoff frequency 1 (units=1/day)
fco2=1/9; %cutoff frequency 2 (units=1/day)
Compute normalized frequencies (normalize by the Nyquist frequency):
fn1=fco1/(fs/2); %normalized cutoff frequency 1
fn2=fco2/(fs/2); %normalized cutoff frequency 2
That should work.
  2 Comments
William Rose
William Rose on 8 Apr 2022
@Tracka Lowd, Please see my comment above, which I wrote before I saw your comment. My initial answer was based on the way cutoff frequencies are specified for filters that are part of Matlab. Then I viewed LanczosFilter.m, and I realized the author (a) does it differently than the Mathworks, and (b) does not do it the way he says he does it. My recommendations at the end of my comment above are still my recommendations. But don't take my word for it - test it yourself, to confirm that it is doing what you expect.

Sign in to comment.

Products


Release

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by