3DOF Decomposition

% HW7

clc,clear,close all

%% Signals

addpath('/Users/jahdielgierboliniserrano/Desktop/INGE5037/HW7/');
S1 = load('S1.txt')'; S2 = load('S2.txt')'; l = length(S1);

% time vector
sf = 100; dt = 1/sf; t = 0:dt:dt*(l-1);

% figures
subplot(211); plot(t,S1); ylabel('clean response');
title('Clean (above) and noisy (below) signals vs time');
subplot(212); plot(t,S2); ylabel('noisy response');

%% Part 1 (SNR)

noise = S2-S1; SNR = sum(S1.^2)/sum(noise.^2);
N = sqrt(var(S1)/SNR);
fprintf('SNR = %2.2f\n',SNR);
fprintf('The level of noise is %0.3f\n\n',N);

%% Part 2 (FFT)

[f1,m1,p1] = osfft(t,S1,0,0,'none'); % clean signal
[f2,m2,p2] = osfft(t,S2,0,0,'none'); % noisy signal

% figures
figure;
subplot(221); plot(f1,m1); ylabel('magnitude'); title('Clean Signal');
xlim([0 10]);
subplot(223); plot(f1,p1); ylabel('phase [rad]'); xlabel('freq. [Hz]');
xlim([0 10]);
subplot(222); plot(f2,m2); ylabel('magnitude'); title('Noisy Signal');
xlim([0 10]);
subplot(224); plot(f2,p2); ylabel('phase [rad]'); xlabel('freq. [Hz]');
xlim([0 10]);

%% Bandpass Filters and Hilbert Transform

% from inspection of the FFT Spectra
lba = 0; uba = 0.6/sf; lbb = uba; ubb = 1.3/sf; lbc = ubb; ubc = 2.7/sf;

% filter coefficients
ba = pbfilter(l,sf,lba,uba,500); bb = pbfilter(l,sf,lbb,ubb,500);
bc = pbfilter(l,sf,lbc,ubc,500);

% filtered signals
S1a = conv(S1,ba,'same'); S1b = conv(S1,bb,'same'); % clean
S1c = conv(S1,bc,'same'); S1f = [S1a' S1b' S1c']';

figure; subplot(211); plot(t,S1f); ylabel('amplitude'); title('Clean');

S2a = conv(S2,ba,'same'); S2b = conv(S2,bb,'same'); % noisy
S2c = conv(S2,bc,'same'); S2f = [S2a' S2b' S2c']';

subplot(212); plot(t,S2f); ylabel('amplitude'); xlabel('time [s]');
title('Noisy');

% Hilbert Transform
[Am1,~,F1] = HilbertRes(S1f,sf);
[Am2,~,F2] = HilbertRes(S2f,sf);

% Response Parameters
[zi1,wn1,A1,phi1,s1] = fr1dof(t,S1f,Am1,F1);
[zi2,wn2,A2,phi2,s2] = fr1dof(t,S2f,Am2,F2);

%% EMD and Hilbert Transform

%% Continuous Wavelet Transform