demo-08b

# filename: demo-08b.m
# writer: won sunggyu
# date: 2025-04-22
# language: octave
# description: 디지탈 Chebyshev Type I 필터(low, high, bandpass, stop)와 파워 스펙트럼

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

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

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

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

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

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] = 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_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] = 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)

ny = 5;  # number of variations
yy = zeros(ny, nx);  # signal values
YY = zeros(ny, nx);  # FFT values
Yd = zeros(ny, nx);  # dB values

fc = [100, 300];  # cut-off frequency [Hz]
border = 4;  # filter order
rp = 1;  # 통과대역 리플 (dB)
rs = 40;  % stopband에서 최소 감쇠량 (dB)
[b_loww, a_loww] = cheby1(border, rp, fc(1) / (fs/2), "low");  # Chebyshev filter Type I
[b_high, a_high] = cheby1(border, rp, fc(2) / (fs/2), "high");  # Chebyshev filter Type I
[b_band, a_band] = cheby1(border, rp, [fc(1), fc(2)] / (fs/2));  # Chebyshev filter Type I
[b_stop, a_stop] = cheby1(border, rp, [fc(1), fc(2)] / (fs/2), "stop");  # Chebyshev filter Type I

yy(1,:) = s1.y_values.values;  # signal values
yy(2,:) = filtfilt(b_loww, a_loww, yy(1,:));  # filtered signal
yy(3,:) = filtfilt(b_high, a_high, yy(1,:));  # filtered signal
yy(4,:) = filtfilt(b_band, a_band, yy(1,:));  # filtered signal
yy(5,:) = filtfilt(b_stop, a_stop, yy(1,:));  # filtered signal

for i = 1:ny
    YY(i,:) = fft(yy(i,:)) / nx;  # FFT
    Yd(i,:) = 20 * log10(abs(YY(i,:)) / log_references{1});  # dB scale
end

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

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

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);

for i = 1:ny
    plot(ax1(i, 1), tt, yy(i, :));
    plot(ax1(i, 2), ffw, Ydw(i, :));

    ylabel(ax1(i, 1), ["Amplitude [" labels{i} "]"]);
    ylabel(ax1(i, 2), ["Magnitude [" labels{i} "]"]);

    set(ax1(i, 1), "Xlabel", "Time [s]", "Ylabel", "Amplitude", "Xlim", [0, Tr]);
    set(ax1(i, 2), "Xlabel", "Frequency [Hz]", "Ylabel", "Magnitude [dB]");
end

for i = 1:ny
  vertical(ax1(i, 2), fc(1), "Color", color_base(2, :), "Linestyle", "--", "Linewidth", 1.5);
  vertical(ax1(i, 2), fc(2), "Color", color_base(2, :), "Linestyle", "--", "Linewidth", 1.5);
end

text_list = {
  "Original Signal"
  "Lowpass Filtered Signal"
  "Highpass Filtered Signal"
  "Bandpass Filtered Signal"
  "Bandstop Filtered Signal"
};

for i = 1:ny
  text(ax1(i, 2), fc(1), get(ax1(i, 2), "YLim")(2), text_list{i}, ...
      "Color", "k", "FontSize", 12, "FontWeight", "bold", ...
      "VerticalAlignment", "top", "HorizontalAlignment", "left");
end