function freq = frequency(x, y, threshold, totaltime)
% function freq = frequency(x, y, threshold,totaltime)
% Counts period of oscillation but does not discriminate between different
% frequency bands, unlike a fast fourier transform (FFT). However, when
% simulation times are short this algorithm is more acccurate than FFT.
% When totaltime is long it performs no better than FFT. It works by
% counting the number of peaks, including partial oscillation periods at
% the beginning and end of the time series
% INPUTS
% x - time vector
% y - amplitude vector
% threshold - minimum level at which oscillations are detected
% totaltime - maximum simulation time
% OUTPUTS
% freq - frequency
[up, down] = peakdet(y(1,:), threshold, x);
if length(up) > 2 || length(down) > 2
if length(up) == length(down)
% count the number of peaks
integer_freq = length(up);
period = mean(diff(up(:,1)));
% add any partial oscillations at the front or end of oscillation.
front_adj = up(1,1)/period;
end_adj = (totaltime - up(end,1))/period;
% Sum and divide by totaltime for frequency.
freq = (integer_freq + front_adj + end_adj-1)/totaltime;
end
if length(up) < length(down)
integer_freq = length(up);
period = mean(diff(up(:,1)));
front_adj = up(1,1)/period;
end_adj = (totaltime - up(end,1))/period;
freq = (integer_freq + front_adj + end_adj-1)/totaltime;
end
if length(up) > length(down)
integer_freq = length(down);
period = mean(diff(down(:,1)));
front_adj = down(1,1)/period;
end_adj = (totaltime - down(end,1))/period;
freq = (integer_freq + front_adj + end_adj-1)/totaltime;
end
else
freq = 0;
end
if isempty(up) || isempty(down)
freq = 0;
end
end |