# Filtering data with matlab

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

## 2 thoughts on “Filtering data with matlab”

1. Josh Walston

Hello 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?

~Josh

Joshua M. Walston