close all;
% clear;
clc;

%% sampling parameters
% 
PSD_pwelch_setup_v2;
% e38_signal_gvw;
% e38_signal_half;
% signal = [E38_gvw; E38_half];
% PSD_pwelch_setup;
% signal = E38_gvw;

%% Evaluation X
[pxxX,fx] = pwelch(g*brr_x,window,overlap,nfft,fs,freqrange,"psd","psd");

% figure Name brr_x
% plot(t,brr_x)        
% xlabel('Time [s]');
% ylabel('Acceleration (g)');
% title('PumpBrac_F_X');

figure Name brr_x
plot(fx,pxxX)
xlabel('Frequency (Hz)')
ylabel('PSD ((m*s^2)^2/Hz)')
axis([5 120 0 0.75])

figure Name brr_x
subplot(2,1,1),plot(t,brr_x)        
xlabel('Time [s]');
ylabel('Acceleration (g)');
title('PumpBrac_F_X');
subplot(2,1,2),plot(fx,pxxX);
xlabel('Frequency (Hz)')
ylabel('PSD ((m*s^2)^2/Hz)')
axis([5 120 0 0.75])

%% Evaluation Y
[pxxY,fy] = pwelch(g*brr_y,window,overlap,nfft,fs,freqrange,"psd");

% figure Name brr_y
% plot(t,brr_y)        
% xlabel('Time [s]');
% ylabel('Acceleration (g)');
% title('PumpBrac_F_Y');

figure Name brr_y
plot(fy,pxxY)
xlabel('Frequency (Hz)')
ylabel('PSD ((m*s^2)^2/Hz)')
axis([5 120 0 0.75])

figure Name brr_y
subplot(2,1,1),plot(t,brr_y)        
xlabel('Time [s]');
ylabel('Acceleration (g)');
title('PumpBrac_F_Y');
subplot(2,1,2),plot(fy,pxxY);
xlabel('Frequency (Hz)')
ylabel('PSD ((m*s^2)^2/Hz)')
axis([5 120 0 0.75])

%% Evaluation Z
[pxxZ,fz] = pwelch(g*brr_z,window,overlap,nfft,fs,freqrange,"psd");

% figure Name brr_z
% plot(t,brr_z)        
% xlabel('Time [s]');
% ylabel('Acceleration (g)');
% title('PumpBrac_F_Z');

figure Name brr_z
plot(fz,pxxZ)
xlabel('Frequency (Hz)')
ylabel('PSD ((m*s^2)^2/Hz)')
axis([5 120 0 0.75])

figure Name brr_z
subplot(2,1,1),plot(t,brr_z)        
xlabel('Time [s]');
ylabel('Acceleration (g)');
title('PumpBrac_F_Z');
subplot(2,1,2),plot(fz,pxxZ);
xlabel('Frequency (Hz)')
ylabel('PSD ((m*s^2)^2/Hz)')
axis([5 120 0 0.75])


%% Graph output
figure
hold on
plot(fz,pxxX,'r')
plot(fz,pxxY,'g')
plot(fz,pxxZ,'b')
grid on
xlabel('Frequency (Hz)')
ylabel('PSD ((m*s^2)^2/Hz)')
title('bracket right')
legend('PSD bracket right x','PSD bracket right y','PSD bracket right z')
axis([5 120 0 0.75])

%% Parseval theorem X

LHSX = var(brr_x);
disp('displaying LHSX (signal)');
disp(LHSX);

RHSX = trapz(fx,pxxX);
disp('displaying RHSX (PSD)');
disp(RHSX);

PTX = abs(LHSX - RHSX);
disp('Parseval X');
disp(PTX);
%% Parseval theorem Y
LHSY = var(brr_y);
disp('displaying LHSY (signal)');
disp(LHSY);

RHSY = trapz(fy,pxxY);
disp('displaying RHSY (PSD)');
disp(RHSY);

PTY = abs(LHSY - RHSY);
disp('Parseval Y');
disp(PTY);
%% Parseval theorem Z
LHSZ = var(brr_z);
disp('displaying LHS (signal)');
disp(LHSZ);

energy_time_domain = (1/length(brr_z)) * sum(abs(brr_z).^2);
disp('displaying LHS (signal) energy');
disp(energy_time_domain);

rms_ev = sqrt(mean(brr_z.^2));
rms2 = rms_ev * rms_ev;

RHSZ = trapz(fz,pxxZ);
disp('displaying RHS (signal) rms');
disp(rms2);

RHSZ = trapz(fz,pxxZ);
disp('displaying RHS (PSD)');
disp(RHSZ);

PTZ = abs(LHSZ - RHSZ);
disp('Parseval Z');
disp(PTZ);

%% RMS

% RMS = rms (pxxX);
% disp('RMS');
% disp(RMS);

% A=fft(brr_z);
% figure
% plot(A)
% 
% figure
% plot(pxxZfft)
