%
% 		 Course:	Biomedical Signal Processing
% Post-Graduate-Studies:        Biomedical Engineering-BΙM
%				University of Patras
%   		   Date:	04-04-2016
% 
%               Part 02:        Narrow and Wide-band Noise 
%				Filters in Time and Frequency Domain
%				Filters representaion in Time, Frequency, z-plane domains
%				Filtering artificial signals
%				Low-pas, band-pass, High-pass filters
%				Examples using ECG signals
%				Examples using EEG signals
%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    Sinusoidal elimination in the time domain using IIR 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 IIR fiter Coefficients

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

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


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

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

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

%%% Estimation of IIR fiter Coefficients

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

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


%%% Estimation of IIR filter Coefficients

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

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


%%% Effect of Coefficients Accuracy

b1 = floor(100*b)/100
a1 = floor(100*a)/100
figure
freqz(b,a,2500,360)
figure
freqz(b1,a1,2500,360)

y = filter(b,a,x) ;
y1 = filter(b1,a1,x) ;
figure
subplot(2,1,1), plot(y)
subplot(2,1,2), plot(y1)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    FIR filters, FFT filtering, FIR filtering using FFT  
%    Examples using artificial signals ...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear all
close all

% Impulse response
h = [2 -3 4 2 ] ;
% Signal
x = [1 1 3 4 1 ];
% Filtering using convolution
y = conv(h,x)

% Filtering using FFT
L=length(h)+length(x)-1;
H = [h zeros(1,L-length(h))]
X = [x zeros(1,L-length(x))]
Hv = fft(H);
Xv = fft(X);
% Multiplication in the frequency domain
Yv=Hv.*Xv ;
yn = real(ifft(Yv)) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    FIR filter design, low, high, stop, pass, bandpass
%    Examples using ECG signals ...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all
fc = 50/180 ;  		% Specifications Low pass at 50 Hz, Sampling frequency at 360 Hz
b = fir1(5, fc); 	% FIR 5th order filter coefficients
freqz(b,1,1200,360);  
figure
zplane(roots(b)) ;       	

b = fir1(10, fc); 	% FIR 10th order  filter coefficients
figure
freqz(b,1,1200,360);         	
figure
zplane(roots(b)) ;       	


b = fir1(20, fc); 	% FIR 20th order filter coefficients
figure
freqz(b,1,1200,360);         	
figure
zplane(roots(b)) ;       	

b = fir1(100, fc); 	% FIR 100th order filter coefficients
figure
freqz(b,1,1200,360);         	
figure
zplane(roots(b)) ;       	


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

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

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



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    FIR filter design, Estimating impulse response 

clear all
close all
fc = 50/180 ;  		% Specifications Low pass at 50 Hz, Sampling frequency at 360 Hz
b = fir1(100, fc); 	% FIR 5th order filter coefficients
freqz(b,1,1200,360);         	

h = filter(b,1,[1 zeros(1,100)]) ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    FIR filtering using FFT  
%    Examples using ECG signals ...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Filtering using convolution
load ECGnoise50

y = conv(h,x) ;
figure
subplot(3,1,1)
plot(x), 
subplot(3,1,2)
plot(y(1:3500),'r' );

% Filtering using FFT
L=length(h)+length(x)-1;
H = [h zeros(1,L-length(h))] ;
X = [x zeros(1,L-length(x))] ;
Hv = fft(H);
Xv = fft(X);
% Multiplication in the frequency domain
Yv=Hv.*Xv ;
yn = real(ifft(Yv)) ;
subplot(3,1,3)
plot(yn(1:3500),'r' );


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    FIR pass-band filtering   
%    Examples using ECG signals ... 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

f1 = 45/180 ;
f2 = 65/180 ;
% FIR 200 coefficients ... 
b = fir1(200, [f1, f2], "pass");
freqz(b,1,1200,360);         	

h = filter(b,1,[1 zeros(1,20)]) ;

% Filtering using convolution
%load ECGnoise50
load ECGVeryHigh50

y = conv(h,x) ;
figure
subplot(3,1,1)
plot(x), 
subplot(3,1,2)
plot(y(1:3500),'r' );
subplot(3,1,3)
plot(y(1:3500)-x,'r' );



