demo-06

# filename: demo-06.m
# writer: won sunggyu
# date: 2025-04-22
# language: octave
# description: 실험 데이터(mat) 읽기, FFT, Sxx, SPL

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

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

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

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

# fieldnames(data)  # 구조체 data의 필드 이름을 확인
# {
#   [1,1] = Signal_0
#   [2,1] = Signal_1
#   [3,1] = Signal_2
#   [4,1] = Signal_3
#   [5,1] = Signal_4
#   [6,1] = Signal_5
#   [7,1] = Signal_6
#   [8,1] = Signal_7
# }

# 데이터를 필요한 형태로 저장
signal = [
  data.Signal_0
  data.Signal_1
  data.Signal_2
  data.Signal_3
  data.Signal_4
  data.Signal_5
  data.Signal_6
  data.Signal_7
];

# 불필요한 변수 삭제
clear("data");

# structuare of single signal
# s1 = signal(1);
# s1.function_record;  # metadata
# s1.x_values.increment;
# s1.x_values.number_of_values;
# s1.x_values.label;
# s1.y_values.values;
# s1.y_values.quantity.label;
# s1.y_values.quantity.unit_transformation.log_reference;

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

# disp(labels);  # display signal labels
# {
#     [1,1] = Pa
#     [1,2] = Pa
#     [1,3] = Pa
#     [1,4] = Pa
#     [1,5] = Pa
#     [1,6] = Pa
#     [1,7] = Pa
#     [1,8] = m/s^2
# }

# log reference 정보 수집
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

# disp(log_references)  # display log reference
# {
#     [1,1] = 2.0000e-05
#     [1,2] = 2.0000e-05
#     [1,3] = 2.0000e-05
#     [1,4] = 2.0000e-05
#     [1,5] = 2.0000e-05
#     [1,6] = 2.0000e-05
#     [1,7] = 2.0000e-05
#     [1,8] = 1.0000e-06
# }

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

s1 = signal(1);
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
Tr = nx * dt;  # total time
df = 1 / Tr;  # frequency resolution
fs = 1 / dt;  # sampling frequency
ff = 0:df:df*(nx-1);  # frequency vector
f2 = ff - df * nx / 2; # frequency vector (shifted)

nchan = length(signal);  # number of channels
yy = zeros(nchan, nx);  # signal values
YY = zeros(nchan, nx);  # FFT
Yd = zeros(nchan, nx);  # dB scale

for i = 1:nchan
    s = signal(i);
    yy(i, :) = s.y_values.values;
    YY(i, :) = fft(yy(i, :)) / nx;
    Yd(i, :) = 20 * log10(abs(YY(i, :)) / log_references{i});
end

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

# 그래프 (Time domain)
figured("Size", [1440, 960], "Move", [0, 0], "Name", mfilename);
ax1 = subplots(4, 2, "Xlabel", "Time (s)", "Xlim", [0, Tr]);

plot(ax1(1, 1), tt, yy(1, :));
plot(ax1(2, 1), tt, yy(2, :));
plot(ax1(3, 1), tt, yy(3, :));
plot(ax1(4, 1), tt, yy(4, :));

plot(ax1(1, 2), tt, yy(5, :));
plot(ax1(2, 2), tt, yy(6, :));
plot(ax1(3, 2), tt, yy(7, :));
plot(ax1(4, 2), tt, yy(8, :));

for i = 1:nchan
    ylabel(ax1(i), ["Magnitude [" labels{i} "]"]);
end

# 그래프 (FFT)
figured("Size", [1440, 960], "Move", [0, 0], "Name", mfilename);
ax2 = subplots(4, 2, "Xlabel", "Frequency (Hz)", "Xlim", [-200, 200]);

for j = 1:2
    for i = 1:4
        idx = (j-1) * 4 + i;
        plot(ax2(i, j), f2, fftshift(YY(idx, :)));
    end
end

for i = 1:nchan
  ylabel(ax2(i), ["Magnitude [" labels{i} "]"]);
end

# 그래프 (Sxx decibell)
figured("Size", [1440, 960], "Move", [0, 0], "Name", mfilename);
ax4 = subplots(1, 1, "Xlabel", "Frequency [Hz]");

xR = 300;  # 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);

for i = 1:nchan
  plot(ax4, ffw, Ydw(1, :));
end

ylabel(ax4, ["Magnitude [dB]"]);