Power spectra of bottom pressure at D5, Digital filtering

57 次查看（过去 30 天）
은진
은진 2024-11-16，5:11
评论： 은진 2024-11-17，13:14
Ran in:
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;

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2024-11-16，6:11
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.

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Spectral Estimation 的更多信息

标签

产品


版本

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by