How can I solve this confidence level error?

I was trying to plot data of power spectra of pressure, u, v. However, the confidence level keep plotting strange.

I think my confidence level is missing something, or reason that I didn't apply code 'butter'.

Attached file is on .txt, original file is on .dat

clear all
close all
clc

data_uv = load('D5_SN111.dat'); % u, v data
data_pressure = load('D5_SN111(1).dat'); % pressure data

u = data_uv(:, 7); % u 성분 (7번째 열)
v = data_uv(:, 8); % v 성분 (8번째 열)
% D5_SN111(1).dat: year month day hour minute sec pressure
pressure = data_pressure(:, 7); % pressure data

pressure = fillmissing(pressure, 'linear'); 
u = fillmissing(u, 'linear'); 
v = fillmissing(v, 'linear'); 

WINDOW = 1024 * 2;
NOVERLAP = 512 * 2;
NFFT = 1024 * 2;
Fs = 1; % (1 sample/hour)

confcen = 10; % 
conf_level = 100; % 
%% pressure data PSD
figure;
[pxx_pressure, f_pressure, Pxxc_pressure] = pwelch(pressure, WINDOW, NOVERLAP, NFFT, Fs, 'ConfidenceLevel', 0.95);
loglog(f_pressure, pxx_pressure, 'k'); % 파워 스펙트럼 밀도 그래프
hold on;

confup = Pxxc_pressure(conf_level, 2) ./ pxx_pressure(conf_level) * confcen;
confdn = Pxxc_pressure(conf_level, 1) ./ pxx_pressure(conf_level) * confcen;
loglog([f_pressure(conf_level), f_pressure(conf_level)], [confdn, confup], 'r', 'LineWidth', 1.5);
title('Power spectra of bottom pressure at D5');
xlabel('Frequency (Hz)');
ylabel('Power/frequency (dbar^2/Hz)');
grid on;
legend('PSD', '95% Confidence Interval');
hold off;
%% u data's PSD
figure;
[pxx_u, f_u, Pxxc_u] = pwelch(u, WINDOW, NOVERLAP, NFFT, Fs, 'ConfidenceLevel', 0.95);
loglog(f_u, pxx_u, 'b'); 
hold on;

confup = Pxxc_u(conf_level, 2) ./ pxx_u(conf_level) * confcen;
confdn = Pxxc_u(conf_level, 1) ./ pxx_u(conf_level) * confcen;
loglog([f_u(conf_level), f_u(conf_level)], [confdn, confup], 'r', 'LineWidth', 1.5);
title('Power spectra of U at D5');
xlabel('Frequency (Hz)');
ylabel('Power/frequency ((cm/s)^2/Hz)');
grid on;
legend('PSD', '95% Confidence Interval');
hold off;
%%v data's PSD
figure;
[pxx_v, f_v, Pxxc_v] = pwelch(v, WINDOW, NOVERLAP, NFFT, Fs, 'ConfidenceLevel', 0.95);
loglog(f_v, pxx_v, 'g'); 
hold on;

confup = Pxxc_v(conf_level, 2) ./ pxx_v(conf_level) * confcen;
confdn = Pxxc_v(conf_level, 1) ./ pxx_v(conf_level) * confcen;
loglog([f_v(conf_level), f_v(conf_level)], [confdn, confup], 'r', 'LineWidth', 1.5);
title('Power spectra of V at D5');
xlabel('Frequency (Hz)');
ylabel('Power/frequency ((cm/s)^2/Hz)');
grid on;
legend('PSD', '95% Confidence Interval');
hold off;

and, this is what I got.

However, I was trying to get such data.

I have no idea what I'm missing with.

Answer 1

[Star Strider]

You are plotting it incorrectly.

Use this:

loglog(f_pressure, Pxxc_pressure, 'r', 'LineWidth', 1.5);

instead of:

loglog([f_pressure(conf_level), f_pressure(conf_level)], [confdn, confup], 'r', 'LineWidth', 1.5);

and other incidents with similar code.

Note that ‘conf_level’ is a scalar with a single value (100), so it will only plot the confidence level at that one point. If that is the only value you want to plot, add that subscript to my code, and change the 'r' to a red marker of some type (perhaps '.r') so it will plot points instead of a line connecting them.

I changed all of them to produce this —

% data_uv = load('D5_SN111.dat'); % u, v data
data_uv = load('D5_SN111.txt'); % u, v data
% data_pressure = load('D5_SN111(1).dat'); % pressure data
data_pressure = load('D5_SN111(1).txt'); % pressure data
u = data_uv(:, 7); % u 성분 (7번째 열)
v = data_uv(:, 8); % v 성분 (8번째 열)
% D5_SN111(1).dat: year month day hour minute sec pressure
pressure = data_pressure(:, 7); % pressure data
pressure = fillmissing(pressure, 'linear'); 
u = fillmissing(u, 'linear'); 
v = fillmissing(v, 'linear'); 
WINDOW = 1024 * 2;
NOVERLAP = 512 * 2;
NFFT = 1024 * 2;
Fs = 1; % (1 sample/hour)
confcen = 10; % 
conf_level = 100; % 
%% pressure data PSD
figure;
[pxx_pressure, f_pressure, Pxxc_pressure] = pwelch(pressure, WINDOW, NOVERLAP, NFFT, Fs, 'ConfidenceLevel', 0.95);

figure
loglog(f_pressure, pxx_pressure, 'k'); % 파워 스펙트럼 밀도 그래프
hold on;
confup = Pxxc_pressure(conf_level, 2) ./ pxx_pressure(conf_level) * confcen

confup = 20.8534

confdn = Pxxc_pressure(conf_level, 1) ./ pxx_pressure(conf_level) * confcen

confdn = 5.8532

loglog(f_pressure, Pxxc_pressure, 'r', 'LineWidth', 1.5);
title('Power spectra of bottom pressure at D5');
xlabel('Frequency (Hz)');
ylabel('Power/frequency (dbar^2/Hz)');
grid on;
legend('PSD', '95% Confidence Interval');
hold off;

%% u data's PSD
figure;
[pxx_u, f_u, Pxxc_u] = pwelch(u, WINDOW, NOVERLAP, NFFT, Fs, 'ConfidenceLevel', 0.95);
loglog(f_u, pxx_u, 'b'); 
hold on;
confup = Pxxc_u(conf_level, 2) ./ pxx_u(conf_level) * confcen;
confdn = Pxxc_u(conf_level, 1) ./ pxx_u(conf_level) * confcen;
loglog(f_u, Pxxc_u, 'r', 'LineWidth', 1.5)
% loglog([f_u(conf_level), f_u(conf_level)], [confdn, confup], 'r', 'LineWidth', 1.5);
title('Power spectra of U at D5');
xlabel('Frequency (Hz)');
ylabel('Power/frequency ((cm/s)^2/Hz)');
grid on;
legend('PSD', '95% Confidence Interval');
hold off;

%%v data's PSD
figure;
[pxx_v, f_v, Pxxc_v] = pwelch(v, WINDOW, NOVERLAP, NFFT, Fs, 'ConfidenceLevel', 0.95);
loglog(f_v, pxx_v, 'g'); 
hold on;
confup = Pxxc_v(conf_level, 2) ./ pxx_v(conf_level) * confcen;
confdn = Pxxc_v(conf_level, 1) ./ pxx_v(conf_level) * confcen;
% loglog([f_v(conf_level), f_v(conf_level)], [confdn, confup], 'r', 'LineWidth', 1.5);
loglog(f_v, Pxxc_v, 'r', 'LineWidth', 1.5);
title('Power spectra of V at D5');
xlabel('Frequency (Hz)');
ylabel('Power/frequency ((cm/s)^2/Hz)');
grid on;
legend('PSD', '95% Confidence Interval');
hold off;

.

[kim]

Always, thanks for your help.

But I'm curious about last picture I attached. As you could see, the 95%confidence interval plotted as bar graph, which plotted on part of the graph. And the code plots whole confidence interval data.

In this case, am I missing something?

[Star Strider]

My pleasure!

I’m not certain that you’re missing anything. I’m also not certain what result you want, although I outlined how to plot the confidence interval only at element 100 in each vector. (My code plots the entiire confidence interval vector.)

To plot the conifdence inteerval only at element 100 in all the plots, try this —

% data_uv = load('D5_SN111.dat'); % u, v data
data_uv = load('D5_SN111.txt'); % u, v data
% data_pressure = load('D5_SN111(1).dat'); % pressure data
data_pressure = load('D5_SN111(1).txt'); % pressure data
u = data_uv(:, 7); % u 성분 (7번째 열)
v = data_uv(:, 8); % v 성분 (8번째 열)
% D5_SN111(1).dat: year month day hour minute sec pressure
pressure = data_pressure(:, 7); % pressure data
pressure = fillmissing(pressure, 'linear'); 
u = fillmissing(u, 'linear'); 
v = fillmissing(v, 'linear'); 
WINDOW = 1024 * 2;
NOVERLAP = 512 * 2;
NFFT = 1024 * 2;
Fs = 1; % (1 sample/hour)
confcen = 10; % 
conf_level = 100; % 
%% pressure data PSD
figure;
[pxx_pressure, f_pressure, Pxxc_pressure] = pwelch(pressure, WINDOW, NOVERLAP, NFFT, Fs, 'ConfidenceLevel', 0.95);

figure
loglog(f_pressure, pxx_pressure, 'k'); % 파워 스펙트럼 밀도 그래프
hold on;
confup = Pxxc_pressure(conf_level, 2) ./ pxx_pressure(conf_level) * confcen

confup = 20.8534

confdn = Pxxc_pressure(conf_level, 1) ./ pxx_pressure(conf_level) * confcen

confdn = 5.8532

loglog(f_pressure(conf_level*[1 1]), Pxxc_pressure(conf_level,:), 'r', 'LineWidth', 1.5);
title('Power spectra of bottom pressure at D5');
xlabel('Frequency (Hz)');
ylabel('Power/frequency (dbar^2/Hz)');
grid on;
legend('PSD', '95% Confidence Interval');
hold off;

%% u data's PSD
figure;
[pxx_u, f_u, Pxxc_u] = pwelch(u, WINDOW, NOVERLAP, NFFT, Fs, 'ConfidenceLevel', 0.95);
loglog(f_u, pxx_u, 'b'); 
hold on;
confup = Pxxc_u(conf_level, 2) ./ pxx_u(conf_level) * confcen;
confdn = Pxxc_u(conf_level, 1) ./ pxx_u(conf_level) * confcen;
loglog(f_u(conf_level*[1 1]), Pxxc_u(conf_level,:), 'r', 'LineWidth', 1.5)
% loglog([f_u(conf_level), f_u(conf_level)], [confdn, confup], 'r', 'LineWidth', 1.5);
title('Power spectra of U at D5');
xlabel('Frequency (Hz)');
ylabel('Power/frequency ((cm/s)^2/Hz)');
grid on;
legend('PSD', '95% Confidence Interval');
hold off;

%%v data's PSD
figure;
[pxx_v, f_v, Pxxc_v] = pwelch(v, WINDOW, NOVERLAP, NFFT, Fs, 'ConfidenceLevel', 0.95);
loglog(f_v, pxx_v, 'g'); 
hold on;
confup = Pxxc_v(conf_level, 2) ./ pxx_v(conf_level) * confcen;
confdn = Pxxc_v(conf_level, 1) ./ pxx_v(conf_level) * confcen;
% loglog([f_v(conf_level), f_v(conf_level)], [confdn, confup], 'r', 'LineWidth', 1.5);
loglog(f_v(conf_level*[1 1]), Pxxc_v(conf_level,:), 'r', 'LineWidth', 1.5);
title('Power spectra of V at D5');
xlabel('Frequency (Hz)');
ylabel('Power/frequency ((cm/s)^2/Hz)');
grid on;
legend('PSD', '95% Confidence Interval');
hold off;

.

[kim]

It solved!

Always thank you for your help, thanks to your kindness.

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.