Octave Encyclopedia

demo-08a

디지탈 버터워스 필터(low, high, bandpass, stop)와 파워 스펙트럼

course/basic/demo-08a.m

코드 인덱스로 돌아가기

카테고리

Course Basic

학습 소스 그룹

코드 길이

150

lines

작성자

won sunggyu

2025-04-22

패키지

none

pkg load 기준

전체 코드

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

# filename: demo-08a.m
# writer: won sunggyu
# date: 2025-04-22
# language: octave
# description: 디지탈 버터워스 필터(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 of digital filter
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] = butter(border, fc(1) / (fs/2), "low");  # Butterworth filter
[b_high, a_high] = butter(border, fc(2) / (fs/2), "high");  # Butterworth filter
[b_band, a_band] = butter(border, [fc(1), fc(2)] / (fs/2), "bandpass");  # Butterworth filter
[b_stop, a_stop] = butter(border, [fc(1), fc(2)] / (fs/2), "stop");  # Butterworth filter

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

코드 해설

목적

디지탈 버터워스 필터(low, high, bandpass, stop)와 파워 스펙트럼

입력

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

출력

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

실행 흐름

  1. 초기화
  2. 데이터 준비
  3. 데이터 연산
  4. 그래프 그리기
  5. 그래프

핵심 함수/주제

yyax1fcbutterfiltfiltlengthsignalzeros

실습 과제

  • 샘플링 주파수나 입력 주파수를 바꿔 스펙트럼 변화를 비교해보세요.
  • 같은 연산을 내장 함수와 사용자 함수 두 방식으로 계산해 오차를 비교해보세요.
  • 질량/감쇠/강성 또는 전달함수 계수를 바꿔 응답 변화를 확인해보세요.

학습 팁

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

같은 카테고리의 다른 코드