Environmental data band-pass filter
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
Alberto CIPPELLETTI
el 30 de Mzo. de 2024
Comentada: Alberto CIPPELLETTI
el 2 de Abr. de 2024
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 comentarios
Mathieu NOE
el 2 de Abr. de 2024
unfortunately I cannot run the entire code as I don't have all required Toolboxes
Respuestas (1)
Mathieu NOE
el 2 de Abr. de 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
Ver también
Categorías
Más información sobre Digital Filtering en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!