Power spectra of bottom pressure at D5, Digital filtering

I want to draw a graph like the picture

But I don't know what the blue line means in subplot(3,2,5) subplot(3,2,6). Please let me know the meaning of the blue line

%%
% 데이터 불러오기
%{
data_prs = load('D5_SN111 prs.dat'); % 해저 압력 데이터 파일
data_tau = load('D5_SN111 tau.dat'); % 염분성분 데이터 파일
data_current=load('D5_SN111.dat'); % u, v성분 데이터 파일
%}
data_prs = load('D5_SN111 prs.txt'); % 해저 압력 데이터 파일
data_tau = load('D5_SN111 tau.txt'); % 염분성분 데이터 파일
data_current=load('D5_SN111.txt'); % u, v성분 데이터 파일

% 해저 압력 데이터 준비
bottom_pressure = data_prs(:, end); % 해저 압력 데이터를 마지막 열로 가정
bottom_pressure=fillmissing(bottom_pressure,'linear'); % NaN 제거

% u, v 성분 데이터 준비
u_component = data_current(:, 7); 
v_component = data_current(:, 8); 
u_component=fillmissing(u_component,'linear'); % NaN 제거
v_component=fillmissing(v_component,'linear'); % NaN 제거

% 파라미터 설정
WINDOW = 1024 * 2;
NOVERLAP = 512 * 2;
NFFT = 1024 * 2;
Fs = 1; % 샘플링 주파수 (1시간 간격 가정)
confidence_level = 0.95; % 95% 신뢰 구간
confcen = 10; % 
conf_level = 100; 

%pressure at D5
[Pxx,F,Pxxc]=pwelch(bottom_pressure,WINDOW,NOVERLAP,NFFT,Fs,'ConfidenceLevel',confidence_level);

confcen=10;
confup=Pxxc(10,2)./Pxx(10)*confcen;
confdn=Pxxc(10,1)./Pxx(10)*confcen;

figure;
loglog(F, Pxx, 'k', 'LineWidth', 1); % PSD
hold on;
loglog(F(conf_level*[1 1]), Pxxc(conf_level,:), 'r', 'LineWidth', 1.5)
xlabel('Frequency (cph)');
ylabel('Power/Frequency ((dbar)^2/cph)');
title('Power Spectrum of Bottom Pressure at D5 with 95% Confidence Interval');
legend('PSD', 'Upper 95% Confidence', 'Lower 95% Confidence');

Warning: Ignoring extra legend entries.

grid on;
hold off;

% u 성분 PSD 계산 및 그래프
[Pxx, F, Pxxc] = pwelch(u_component, WINDOW, NOVERLAP, NFFT, Fs, 'ConfidenceLevel', confidence_level);

figure;
loglog(F, Pxx, 'k', 'LineWidth', 1); % PSD
hold on;
loglog(F(conf_level*[1 1]), Pxxc(conf_level,:), 'r', 'LineWidth', 1.5)
xlabel('Frequency (cph)');
ylabel('Power/Frequency ((cm/s)^2/cph)');
title('Power Spectrum of U Component at D5 with 95% Confidence Interval');
legend('PSD', 'Upper 95% Confidence', 'Lower 95% Confidence');

Warning: Ignoring extra legend entries.

grid on;
hold off;

% v 성분 PSD 계산 및 그래프
[Pxx, F, Pxxc] = pwelch(v_component, WINDOW, NOVERLAP, NFFT, Fs, 'ConfidenceLevel', confidence_level);

figure;
loglog(F, Pxx, 'k', 'LineWidth', 1); % PSD
hold on;
loglog(F(conf_level*[1 1]), Pxxc(conf_level,:), 'r', 'LineWidth', 1.5)
xlabel('Frequency (cph)');
ylabel('Power/Frequency ((cm/s)^2/cph)');
title('Power Spectrum of V Component at D5 with 95% Confidence Interval');
legend('PSD', 'Upper 95% Confidence', 'Lower 95% Confidence');

Warning: Ignoring extra legend entries.

grid on;
hold off;

%%
% 샘플 데이터 생성 (또는 실제 데이터로 대체)
dt = 1; % 샘플링 간격 (시간 단위: 1시간 간격)
x = bottom_pressure; % 예제 데이터 (해당하는 데이터로 대체)

% 1) 대역 통과 필터 설정 (컷오프 주기 11-15시간, 6차 Butterworth 필터)
tlow = 15; % 하한 주기 (시간 단위)
thigh = 11; % 상한 주기 (시간 단위)
n_bandpass = 6; % 필터 차수

flow = 2 * dt / thigh; % 하한 주파수
fhigh = 2 * dt / tlow; % 상한 주파수
[b_bandpass, a_bandpass] = butter(n_bandpass, flow); % 대역 통과 필터 설계
y_bandpass = filtfilt(b_bandpass, a_bandpass, x); % 대역 통과 필터 적용

% 2) 저역 필터 설정 (컷오프 주기 48시간, 3차 Butterworth 필터)
tlow_lowpass = 48; % 저역 필터 컷오프 주기 (시간 단위)
n_lowpass = 3; % 저역 필터 차수

flow_lowpass = 2 * dt / tlow_lowpass; % 저역 필터 주파수
[b_lowpass, a_lowpass] = butter(n_lowpass, flow_lowpass, 'low'); % 저역 필터 설계
y_lowpass = filtfilt(b_lowpass, a_lowpass, x); % 저역 필터 적용

% 결과 그래프
figure;
subplot(3, 1, 1);
plot(x);
title('Original Data');
xlabel('Time (hours)');
ylabel('Amplitude');

subplot(3, 1, 2);
plot(y_bandpass);
title('Band-pass Filtered Data (11-15 hours)');
xlabel('Time (hours)');
ylabel('Amplitude');

subplot(3, 1, 3);
plot(y_lowpass);
title('Low-pass Filtered Data (48 hours)');
xlabel('Time (hours)');
ylabel('Amplitude');
%%
subplot(3,1,1);
plot(y_bandpass);
title('Band-pass Filtered Data (11-15 hours)');
xlabel('Time (hours)');
ylabel('Amplitude');

subplot(3,1,2);
plot(y_lowpass);
title('Low-pass Filtered Data (48 hours)');
xlabel('Time (hours)');
ylabel('Amplitude');

subplot(3,2,5);
% Power Spectral Density (PSD) 및 신뢰구간 계산
[Pxx, F, Pxxc] = pwelch(bottom_pressure, WINDOW, NOVERLAP, NFFT, Fs, 'ConfidenceLevel', confidence_level);

% 그래프 1: PSD와 신뢰구간 상한선
loglog(F, Pxx, 'k', 'LineWidth', 1); % PSD (검은색 실선)
hold on;
loglog(F(conf_level*[1 1]), Pxxc(conf_level,:), 'r', 'LineWidth', 1.5)
hold on;
loglog(F, Pxxc(:, 2), 'b--', 'LineWidth', 1); % 신뢰구간 상한선 (파란색 점선)
xlabel('Frequency (cph)');
ylabel('Power/Frequency ((dbar)^2/cph)');
title('PSD and Upper 95% Confidence Interval');
legend('PSD', 'Upper 95% Confidence');
grid on;
hold off;

subplot(3,2,6);
loglog(F, Pxx, 'k', 'LineWidth', 1); % PSD (검은색 실선)
hold on;
loglog(F(conf_level*[1 1]), Pxxc(conf_level,:), 'r', 'LineWidth', 1.5)
loglog(F, Pxxc(:, 1), 'b--', 'LineWidth', 1); % 신뢰구간 하한선 (파란색 점선)
xlabel('Frequency (cph)');
ylabel('Power/Frequency ((dbar)^2/cph)');
title('PSD and Lower 95% Confidence Interval');
legend('PSD', 'Lower 95% Confidence');
grid on;
hold off;

Answer 1

[Walter Roberson]

subplot(3,2,5);
% Power Spectral Density (PSD) 및 신뢰구간 계산
[Pxx, F, Pxxc] = pwelch(bottom_pressure, WINDOW, NOVERLAP, NFFT, Fs, 'ConfidenceLevel', confidence_level);

So Pxxc is going to be,

Confidence bounds, returned as a matrix with real-valued elements. The row size of the matrix is equal to the length of the PSD estimate, pxx. pxxc has twice as many columns as pxx. Odd-numbered columns contain the lower bounds of the confidence intervals, and even-numbered columns contain the upper bounds. Thus, pxxc(m,2*n-1) is the lower confidence bound and pxxc(m,2*n) is the upper confidence bound corresponding to the estimate pxx(m,n). The coverage probability of the confidence intervals is determined by the value of the probability input.

loglog(F, Pxxc(:, 2), 'b--', 'LineWidth', 1); % 신뢰구간 상한선 (파란색 점선)

So you are plotting the upper bounds of the confidence interval as the blue line.

[은진]

Thank you so much!!!