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?

    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

    Reply
    1. Patrick Martineau Post author

      Hi 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

      Reply

Leave a Reply

Your email address will not be published. Required fields are marked *