%% Feature Extraction and Detection - Part 1
% Goal: Identify spectral features in signals using wavelet techniques
%% Steps
% Load and visualize the signal
% Preprocess the signal to remove artifacts
% Perform wavelet based time-frequency analysis to identify features
close all
%% Load the signal
load('ecgSample.mat');
%% Visualize the signal
figure;
plot(t,ecgSample,'r');
grid on; title('ECG signal');
xlabel('time (sec)');
% OR you can use
%%
signalAnalyzer
%% Visualize the frequency components in the signal
figure;
pwelch(ecgSample,[],[],[],Fs); title('Power spectrum of signal');
%% Filter Signal
notchFilt = designfilt('bandstopiir', 'FilterOrder', 6, 'HalfPowerFrequency1', 59, ...
'HalfPowerFrequency2', 61, 'SampleRate', Fs);
ecgSampleFiltered = filtfilt(notchFilt, ecgSample);
%% Filtered Signal -Comparison
figure;
subplot(2,2,1)
helperTimeDomain(t,ecgSample, 'Original Signal',10,'r')
subplot(2,2,2)
helperFrequencyDomain(ecgSample,Fs,'Frequency Spectrum','r');
subplot(2,2,3)
helperTimeDomain(t,ecgSampleFiltered, 'Filtered Signal',10,'b')
subplot(2,2,4)
helperFrequencyDomain(ecgSampleFiltered,Fs,'Frequency Spectrum','b');
%% Further Exploration
% Let us see if we can use signal processing techniques to identify the QRS
% wave. To do this we can try and use the spectrogram technique which is a
% well known technique to examine the time frequency view of the
% signal.
%% %% Time-Frequency using Spectrogram
close all;
figure;
spectrogram(ecgSampleFiltered,[],[],[],Fs,'yaxis');
title('Time Frequency using Spectrogram');
figure;
plot(t,ecgSampleFiltered,'r');
grid on; title('Filtered ECG signal');
xlabel('time (sec)');
% <BACK TO PPT>
%% Time Frequency using CWT
figure;
ax1 = subplot(211);
helperTimeDomain(t./60,ecgSampleFiltered, 'Filtered Signal',[],'b')
ax2 = subplot(212);
cwt(ecgSampleFiltered,Fs);
% linkaxes([ax1,ax2],'x');
%% Exploring with CWT
figure;
cwt(ecgSampleFiltered,Fs);
xlim([0.02,0.08]);
[~,y1] = ginput(1);
[~,y2] = ginput(1);
fprintf('Frequency range of QRS wave is %0.2f Hz - %0.2f Hz\n',2.^min(y1,y2),2.^max(y1,y2));