seabirdFilter

SEABIRDFILTER  Applies a low-pass filter to an input signal

  Syntax:
    outputSignal = seabirdFilter(inputSignal, timeConstant, samplingPeriod)

  Description:
    The filter runs a low-pass filter on an input signal. A low-pass filter
    smoothes high frequency (rapidly changing) data. To produce zero phase
    (no time shift), the filter is first run forward through the data and then
    run backward through the data, so that no delays are caused by the filter.
    This filter is described in SeaBird Data Processing Manual (page 91).

  Inputs:
    inputSignal    - Signal to be filtered
    timeConstant   - Duration of the filter (in s)
    samplingPeriod - Sampling period of the input signal (in s)

  Outputs:
    outputSignal - inputSignal filtered with desired parameters

  Example:
    dataSize = 500;
    samplingPeriod = 2;
    time = (0:(dataSize-1))' * samplingPeriod;
    originalSignal = 100 * (1 + sin(linspace(0, 2 * pi, dataSize)'));
    someNoise      =  10 * (rand(dataSize, 1) - 0.5);
    inputSignal    = originalSignal + someNoise;
    outputSignal   = seabirdFilter(inputSignal, 4, samplingPeriod);
    figure(1); clf;
    plotHandle = ...
       plot(time, originalSignal, 'k', ...
            time, inputSignal,    'b', ...
            time, outputSignal,   'r');
    legend(plotHandle, 'original signal', 'noisy signal', 'filtered signal');
    xlabel('Time (secs)');
    ylabel('Depth (meters)');
    title('Glider depth');

  See also:
    BUTTER
    BUTTORD
    CHEBY1
    FILTER

  Other m-files required: none
  Subfunctions: none
  MAT-files required: none
  Dec 2010; Last revision: 23-Dec-2010

 Authors:
   Bartolome Garau  <tomeu.garau@socib.es>

0001 function outputSignal = seabirdFilter(inputSignal, timeConstant, samplingPeriod)
0002 %SEABIRDFILTER  Applies a low-pass filter to an input signal
0003 %
0004 %  Syntax:
0005 %    outputSignal = seabirdFilter(inputSignal, timeConstant, samplingPeriod)
0006 %
0007 %  Description:
0008 %    The filter runs a low-pass filter on an input signal. A low-pass filter
0009 %    smoothes high frequency (rapidly changing) data. To produce zero phase
0010 %    (no time shift), the filter is first run forward through the data and then
0011 %    run backward through the data, so that no delays are caused by the filter.
0012 %    This filter is described in SeaBird Data Processing Manual (page 91).
0013 %
0014 %  Inputs:
0015 %    inputSignal    - Signal to be filtered
0016 %    timeConstant   - Duration of the filter (in s)
0017 %    samplingPeriod - Sampling period of the input signal (in s)
0018 %
0019 %  Outputs:
0020 %    outputSignal - inputSignal filtered with desired parameters
0021 %
0022 %  Example:
0023 %    dataSize = 500;
0024 %    samplingPeriod = 2;
0025 %    time = (0:(dataSize-1))' * samplingPeriod;
0026 %    originalSignal = 100 * (1 + sin(linspace(0, 2 * pi, dataSize)'));
0027 %    someNoise      =  10 * (rand(dataSize, 1) - 0.5);
0028 %    inputSignal    = originalSignal + someNoise;
0029 %    outputSignal   = seabirdFilter(inputSignal, 4, samplingPeriod);
0030 %    figure(1); clf;
0031 %    plotHandle = ...
0032 %       plot(time, originalSignal, 'k', ...
0033 %            time, inputSignal,    'b', ...
0034 %            time, outputSignal,   'r');
0035 %    legend(plotHandle, 'original signal', 'noisy signal', 'filtered signal');
0036 %    xlabel('Time (secs)');
0037 %    ylabel('Depth (meters)');
0038 %    title('Glider depth');
0039 %
0040 %  See also:
0041 %    BUTTER
0042 %    BUTTORD
0043 %    CHEBY1
0044 %    FILTER
0045 %
0046 %  Other m-files required: none
0047 %  Subfunctions: none
0048 %  MAT-files required: none
0049 %  Dec 2010; Last revision: 23-Dec-2010
0050 %
0051 % Authors:
0052 %   Bartolome Garau  <tomeu.garau@socib.es>
0053 
0054 %  Copyright (C) 2013-2016
0055 %  ICTS SOCIB - Servei d'observacio i prediccio costaner de les Illes Balears
0056 %  <http://www.socib.es>
0057 %
0058 %  This program is free software: you can redistribute it and/or modify
0059 %  it under the terms of the GNU General Public License as published by
0060 %  the Free Software Foundation, either version 3 of the License, or
0061 %  (at your option) any later version.
0062 %
0063 %  This program is distributed in the hope that it will be useful,
0064 %  but WITHOUT ANY WARRANTY; without even the implied warranty of
0065 %  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0066 %  GNU General Public License for more details.
0067 %
0068 %  You should have received a copy of the GNU General Public License
0069 %  along with this program.  If not, see <http://www.gnu.org/licenses/>."
0070 
0071   % Store original size of input signal and
0072   % make sure it is a column vector (later will be flipped ud)
0073   origSize    = size(inputSignal);
0074   inputSignal = inputSignal(:);
0075 
0076   % Precompute some filter constants
0077   magicNumber = 2 * timeConstant / samplingPeriod;
0078   A = 1 / ( 1 + magicNumber);
0079   B = (1 - magicNumber) * A;
0080 
0081   % Pass the filter twice, first forward and then backward
0082   % to produce zero phase and thus avoid delays
0083   theSignal = inputSignal;
0084   for filterPass = 1:2
0085     outputSignal = theSignal;
0086     % Loop through the scans, recursive filter
0087     for scanIdx = 2:length(theSignal)
0088         outputSignal(scanIdx) = ...
0089             A .* (theSignal(scanIdx) + theSignal(scanIdx - 1)) -...
0090             B .* outputSignal(scanIdx - 1);
0091     end
0092     % Make the signal to be filtered the output
0093     % from the first filter pass, but turned upside down
0094     theSignal = flipud(outputSignal);
0095   end
0096 
0097   % The result from the second filter pass is in theSignal,
0098   % so reshape it to be same size as input signal
0099   outputSignal = reshape(theSignal, origSize);
0100 
0101 end