main-03

sinusoid signal and auto-power spectrum

ex-03/main-03.m

코드 인덱스로 돌아가기

카테고리

Exercise 03

학습 소스 그룹

코드 길이

138

lines

작성자

won sunggyu

2025-03-26

패키지

none

pkg load 기준

전체 코드

# filename: main-03.m
# writer: won sunggyu
# date: 2025-03-26
# language: octave
# description: sinusoid signal and auto-power spectrum

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

# Save sinusoid parameters to a text file
# save_struct_txt("s1.txt", ss); # s1, s2

# Load sinusoid parameters from a text file as a struct
ss = load_struct_txt("s1.txt"); # s1, s2

# Sinusoid
[~, x1] = generate_sinusoid(ss.frequency, ss.phase, ss.amplitude, ss.fs, ss.duration); % 1xN, 1xN

# Sampling
[fs, duration, dt, df, nn] = sampling_settings(ss.fs, ss.duration);

# Create Domain
[tt, ff] = make_axes(ss.fs, ss.duration); # 1xN, 1xN

# hanning window coefficient
# - magnitude equivalent
# - power equivalent
[hann_m, hann_p] = deal(2.0, 1.633);

# Hann Window
hann_c = hann_p;
hann_w = hanning(nn).' * hann_c; % 1xN [0, hann_c]

# Windowed Signal
xw = x1 .* hann_w; % 1xN

# FFT and Auto-Power Spectrum
X1 = fft(x1) / nn;
Xw = fft(xw) / nn;
Sxx1 = real(conj(X1) .* X1); % 1xN
Sxxw = real(conj(Xw) .* Xw); % 1xN

# FFT Shift
ff = ff - df * nn / 2;
Sxx1 = fftshift(Sxx1); % 1xN
Sxxw = fftshift(Sxxw); % 1xN

# Power sum
Psum1 = power_sum(x1); # integration of continuous periodic time signal 
Psumw = power_sum(xw); # integration of continuous periodic time signal 
PSxx1 = sum(Sxx1); # sum of discrete series
PSxxw = sum(Sxxw); # sum of discrete series

pow_ratio_t = Psumw / Psum1;
pow_ratio_f = PSxxw / PSxx1;

# printf(fmt("power sum\n", "#69A1FA"))
# printf(fmt("T axis: Psumw = {Psumw}, Psum1 = {Psum1}, ratio = {pow_ratio_t}\n"))
# printf(fmt("F axis: PSxxw = {PSxxw}, PSxx1 = {PSxx1}, ratio = {pow_ratio_f}\n"))

printf(fmt("power sum\n", "#69A1FA"))
printf(fmt("T axis----------, F axis----------,\n"))
printf(fmt("Psumw = {Psumw:8.4f}, PSxxw = {PSxxw:8.4f},\n"))
printf(fmt("Psum1 = {Psum1:8.4f}, PSxx1 = {PSxx1:8.4f},\n"))
printf(fmt("ratio = {pow_ratio_t:8.4f}, ratio = {pow_ratio_f:8.4f},\n"))

# Decibell
aref = 1e-6;
X1h = abs(X1)(2: nn / 2) * 2 / sqrt(2); # one-sided FFT abs
Xwh = abs(Xw)(2: nn / 2) * 2 / sqrt(2); # one-sided FFT abs
decibell1 = decibell_overall(X1h, aref);
decibellw = decibell_overall(Xwh, aref);
db_ratio = decibellw - decibell1;

printf(fmt("decibel overall\n", "#69A1FA"))
printf(fmt("decibellw = {decibellw} dB\n"))
printf(fmt("decibell1 = {decibell1} dB\n"))
printf(fmt("ratio     = {db_ratio} dB\n"))

# Visualization
param_f = {"Size", [960, 960], "Name", "Sinusoid"};
param_at = {"Xlabel", "Time [sec]", "Ylabel", "Magnitude [m/s^2]", "Xlim", [0, ss.duration]};
param_af = {"Xlabel", "Frequency [Hz]", "Ylabel", "Auto-Power [(m/s^2)^2]", "Xlim", ss.xrng, "Ylim", ss.yrng};

figured(param_f);
ax = subplots(2, 1);

plotd(ax(1), tt, x1, ";x1;");
plotd(ax(1), tt, xw, ";xw;");
plotd(ax(2), ff, Sxx1, ";Sxx1;");
plotd(ax(2), ff, Sxxw, ";Sxxw;");

textd(ax(1), 0.05, 0.90, "transient", "Color", "#69A1FA");
textd(ax(1), 0.05, 0.80, fmt("transient-hanning [ c={hann_c:.3f} ]"), "Color", "#CF87DA");
textd(ax(2), 0.05, 0.90, "two-sided", "Color", "#69A1FA");
textd(ax(2), 0.05, 0.80, fmt("two-sided-hanning [ c={hann_c:.3f} ]"), "Color", "#CF87DA");

set(ax(1), param_at{:});
set(ax(2), param_af{:});
set(ax(2), "Ylim", [0, 25]);

# 
# 소수가 원하는 답변을 못하면 개별 지도하라.
# 다수가 원하는 답변을 못하면 문제를 다시 작성하라.
# 신호를 생성하고 퍼시벌의 정리가 성립함을 보여라.
# 신호를 생성하고 윈도우(hanning)을 적용하고 퍼시벌의 정리가 성립함을 보여라.
# 
# 상황
# - 값을 구하고 비교하라는 주문을 했다.
# - 피주문자들은 주문에 따라 값을 구하였지만 그 값들을 비교하라는 정확한 표현이 없어서 비교를 하지 않는다.
# - 해당 값들은 퍼시벌 정리에 해당하며 항상 머리 속에 있는 상식 수준이므로 당연히 비교라는 단어에 반응했어야 한다.
# - 피주문자들은 어느 값들을 비교하라는 주문을 명시하지 않은 이유로 비교하지 않았다.
# 
# 문제
# - create sinusoidal signal (acceleration, m/s^2)
# - power sum
# - calculate decibell
# - FFT on the signal 
# - Auto-Power Spectrum
# - power sum
# - calculate decibell
# - apply hanning window to the signal
# - power sum
# - calculate decibell
# - FFT on the windowed signal
# - Auto-Power Spectrum
# - power sum
# - calculate decibell
# 
# 주문
# - signal과 windowed signal을 비교하라.
# - signal의 auto-power spectrum와 windowed signal의 auto-power spectrum을 비교하라.
# - power sum 네 값을 비교하라. 네 개의 power sum은 같아야 한다.
# - decibell를 비교하라. 네 개의 decibell은 같아야 한다. 가속도 레퍼런스는 1e-6이다.
#

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

# filename: main-03.m
# writer: won sunggyu
# date: 2025-03-26
# language: octave
# description: sinusoid signal and auto-power spectrum

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

# Save sinusoid parameters to a text file
# save_struct_txt("s1.txt", ss); # s1, s2

# Load sinusoid parameters from a text file as a struct
ss = load_struct_txt("s1.txt"); # s1, s2

# Sinusoid
[~, x1] = generate_sinusoid(ss.frequency, ss.phase, ss.amplitude, ss.fs, ss.duration); % 1xN, 1xN

# Sampling
[fs, duration, dt, df, nn] = sampling_settings(ss.fs, ss.duration);

# Create Domain
[tt, ff] = make_axes(ss.fs, ss.duration); # 1xN, 1xN

# hanning window coefficient
# - magnitude equivalent
# - power equivalent
[hann_m, hann_p] = deal(2.0, 1.633);

# Hann Window
hann_c = hann_p;
hann_w = hanning(nn).' * hann_c; % 1xN [0, hann_c]

# Windowed Signal
xw = x1 .* hann_w; % 1xN

# FFT and Auto-Power Spectrum
X1 = fft(x1) / nn;
Xw = fft(xw) / nn;
Sxx1 = real(conj(X1) .* X1); % 1xN
Sxxw = real(conj(Xw) .* Xw); % 1xN

# FFT Shift
ff = ff - df * nn / 2;
Sxx1 = fftshift(Sxx1); % 1xN
Sxxw = fftshift(Sxxw); % 1xN

# Power sum
Psum1 = power_sum(x1); # integration of continuous periodic time signal 
Psumw = power_sum(xw); # integration of continuous periodic time signal 
PSxx1 = sum(Sxx1); # sum of discrete series
PSxxw = sum(Sxxw); # sum of discrete series

pow_ratio_t = Psumw / Psum1;
pow_ratio_f = PSxxw / PSxx1;

# printf(fmt("power sum\n", "#69A1FA"))
# printf(fmt("T axis: Psumw = {Psumw}, Psum1 = {Psum1}, ratio = {pow_ratio_t}\n"))
# printf(fmt("F axis: PSxxw = {PSxxw}, PSxx1 = {PSxx1}, ratio = {pow_ratio_f}\n"))

printf(fmt("power sum\n", "#69A1FA"))
printf(fmt("T axis----------, F axis----------,\n"))
printf(fmt("Psumw = {Psumw:8.4f}, PSxxw = {PSxxw:8.4f},\n"))
printf(fmt("Psum1 = {Psum1:8.4f}, PSxx1 = {PSxx1:8.4f},\n"))
printf(fmt("ratio = {pow_ratio_t:8.4f}, ratio = {pow_ratio_f:8.4f},\n"))

# Decibell
aref = 1e-6;
X1h = abs(X1)(2: nn / 2) * 2 / sqrt(2); # one-sided FFT abs
Xwh = abs(Xw)(2: nn / 2) * 2 / sqrt(2); # one-sided FFT abs
decibell1 = decibell_overall(X1h, aref);
decibellw = decibell_overall(Xwh, aref);
db_ratio = decibellw - decibell1;

printf(fmt("decibel overall\n", "#69A1FA"))
printf(fmt("decibellw = {decibellw} dB\n"))
printf(fmt("decibell1 = {decibell1} dB\n"))
printf(fmt("ratio     = {db_ratio} dB\n"))

# Visualization
param_f = {"Size", [960, 960], "Name", "Sinusoid"};
param_at = {"Xlabel", "Time [sec]", "Ylabel", "Magnitude [m/s^2]", "Xlim", [0, ss.duration]};
param_af = {"Xlabel", "Frequency [Hz]", "Ylabel", "Auto-Power [(m/s^2)^2]", "Xlim", ss.xrng, "Ylim", ss.yrng};

figured(param_f);
ax = subplots(2, 1);

plotd(ax(1), tt, x1, ";x1;");
plotd(ax(1), tt, xw, ";xw;");
plotd(ax(2), ff, Sxx1, ";Sxx1;");
plotd(ax(2), ff, Sxxw, ";Sxxw;");

textd(ax(1), 0.05, 0.90, "transient", "Color", "#69A1FA");
textd(ax(1), 0.05, 0.80, fmt("transient-hanning [ c={hann_c:.3f} ]"), "Color", "#CF87DA");
textd(ax(2), 0.05, 0.90, "two-sided", "Color", "#69A1FA");
textd(ax(2), 0.05, 0.80, fmt("two-sided-hanning [ c={hann_c:.3f} ]"), "Color", "#CF87DA");

set(ax(1), param_at{:});
set(ax(2), param_af{:});
set(ax(2), "Ylim", [0, 25]);

# 
# 소수가 원하는 답변을 못하면 개별 지도하라.
# 다수가 원하는 답변을 못하면 문제를 다시 작성하라.
# 신호를 생성하고 퍼시벌의 정리가 성립함을 보여라.
# 신호를 생성하고 윈도우(hanning)을 적용하고 퍼시벌의 정리가 성립함을 보여라.
# 
# 상황
# - 값을 구하고 비교하라는 주문을 했다.
# - 피주문자들은 주문에 따라 값을 구하였지만 그 값들을 비교하라는 정확한 표현이 없어서 비교를 하지 않는다.
# - 해당 값들은 퍼시벌 정리에 해당하며 항상 머리 속에 있는 상식 수준이므로 당연히 비교라는 단어에 반응했어야 한다.
# - 피주문자들은 어느 값들을 비교하라는 주문을 명시하지 않은 이유로 비교하지 않았다.
# 
# 문제
# - create sinusoidal signal (acceleration, m/s^2)
# - power sum
# - calculate decibell
# - FFT on the signal 
# - Auto-Power Spectrum
# - power sum
# - calculate decibell
# - apply hanning window to the signal
# - power sum
# - calculate decibell
# - FFT on the windowed signal
# - Auto-Power Spectrum
# - power sum
# - calculate decibell
# 
# 주문
# - signal과 windowed signal을 비교하라.
# - signal의 auto-power spectrum와 windowed signal의 auto-power spectrum을 비교하라.
# - power sum 네 값을 비교하라. 네 개의 power sum은 같아야 한다.
# - decibell를 비교하라. 네 개의 decibell은 같아야 한다. 가속도 레퍼런스는 1e-6이다.
#

코드 해설

목적

sinusoid signal and auto-power spectrum

입력

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

출력

콘솔 텍스트 출력

실행 흐름

Sampling
hanning window coefficient
Hann Window
Windowed Signal
FFT and Auto-Power Spectrum
FFT Shift
Visualization
- create sinusoidal signal (acceleration, m/s^2)
- FFT on the signal
- apply hanning window to the signal

핵심 함수/주제

fmtprintfaxplotdtextdsetabsconj

실습 과제

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

학습 팁

FFT 결과는 샘플링 주파수(fs)와 길이(nn) 설정에 민감하므로 먼저 축 정의를 확인하세요.
입력 파일 경로가 현재 작업 디렉터리 기준인지 먼저 확인하세요.

같은 카테고리의 다른 코드

main_splecgram
ex-03/main_splecgram.m
main-03
ex-03/main_spwm_plotyy.m
main_spwm
ex-03/main_spwm.m