## Purpose

This matlab code can be useful in the atmospheric sciences but also for other geosciences. It can act as a low-pass, high-pass and band-pass filter. For example, it can be used to filter time-series to remove undesirable frequencies. The main advantage is that it can be easily applied to multi-dimensional arrays. It is based on the Discrete Fourier Transform (DFT)

## Code

function [reconstructed] = fft_filter(data,dimension,bandlimits) %The filter is designed to filter gridded data in an array of up to 5 dimensions (DATA) %with the filtering done only in the dimension specified by DIMENSION. %The signal in data is first decomposed in the frequency domain (amplitude vs frequency). %Then the amplitudes associated with frequencies outside of the BANDLIMITS are set to zero. %Then an inverse Fourier transform is performed resulting in the reconstructed signal. % DATA : array to filter of up to 5 dimensions, the dimension to filter must be evenly spaced and be a multiple of 2. % if the dimension is spaced unevenly, interpolation to regular interval is advised. % DIMENSION : 1 integer indication the dimension to filter % BANDLIMITS : 2 numbers indicating the boundary of the selected frequency band. The bandlimit is actually % expressed as the period. Bandlimits have the same unit as the spacing between measurements % ex: if filtering with measurements every day, bandlimit of [5 10] means only oscillation with periods % of 5 to 10 days are kept in the signal. To construct a low-pass or high-pass, on can use [5 Inf] or [0 5] % RECONSTRUCTED : array of same size as data containing the filtered signal. % ATTENTION : Fourier transform are accurate for periodic signals. If the filtered values are not periodic, % it is advised to make them periodic. Ex signal= 1 2 5 coult be modified to 1 2 5 5 2 1 for periodicity. % Only the first half of reconstructed should be used in this case. %Bring dimension to filter to dimension 1 si_data=size(data); order=1:length(si_data); otherdim=order(ismember(order,dimension)==0); neworder=[dimension,otherdim]; data=permute(data,[dimension,otherdim]); %Fourier transform and associated frequencies amplitudes = fft(data,[],1); n1=size(data); n=size(amplitudes,1); %replaced dimension by 1 nyquist=1/2; nyquist_location=n/2+1; %Generate a list of the frequencies if rem(size(data,1),2)==0; frequencies=[0,(1:n/2-1)/n,n/2/n,(n/2-1:-1:1)/n]; end if rem(size(data,1),2)==1; frequencies=[0,(1:(n-1)/2)/n,((n-1)/2:-1:1)/n]; end %Computes periods and decide of the ones to keep periods=1./frequencies; periodstokeep=(periods>=bandlimits(1) & periods<=bandlimits(2)); periods(periodstokeep); filteredamps=amplitudes*0; %Only keep desired periods (amplitudes) if ndims(data)==1 filteredamps(periodstokeep)=amplitudes(periodstokeep); end if ndims(data)==2 filteredamps(periodstokeep,:)=amplitudes(periodstokeep,:); end if ndims(data)==3 filteredamps(periodstokeep,:,:)=amplitudes(periodstokeep,:,:); end if ndims(data)==4 filteredamps(periodstokeep,:,:,:)=amplitudes(periodstokeep,:,:,:); end if ndims(data)==5 filteredamps(periodstokeep,:,:,:,:)=amplitudes(periodstokeep,:,:,:,:); end %Reconstruct the signal reconstructed=ifft(filteredamps,[],1); reconstructedr=real(reconstructed); reconstructedi=imag(reconstructed); %Return matrix in original order ordereturn=zeros([1,length(si_data)]); ordereturn(dimension)=1; iszero=find(ordereturn==0); for i=1:length(iszero) ordereturn(iszero(i))=i+1; end reconstructed=permute(reconstructed,ordereturn); end |

Josh WalstonHello Patrick,

My name is Josh and I am a graduate student at the University of Alaska Fairbanks. I have 55 years of daily model output that I am attempting to filter and I found your matlab function fft_filter(data,dimension,bandlimits) yesterday. First, I would like to say thank you for making this publicly available, I think it is a great function and saves me a lot of time! Second, I was wondering if you could help me understand something. I would like to use this function as a band pass filter where [7 90] day oscillations are kept in the signal. I am attempting this filter on both reference height temperature and sea level pressure vectors that are [1,20446] in length and when I do this, my filtered vector is phase shifted by one whole period. In addition the filtered output looks more convoluted than the original. Do you know why this is happening?

Thank you for your help in advance,

~Josh

—

Joshua M. Walston

Graduate Research Assistant

University of Alaska Fairbanks

Department of Atmospheric Sciences

International Arctic Research Center 408 C2

Office: (907) 474-2681

Patrick MartineauPost authorHi Josh,

It is not easy for me to understand where the problem comes from but I may have a hint. You have not mentioned if the data is periodic. Fourier transforms work for periodic functions. If the data is not periodic, ie doesn’t start and end with the same value, I think one could pad the edges with 0s so it becomes periodic (I never tried that). I usually duplicate the data to make it periodic. E.g. if my time series is 1 2 3 4 5 6 7 which is not periodic, I will do the filtering on 1 2 3 4 5 6 7 6 5 4 3 2 1 which is periodic. I then take only the first half of the filtered signal. I obtain satisfactory results with this method. I have never observed a phase shift, what do you mean by “whole period”?

Patrick