Environmental data band-pass filter
1 view (last 30 days)
Show older comments
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')
8 Comments
Mathieu NOE
on 2 Apr 2024
unfortunately I cannot run the entire code as I don't have all required Toolboxes
Answers (1)
Mathieu NOE
on 2 Apr 2024
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
See Also
Categories
Find more on Digital Filtering in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!