demo-09

# filename: demo-09.m
# writer: won sunggyu
# date: 2025-04-22
# language: octave
# description: 데이터의 봉우리 찾기 (findpeaks)

#------------------------------------------------------------------------------
# 초기화
#------------------------------------------------------------------------------

run("startup.m");
printf(fmt("{mfilename}\n", "#FF5733"));

#------------------------------------------------------------------------------
# 데이터 준비
#------------------------------------------------------------------------------

filename = "comp_2k_55Hz.mat";
data = load(filename);

# fieldnames(data)  # 구조체 data의 필드 이름을 확인
# {
#   [1,1] = Signal_0
#   [2,1] = Signal_1
# }

signal = [
  data.Signal_0
  data.Signal_1
];
clear("data");

labels = cell(1, length(signal));  # signal labels
for i = 1:length(signal)
  s = signal(i);
  labels{i} = s.y_values.quantity.label;
end

# {
#     [1,1] = A
#     [1,2] = m/s^2
# }

log_references = cell(1, length(signal));  # log reference 
for i = 1:length(signal)
  s = signal(i);
  log_references{i} = s.y_values.quantity.unit_transformation.log_reference;
end

# {
#   [1,1] = 1
#   [1,2] = 1.0000e-06
# }

#------------------------------------------------------------------------------
# 데이터 연산
#------------------------------------------------------------------------------

s1 = signal(2);
dt = s1.x_values.increment;  # sampling time
nx = s1.x_values.number_of_values;  # number of samples
tt = 0:dt:dt*(nx-1);  # time vector

[idx, tt, yy, ff, Tr, dt, fs, df, nx] = transient_cut(tt, s1.y_values.values(:, 1)', 1.0, 3.0);

YY = fft(yy) / nx;  # FFT
Yd = 20 * log10(abs(YY) / log_references{2});  # dB scale

# Peak 찾기 (대상: dB scale Yd)
mpheight = 66.5;  # dB
mpdist = 9;  # number of points
[pks, locs] = findpeaks(Yd, "MinPeakHeight", mpheight, "MinPeakDistance", mpdist);

#------------------------------------------------------------------------------
# 그래프 그리기
#------------------------------------------------------------------------------

# 그래프
figured("Size", [1440, 960], "Move", [-1280, 0], "Name", mfilename);
ax1 = subplots(2, 1);

xR = 600;  # xR < fs/2 [Hz]
nw = find(ff > xR, 1, "first");
nw = nw - 1;  # number of points in the range [0, xR]
ffw = ff(1:nw);
Ydw = Yd(1:nw);

plot(ax1(1), tt, yy);
plot(ax1(2), ffw, Ydw);
scatter(ax1(2), ff(locs), pks, "r", "filled", "MarkerEdgeColor", "k", "MarkerFaceColor", "r");
horizontal(ax1(2), mpheight, "Color", "k", "LineWidth", 0.5);

set(ax1(1), "Xlabel", "Time [s]", "Ylabel", ["Amplitude [" labels{2} "]"], "Xlim", [0, Tr]);
set(ax1(2), "Xlabel", "Frequency [Hz]", "Ylabel", ["Magnitude [" labels{2} "]"], "Xlim", [0, xR]);