%
% 		 Course:	Biomedical Signal Processing
% Post-Graduate-Studies:        Biomedical Engineering-BΙM
%				University of Patras
%   		   Date:	28-03-2016
% 
%               Part 01:        Noise and interference elimination 
%				Filters in Time and Frequency Domain
%				Filters representaion in Time, Frequency, z-plane domains
%				Filtering artificial signals
%				Examples using ECG signals
%				Examples using EEG signals
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    Sinusoidal elimination in the Time domain 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear all
close all

% two Simusoidal signal at 5 and 60Hz. Sampling rate at 1 KHz 	
x = sin(2*pi*(1:5000)*50/2048) ;
%figure, plot(x);

% Sinusoidal elimination using FIR Filters of two zeros
% Two zeros are defined at 60 Hz.
b = poly( [exp(2*pi*i*50/2048) exp(-2*pi*i*50/2048)]) 
a = poly( [0.99*exp(2*pi*i*50/2048) 0.99*exp(-2*pi*i*50/2048)]) 

figure;
freqz(b,a,2000,1000) ;
x1 = filter(b,a,x) ;

figure ;
plot(x)
hold on
plot(x1,'r')



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    DC and very low frequencies elimination in the Time domain, ECG examples  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all

% Load ECG signal x. Sampling rate at 360 Hz 	

load ECGx1

% Compute DC and remove it.
x1 = x - mean(x) ;
figure

plot(x);
hold on 
plot(x1,'r' );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    DC elimination by subtracting two sequencial values  

figure
freqz( [1 -1], 1, 2000, 360);

y = filter( [1 -1], 1, x) ;

figure
plot(x);
hold on 
plot(y,'r' );


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    DC and very low frequencies elimination in the Frequency domain, ECG examples 


% DC using Discete Fourier Transform
fx = fft(x) ;
figure, semilogy(abs(fx)) ;

% The DC componenet appears in the first position
fx(1)=0;

% Apply IFFT to the fx and plot the real part
x1=real(ifft(fx)) ;
figure ;
subplot(2,1,1), plot(x)
subplot(2,1,2), plot(x1)

% Elimination of the first 5 AC components i.e. from DC up to 5*360/3500= 0.5Hz
fx(2:6)=0;
fx((3500-4):3500)=0;
x2=real(ifft(fx)) ;
figure ;
subplot(3,1,1), plot(x)
subplot(3,1,2), plot(x1)
subplot(3,1,3), plot(x2)




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    Sinusoidal elimination in the frequency domain 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear all
close all

% two Simusoidal signal at 5 and 60Hz. Sampling rate at 1 KHz 	
x = sin(2*pi*(1:2000)*5/1000) + sin(2*pi*(1:2000)*60/1000) ;
figure, plot(x);
% pause

% Sinusoidal elimination using Discete Fourier Transform
fx = fft(x) ;
figure, semilogy(abs(fx)) ;
% pause

% Estimation of 60Hz position in the fx array: Freq * SignalSize / SamplingFrequency + 1 = 60*2000/1000 = 120
% Elimination of 60Hz components in two positions due to FFT symmetry
% The DC componenet appears in the first position ONLY, so 120+1 = 121
fx(121)=0;
% The DC componenet appears in the first position, so 120-1 = 119
fx(2000-119)=0;
figure, semilogy(abs(fx)) ;
% Apply IFFT to the fx and plot the real part
x1=real(ifft(fx)) ;
figure ;
subplot(2,1,1), plot(x)
subplot(2,1,2), plot(x1)





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    Sinusoidal elimination in the frequency domain 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear all
close all

% Load ECG signal x. Sampling rate at 360 Hz 	
load ECGnoise50

% Compute DC and remove it.
% x = x - mean(x) ;

xnoise=x ;
% Sinusoidal elimination using Discete Fourier Transform
fx = fft(x) ;
figure ;
axis([1 3500 10 100000]) ;
semilogy(abs(fx)) ;
% Two high values appears in 487 and symmetrically in 3015 array position of the FFT
% which correspond to (487-1)*360/3500 = 49.95Hz and (3500-3015-1)*360/3500=49.9Hz frequency
fx(487)=0;
fx(3015)=0;
x1=real(ifft(fx)) ;
figure ;
%subplot(2,1,1), plot(x)
%subplot(2,1,2), plot(x1)
% In the window limits small components of 50 Hz remains, due to windowing effect
hold on
load ECGx1
plot(x,'r') ;

% If greater number of frequencies are set to zero the 50 Hz interference is eliminated
fx((487-3):(487+3))=0;
fx((3015-3):(3015+3))=0;
x1=real(ifft(fx)) ;
figure ;
subplot(2,1,1), plot(x)
subplot(2,1,2), plot(x1)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    Sinusoidal elimination in the frequency domain, Windowing effect

% Improve filtering by windowing ....
clear all
close all

% Load ECG signal x. Sampling rate at 360 Hz 	
load ECGVeryHigh50

% Compute DC and remove it.
x = x - mean(x) ;
w = hamming(length(x))' ;
x1 = w.*x;
% Sinusoidal elimination using Discete Fourier Transform
fx1 = fft(x1) ;
axis([1 length(x1) 10 100000]) ;
semilogy(abs(fx1)) ;

fx1(475:500)=0;
fx1((3500-500-2):(3500-475-2))=0;
x11=real(ifft(fx1)) ;
x12 = x11 ./ w ;
figure ;
subplot(2,1,1), plot(x)
subplot(2,1,2), plot(x12)



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    Sinusoidal elimination in the time domain using FIR filters 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear all
close all

% Load ECG signal x. Sampling rate at 360 Hz 	
load ECGnoise50

% Compute DC and remove it.
x = x - mean(x) ;

%%% Estimation of FIR fiter Coefficients

b = poly( [ exp(2*pi*i*50/360) exp(-2*pi*i*50/360)] )
freqz(b,1,2500,360)

y = filter(b,1,x) ;
figure
subplot(2,1,1), plot(x)
subplot(2,1,2), plot(y)

