Homec4science
Diffusion yearII_reports (master) Browse

rSTICAZZI/labo/A4/battement.m
851e56311b98master

Last Change
Raw File

battement.m
View Options

file = 'rawdata/force4.8_f90.dat';
[data, N] = load_stable(file, 8, 50);
t = data(:,1);
f = data(:,3);
[thp, hp] = real_peaks(t, f, 3);
[tlp, lp] = real_peaks(t, -f, 3);
lp = -lp;
ob = omega_b(thp, hp);
O = Omega_factor(thp, 20, 0);
omega = 4.02;
fprintf('file: %s, omega_b = %g, Omega = %g, |omega - O|=%g\n', file, ob, O, abs(omega - O));
figure
plot(t, f, 'k+', thp, hp, 'r*', tlp, lp, 'g*')
xlabel('t')
ylabel(sprintf('V'))
grid on
%% UTILS %%
function o = Omega_factor(tp, first, last)
n = length(tp) - last - (first-1);
diff = zeros(n-1, 1);
T = 0;
for i=1:(n-1)
diff(i) = tp((first-1) + i+1) - tp((first-1) + i);
T = T + diff(i);
end
T = T / (n-1);
o = 2 * pi / T;
end
function o = omega_b(tp, hp)
n = length(hp);
s = zeros(n, 1);
for i=1:n
s(i) = hp(i);
end
[~, first] = min(s);
s(first) = [];
[~, second] = min(s);
o = 2 * pi / (tp(first) - tp(second));
end
function [data, N] = load_stable(file, lowfilter, highfilter)
data = load(file);
N = length(data(:,1));
% filter out zeros
zero = N;
while abs(data(zero,1)) < 1e-14 || data(zero, 2) > highfilter
zero = zero - 1;
end
% filter out initial domain
filt = 1;
while data(filt,2) < lowfilter
filt = filt + 1;
end
data = data(filt:zero,2:end);
N = length(data);
end
function [T, M] = real_peaks(t, f, a)
smoothed = smooth(f);
[~, inds] = findpeaks(smoothed);
s = length(inds);
T = zeros(s, 1);
M = zeros(s, 1);
i = 1;
while i <= s
try
fit = polyfit(t(inds(i)-a : inds(i)+a), f(inds(i)-a : inds(i)+a), 2);
T(i) = - fit(2) / (2*fit(1));
M(i) = fit(3) - fit(2)^2 / (4*fit(1));
i = i + 1;
catch
T(i) = [];
M(i) = [];
inds(i) = [];
s = s - 1;
end
end
end
c4science · Help