Octave Encyclopedia

demo-06

실험 데이터(mat) 읽기, FFT, Sxx, SPL

course/basic/demo-06.m

코드 인덱스로 돌아가기

카테고리

Course Basic

학습 소스 그룹

코드 길이

174

lines

작성자

won sunggyu

2025-04-22

패키지

none

pkg load 기준

전체 코드

전체 코드를 복사해서 Octave에서 바로 실행할 수 있습니다.

# 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]"]);

코드 해설

목적

실험 데이터(mat) 읽기, FFT, Sxx, SPL

입력

  • 스크립트 상단에서 정의한 파라미터/입력 데이터를 사용합니다.

출력

  • 그래프/figure 출력
  • 콘솔 텍스트 출력

실행 흐름

  1. 초기화
  2. 데이터 준비
  3. 데이터 연산
  4. 그래프 그리기
  5. 그래프 (Time domain)
  6. [1,1] = Signal_0
  7. [2,1] = Signal_1
  8. [3,1] = Signal_2
  9. [4,1] = Signal_3
  10. [5,1] = Signal_4

핵심 함수/주제

plotyyax1lengthsignalfiguredsubplotsylabel

실습 과제

  • 샘플링 주파수나 입력 주파수를 바꿔 스펙트럼 변화를 비교해보세요.
  • 질량/감쇠/강성 또는 전달함수 계수를 바꿔 응답 변화를 확인해보세요.
  • 축 범위와 라벨을 바꿔 그래프 해석성이 어떻게 달라지는지 확인해보세요.

학습 팁

  • FFT 결과는 샘플링 주파수(fs)와 길이(nn) 설정에 민감하므로 먼저 축 정의를 확인하세요.
  • 그래프 비교 시 축 범위(XLim/YLim)와 단위를 먼저 고정하면 해석 오류를 줄일 수 있습니다.
  • 입력 파일 경로가 현재 작업 디렉터리 기준인지 먼저 확인하세요.

같은 카테고리의 다른 코드