Environmental data band-pass filter

hello everyone!

I am working with environmental data, in this case I am working with CO data.

I have to perform the kriging operation but it gives me problems because I have drift in the data I acquire via the sensors. For this reason I have to apply a band pass filter to the original data, to modify the variogram curve.

As this is the first time I am working with these things, can you give me some help with this?

if it helps I can share csv files

thank you very much

ALBERTO

This is the code:

clear;
data_gps = load('C:\Users\alber\OneDrive\Desktop\TESI\PROVE_CAMPUS\3-03 (casa)\output_gga.csv');
data_datalog = load('C:\Users\alber\OneDrive\Desktop\TESI\PROVE_CAMPUS\3-03 (casa)\DATALOG.csv');
latitudine = data_gps(:,1);
longitudine = data_gps(:,2);
altitudine = data_gps(:,3);
% converto le coord in coord locali
origine = [latitudine(1) longitudine(1) altitudine(1)];
[xEast, yNorth, zUp] = latlon2local(latitudine, longitudine, altitudine, origine);

[X, Y] = meshgrid(linspace(min(xEast), max(xEast), 166), linspace(min(yNorth), max(yNorth), 166));
n = length(data_datalog);
x = xEast; 
y = yNorth;

co = data_datalog(:,1);
z = co;

% ----------------------- BANDPASS FILTER ----------------------------
% Defining band-pass filter cut-off frequencies
lowFreq = 0.04; % Low cut-off frequency (Hz)
highFreq = 0.2; % High cut-off frequency (Hz)

% Apply bandpass filter to 'co' data
Fs = 1; % Sample rate,  I sample at 1Hz
co_filtered = bandpass(z, [lowFreq highFreq], Fs);
% ----------------------------------------------------------------------
subplot(2,2,1)
scatter(x, y, 10, co_filtered, 'filled'); axis image; axis xy
% geoscatter(lat,lon,10, z, 'filled')
% geobasemap openstreetmap
colorbar
title('track with pollution samples CO')

v = variogram([x y], co_filtered, 'plotit', false, 'maxdist', 100);
% and fit a spherical variogram
subplot(2,2,2)
[~,~,~,vstruct] = variogramfit(v.distance, v.val, [], [], [], 'model', 'stable');
xlabel('lag distance h [m]')
ylabel('\gamma(h) [ppm]')
title('variogram')

% now use the sampled locations in a kriging
[Zhat, Zvar] = kriging(vstruct, x, y, co_filtered, X, Y);
subplot(2,2,3)
imagesc(X(1,:),Y(:,1),Zhat); axis image; axis xy
title('kriging predictions')
subplot(2,2,4)
contour(X, Y, Zvar); axis image
title('kriging variance')

[Mathieu NOE]

hello

so did you get what you wanted or ?

[Alberto CIPPELLETTI]

no, unfortunately the bandpass filter does not seem to work properly.

I don't know if I'm doing something wrong with the bandpass filter or if I'm doing something else wrong.

The kriging does not work correctly.

[Taylor]

It would definitely helped if you shared the CSV files

[Alberto CIPPELLETTI]

OK, I will share the csv file with you.

The DATALOG.csv file contains several environmental data values, the CO ones are in the first column.

The output_gga.csv file contains latitude and longitude coordinates obtained with a gps.

[Mathieu NOE]

so this is how your CO data looks like before and after the bandpass filter (seems to do what it's designed for)

so how would you like your CO data look like ? I wonder if what you want is simply a baseline correction (so all the CO values would reamin positive but the lower peaks would be aligned on a y = 0 line)

clear;
data_datalog = load('DATALOG.csv');
co = data_datalog(:,1);
z = co;
% ----------------------- BANDPASS FILTER ----------------------------
% Defining band-pass filter cut-off frequencies
lowFreq = 0.04; % Low cut-off frequency (Hz)
highFreq = 0.2; % High cut-off frequency (Hz)
% Apply bandpass filter to 'co' data
Fs = 1; % Sample rate,  I sample at 1Hz
co_filtered = bandpass(z, [lowFreq highFreq], Fs);

plot(co)
hold on
plot(co_filtered);
legend('co','co filtered');

[Alberto CIPPELLETTI]

yes I saw that too.

my problem is that, even with filtered CO data, the kriging operation does not work well.

If you try all the code I uploaded, in the last part, when kriging, you will see that in its plot there is no correct extrapolation.

[Mathieu NOE]

unfortunately I cannot run the entire code as I don't have all required Toolboxes

[Alberto CIPPELLETTI]

they are not too heavy and useful toolboxes, if you can download them. you can always remove them.

Answer 1

[Mathieu NOE]

it's me again :) !

I'm just guessing, maybe this is what you wanted (baseline correction)

clc
clear;
data_datalog = load('DATALOG.csv');
co = data_datalog(:,1);
z = co;


%% Run the algorithm
Base = env_secant((1:numel(z))', z, 10, 'bottom'); % option 3 with env_secant (see function attached)
Corrected_z = z - Base;

subplot(2,1,1),plot(z)
hold on
plot(Base,'--')
hold off
legend('raw data','baseline');
subplot(2,1,2),plot(Corrected_z)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ù
function [env] = env_secant(x_data, y_data, view, side)
% Function call: env_secant(x_data, y_data, view, side)
% Calculates the top envelope of data <y_data> over <x_data>.
% Method used: 'secant-method'
% env_secant() observates the max. slope of about <view> points,
% and joints them to the resulting envelope.
% An interpolation over original x-values is done finally.
% <side> ('top' or 'bottom') defines which side to evolve.
% Author: Andreas Martin, Volkswagen AG, Germany

    if nargin == 0
        test( false );
        test( true );
        return
    end

    if nargin < 4
        error( '%s needs at least 4 input arguments!', mfilename );
    end

    assert( isnumeric(view) && isscalar(view) && view > 1, ...
           'Parameter <view> must be a value greater than 1!' );
    assert( isvector(x_data) && isnumeric(x_data) && all( isfinite( x_data ) ), ...
           'Parameter <x_data> has to be of vector type, holding finite numeric values!' );
    assert( isvector(y_data) && isnumeric(y_data) && all( isfinite( y_data ) ), ...
           'Parameter <y_data> has to be of vector type, holding finite numeric values!' );
    assert( isequal( size(x_data), size(y_data) ), ...
           'Parameters <x_data> and <y_data> must have same size and dimension!' );
    assert( ischar(side), ...
            'Parameter <side> must be ''top'' or ''bottom''!' );

    switch lower( side )
        case 'top'
            side = 1;
        case 'bottom'
            side = -1;
        otherwise
            error( 'Parameter <side> must be ''top'' or ''bottom''!' );
    end

    sz = size( x_data );

    x_data = x_data(:);
    x_diff = diff( x_data );
    x_diff = [min(x_diff), max(x_diff)];
    assert( x_diff(1) > 0, '<x_data> must be monotonic increasing!' );

    y_data = y_data(:);
    data_len = length( y_data );
    x_new = [];
    y_new = [];

    if diff( x_diff ) < eps( max(x_data) ) + eps( min(x_data) )
        % x_data is equidistant
        search_fcn = @( y_data, ii, i ) ...
                     max( ( y_data(ii) - y_data(i) ) ./ (ii-i)' * side );
    else
        % x_data is not equidistant
        search_fcn = @( y_data, ii, i ) ...
                     max( ( y_data(ii) - y_data(i) ) ./ ( x_data(ii) - x_data(i) ) * side );
    end


    i = 1;
    while i < data_len;
        ii = i+1:min( i + view, data_len );
        [ m, idx ] = search_fcn( y_data, ii, i );

        % New max. slope: store new "observation point"
        i = i + idx;
        x_new(end+1) = x_data(i);
        y_new(end+1) = y_data(i);
    end;

    env = interp1( x_new, y_new, x_data, 'linear', 'extrap' );
    env = reshape( env, sz );




function test( flagMonotonic )
    npts = 100000;
    y_data = cumsum( randn( npts, 1 ) ) .* cos( (1:npts)/50 )' + 100 * cos( (1:npts)/6000 )';
    if flagMonotonic
        x_data = (1:npts)' + 10;
    else
        x_diff = rand( size( y_data ) );
        x_data = cumsum( x_diff );
    end

    view = ceil( npts * 0.01 ); % 1 Percent of total length
    env_up = env_secant( x_data, y_data, view, 'top' );
    env_lo = env_secant( x_data, y_data, view, 'bottom' );

    figure
    plot( x_data, y_data, '-', 'Color', 0.8 * ones(3,1) );
    hold all
    h(1) = plot( x_data, env_up, 'b-', 'DisplayName', 'top' );
    h(2) = plot( x_data, env_lo, 'g-', 'DisplayName', 'bottom' );

    grid
    legend( 'show', h )
end

end

[Alberto CIPPELLETTI]

thank you very much for the help you are giving me, but I have a different problem.

If you see the code I shared, you will see that there is an operation where the variogram is calculated and then the kriging operation is done.

I have problems with this part: the kriging is not done correctly because the variogram curve grows too quickly.