%
% 		 Course:	Biomedical Signal Processing
% Post-Graduate-Studies:        Biomedical Engineering-BΙM
%				University of Patras
%   		   Date:	17-04-2016
% 
%               Part 03:	Design Filters in Time and Frequency Domain
%				FIR - Low-pass, band-pass, High-pass filters
%				IIR - Low-pass, band-pass, High-pass filters
%                               Estimating filter order
%				Examples using ECG signals
%				Examples using EEG signals
%



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

clear all
close all

f1 = 45/180 ;
f2 = 65/180 ;
b = fir1(100, [f1, f2], "stop");
freqz(b,1,1200,360);         	

%load ECGnoise50
load ECGx1
subplot(3,1,1)
plot(x);
x = x + 50 * cos( (1:3500)*2*pi*50/360) + 50 * cos( (1:3500)*2*pi*60/360) ;
subplot(3,1,2)
plot(x);
y = filter(b,1,x) ;
subplot(3,1,3)
plot(y) 



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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear all
close all

f = [0 3 40 45 65 70 80 100 180]/180 ;
b = fir2(100, f, [0 1 1 0 0 1 1 0 0] );

freqz(b,1,1200,360);         	
figure 
zplane(b,1);
figure

load ECGx1
subplot(3,1,1)
plot(x);
x = x + 50 * cos( (1:3500)*2*pi*50/360) + 50 * cos( (1:3500)*2*pi*60/360) ;
subplot(3,1,2)
plot(x);
y = filter(b,1,x) ;
subplot(3,1,3)
plot(y) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear all
close all

f = [0 3 40 45 65 70 80 100 180]/180 ;
b = fir2(100, f, [0 1 1 0.1 0.1 1 0.9 0.1 0] );

freqz(b,1,1200,360);         	
figure 
zplane(b,1);

figure
load ECGx1
subplot(3,1,1)
plot(x);
x = x + 50 * cos( (1:3500)*2*pi*50/360) + 50 * cos( (1:3500)*2*pi*60/360) ;
subplot(3,1,2)
plot(x);
y = filter(b,1,x) ;
subplot(3,1,3)
plot(y) 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% FIR filter design using least squares method. 
% Estimating FIR such that the integral of the weighted mean squared error in the specified bands is minimized. 
% I. Selesnick, "Linear-Phase FIR Filter Design by Least Squares," http://cnx.org/content/m10577 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all
   
f = [0    3  4 40  45    65     70 100 110 180]/180 ;
g = [0.01 1  1 1  0.001  0.001   1  1  0.1 0.1] ;
% w = [  1     1       2           1     1  ] ;
w = [   100    1       2           1     1  ] ;
b = firls(100, f,g,w );

freqz(b,1,1200,360);         	
figure 
plot(roots(b),'or');

figure
load ECGx1
subplot(3,1,1)
plot(x);
x = x + 50 * cos( (1:3500)*2*pi*50/360) + 50 * cos( (1:3500)*2*pi*60/360) ;
subplot(3,1,2)
plot(x);
y = filter(b,1,x) ;
subplot(3,1,3)
plot(y) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Estimate the minimum order FIR filter with the desired response characteristics. 
% [n,Wn,beta,ftype] = kaiserord(f,a,dev)
% kaiserord returns a filter order n and beta parameter to specify a Kaiser window for use with the fir1 function.
% Given a set of specifications in the frequency domain, kaiserord estimates the minimum FIR filter order that will
% approximately meet the specifications. 
% kaiserord converts the given filter specifications into passband and stopband ripples and converts cutoff frequencies
% into the form needed for windowed FIR filter design.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


clear all
close all

% Pass-Band [0 40]Hz, cut-off [45 180]
% Specify a passband ripple of 5% and a stopband attenuation of 40 dB (0.01).
[n,Wn,beta,ftype] = kaiserord([40 45],[1 0],[0.05 0.01],360);
b = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
freqz(b,1,1200,360);         	

% Design an odd-length band-pass filter. 
% Note that odd length means even order, so the input to fir1 must be an even integer.
% Pass-Band [2 40]Hz, cut-off less than 1 and greater than 45
% Passband ripple of 5% and a stopband attenuation of 40 dB.
[n,Wn,beta,ftype] = kaiserord([1 2 40 45],[0 1 0],[0.01 0.05 0.01],360);
b = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
freqz(b,1,1200,360);         	




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% IIR Butterworth filter design. 
% First-order (2) Butterworth Low-pass.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all

[b,a]=butter(2, 45/180);
%[b,a]=butter(10, 45/180);
%[b,a]=butter(10, 40/180);

freqz(b,a,1200,360);         	
figure 
zplane(b,a)

figure
load ECGx1
subplot(3,1,1)
plot(x);
x = x + 50 * cos( (1:3500)*2*pi*50/360) + 50 * cos( (1:3500)*2*pi*60/360) ;
subplot(3,1,2)
plot(x);
y = filter(b,a,x) ;
subplot(3,1,3)
plot(y) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% IIR Butterworth filter design. 
% Butterworth High-pass.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all

[b,a]=butter(2, 1/180, "high" );
%[b,a]=butter(10, 45/180);
%[b,a]=butter(10, 40/180);

freqz(b,a,1200,360);         	
figure 
zplane(b,a) ;

figure
load ECGx1
x = x + 50 * cos( (1:3500)*2*pi*50/360) + 50 * cos( (1:3500)*2*pi*60/360) ;
subplot(3,1,1)
plot(x);
y = filter(b,a,x) ;
subplot(3,1,2)
plot(y) 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Phase Correction ....
% Correction for phase distortion introduced by a one-pass filter
% Does square the magnitude response in the process.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

y1 = filtfilt(b,a,x); 
subplot(3,1,3)
plot(y1) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% What happend during phase correction ?
% After filtering the data in the forward direction, 
% filtfilt reverses the filtered sequence and runs it back through the filter. 
%The result has the following characteristics:
%    Zero-phase distortion
%    filter transfer function, which equals the squared magnitude of the original filter transfer function
%    filter order that is double the order of the filter specified by b and a
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% IIR Butterworth filter design. 
% Butterworth Stop-band.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all

[b,a]=butter(10, [45 65]/180, "stop" );

freqz(b,a,1200,360);         	
figure 
zplane(b,a) ;

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

% Phase Correction ....
y1 = filtfilt(b,a,x); 
subplot(3,1,3)
plot(y1) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% IIR Butterworth filter design. 
% Combining filters - Serial connection ...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all

% high-Pass filter
[b2,a2]=butter(4, 1/180, "high" );
% Stop-Band filter
[b1,a1]=butter(10, [45 65]/180, "stop" );

b = conv(b2,b1) ;
a = conv(a2,a1) ;

freqz(b,a,1200,360);         	
figure 
zplane(b,a) ;

figure
load ECGx1
x = x + 50 * cos( (1:3500)*2*pi*50/360) + 50 * cos( (1:3500)*2*pi*60/360) ;
subplot(3,1,1)
plot(x);
y = filter(b,a,x) ;
subplot(3,1,2)
plot(y) 
% Phase Correction ....
y1 = filtfilt(b,a,x); 
subplot(3,1,3)
plot(y1) 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Compute the minimum filter order of an butter filter with the desired response characteristics. 
% n = buttord (wp, ws, rp, rs)
% If wp and ws are scalars, then wp is the passband cutoff frequency and ws is the stopband edge frequency. 
% If ws is greater than wp, the filter is a low-pass filter. If wp is greater than ws, the filter is a high-pass filter.

% If wp and ws are vectors of length 2, then wp defines the passband interval and ws defines the stopband interval.
% If wp is contained within ws (ws1 < wp1 < wp2 < ws2), the filter is a band-pass filter.
% If ws is contained within wp (wp1 < ws1 < ws2 < wp2), the filter is a band-stop or band-reject filter. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

close all
clear all

% low-pass
Ord = buttord( 40/180, 45/180, 6, 40 ) 

% low-pass
Ord = buttord( 40/180, 41/180, 6, 40 ) 

% low-pass
Ord = buttord( 40/180, 41/180, 0.01, 40 ) 

% high-pass
Ord = buttord( 41/180, 40/180, 0.01, 40 ) 

% Band-pass
Ord = buttord( [45 65]/180, [40 70]/180, 6, 40 ) 

% Band-reject
Ord = buttord( [40 70]/180, [45 65]/180, 6, 40 ) 




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% IIR Chebyshev type I filter . 
% Low-pass design.
% Generate an Chebyshev type I filter with Rp dB of pass band ripple.
% 
% [b, a] = cheby1(n, Rp, Wc)
%    low pass filter with cutoff pi*Wc radians
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all

[b,a]=cheby1(5, 6, 40/180 );

freqz(b,a,1200,360);         	
figure 
zplane(b,a) ;

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

% Phase Correction ....
y1 = filtfilt(b,a,x); 
subplot(3,1,3)
plot(y1) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% IIR Chebyshev type I filter . 
% High-pass design.
% Generate an Chebyshev type I filter with Rp dB of pass band ripple.
% 
% [b, a] = cheby1(n, Rp, Wc, "high" )
%    High pass filter with cutoff pi*Wc radians
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all

[b,a]=cheby1(5, 6, 2/180, "high" );

freqz(b,a,1200,360);         	
figure 
zplane(b,a) ;

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

% Phase Correction ....
y1 = filtfilt(b,a,x); 
subplot(3,1,3)
plot(y1) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% IIR Chebyshev type I filter . 
% Pass-band design.
% Generate an Chebyshev type I filter with Rp dB of pass band ripple.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all

[b,a]=cheby1(5, 10, [2 40]/180 );

freqz(b,a,1200,360);         	
figure 
zplane(b,a) ;

figure
load ECGx1
x = x + 50 * cos( (1:3500)*2*pi*50/360) + 50 * cos( (1:3500)*2*pi*60/360) ;
subplot(3,1,1)
plot(x);
y = filter(b,a,x) ;
subplot(3,1,2)
plot(y) 

% Phase Correction ....
y1 = filtfilt(b,a,x); 
subplot(3,1,3)
plot(y1) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% IIR Chebyshev type I filter . 
% Stop-band design.
% Generate an Chebyshev type I filter with Rp dB of pass band ripple.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all

[b,a]=cheby1(5, 10, [40 60]/180, "stop" );

freqz(b,a,1200,360);         	
figure 
zplane(b,a) ;

figure
load ECGx1
x = x + 50 * cos( (1:3500)*2*pi*50/360) + 50 * cos( (1:3500)*2*pi*60/360) ;
subplot(3,1,1)
plot(x);
y = filter(b,a,x) ;
subplot(3,1,2)
plot(y) 

% Phase Correction ....
y1 = filtfilt(b,a,x); 
subplot(3,1,3)
plot(y1) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Compute the minimum filter order of an cheby1 filter with the desired response characteristics. 
% n = cheb1ord (wp, ws, rp, rs)
% If wp and ws are scalars, then wp is the passband cutoff frequency and ws is the stopband edge frequency. 
% If ws is greater than wp, the filter is a low-pass filter. If wp is greater than ws, the filter is a high-pass filter.

% If wp and ws are vectors of length 2, then wp defines the passband interval and ws defines the stopband interval.
% If wp is contained within ws (ws1 < wp1 < wp2 < ws2), the filter is a band-pass filter.
% If ws is contained within wp (wp1 < ws1 < ws2 < wp2), the filter is a band-stop or band-reject filter. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

close all
clear all

% low-pass
Ord = cheb1ord( 40/180, 45/180, 6, 40 ) 

% low-pass
Ord = cheb1ord( 40/180, 41/180, 6, 40 ) 

% low-pass
Ord = cheb1ord( 40/180, 41/180, 0.01, 40 ) 

% high-pass
Ord = cheb1ord( 41/180, 40/180, 0.01, 40 ) 

% Band-pass
Ord = cheb1ord( [41 69]/180, [40 70]/180, 6, 40 ) 

% Band-reject
Ord = cheb1ord( [40 70]/180, [45 65]/180, 6, 40 ) 






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% IIR Chebyshev type II filter . 
% Low-pass design.
% Generate a Chebyshev type II filter with rs dB of stopband attenuation
% [b, a] = cheby2 (n, rs, wc)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all

[b,a]=cheby2(5, 40, 40/180 );

freqz(b,a,1200,360);         	
figure 
zplane(b,a) ;

figure
load ECGx1
x = x + 50 * cos( (1:3500)*2*pi*50/360) + 50 * cos( (1:3500)*2*pi*60/360) ;
subplot(3,1,1)
plot(x);
y = filter(b,a,x) ;
subplot(3,1,2)
plot(y) 

% Phase Correction ....
y1 = filtfilt(b,a,x); 
subplot(3,1,3)
plot(y1) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% IIR Chebyshev type II filter . 
% Pass-Band design.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all

[b,a]=cheby2(5, 40, [2 40]/180 );

freqz(b,a,1200,360);         	
figure 
zplane(b,a) ;

figure
load ECGx1
x = x + 50 * cos( (1:3500)*2*pi*50/360) + 50 * cos( (1:3500)*2*pi*60/360) ;
subplot(3,1,1)
plot(x);
y = filter(b,a,x) ;
subplot(3,1,2)
plot(y) 

% Phase Correction ....
y1 = filtfilt(b,a,x); 
subplot(3,1,3)
plot(y1) 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% IIR Chebyshev type II filter . 
% Stop-Band design.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all

[b,a]=cheby2(5, 40, [40 60]/180, "stop" );

freqz(b,a,1200,360);         	
figure 
zplane(b,a) ;

figure
load ECGx1
x = x + 50 * cos( (1:3500)*2*pi*50/360) + 50 * cos( (1:3500)*2*pi*60/360) ;
subplot(3,1,1)
plot(x);
y = filter(b,a,x) ;
subplot(3,1,2)
plot(y) 

% Phase Correction ....
y1 = filtfilt(b,a,x); 
subplot(3,1,3)
plot(y1) 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Compute the minimum filter order of an cheby2 filter with the desired response characteristics. 
% n = cheb2ord (wp, ws, rp, rs)
% If wp and ws are scalars, then wp is the passband cutoff frequency and ws is the stopband edge frequency. 
% If ws is greater than wp, the filter is a low-pass filter. If wp is greater than ws, the filter is a high-pass filter.

% If wp and ws are vectors of length 2, then wp defines the passband interval and ws defines the stopband interval.
% If wp is contained within ws (ws1 < wp1 < wp2 < ws2), the filter is a band-pass filter.
% If ws is contained within wp (wp1 < ws1 < ws2 < wp2), the filter is a band-stop or band-reject filter. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

close all
clear all

% low-pass
Ord = cheb2ord( 40/180, 45/180, 6, 40 ) 

% low-pass
Ord = cheb2ord( 40/180, 41/180, 6, 40 ) 

% low-pass
Ord = cheb2ord( 40/180, 41/180, 0.01, 40 ) 

% high-pass
Ord = cheb2ord( 41/180, 40/180, 0.01, 40 ) 

% Band-pass
Ord = cheb2ord( [41 69]/180, [40 70]/180, 6, 40 ) 

% Band-reject
Ord = cheb2ord( [40 70]/180, [45 65]/180, 6, 40 ) 




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% IIR Elliptic or Cauer filter 
% Generate an elliptic filter with RP dB of passband ripple and RS dB of stopband attenuation.
% low-pass design.
% [b, a] = ellip (n, rp, rs, wp)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


clear all
close all

[b,a]=ellip(5, 6, 40, 40/180 );

freqz(b,a,1200,360);         	
figure 
zplane(b,a) ;

figure
load ECGx1
x = x + 50 * cos( (1:3500)*2*pi*50/360) + 50 * cos( (1:3500)*2*pi*60/360) ;
subplot(3,1,1)
plot(x);
y = filter(b,a,x) ;
subplot(3,1,2)
plot(y) 

% Phase Correction ....
y1 = filtfilt(b,a,x); 
subplot(3,1,3)
plot(y1) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% IIR Elliptic or Cauer filter 
% Generate an elliptic filter with RP dB of passband ripple and RS dB of stopband attenuation.
% Band-pass design.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


clear all
close all

[b,a]=ellip(5, 6, 40, [2 40]/180 );

freqz(b,a,1200,360);         	
figure 
zplane(b,a) ;

figure
load ECGx1
x = x + 50 * cos( (1:3500)*2*pi*50/360) + 50 * cos( (1:3500)*2*pi*60/360) ;
subplot(3,1,1)
plot(x);
y = filter(b,a,x) ;
subplot(3,1,2)
plot(y) 

% Phase Correction ....
y1 = filtfilt(b,a,x); 
subplot(3,1,3)
plot(y1) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% IIR Elliptic or Cauer filter 
% Generate an elliptic filter with RP dB of passband ripple and RS dB of stopband attenuation.
% Stop-Band design.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


clear all
close all

[b,a]=ellip(5, 6, 40, [40 60]/180, "stop" );

freqz(b,a,1200,360);         	
figure 
zplane(b,a) ;

figure
load ECGx1
x = x + 50 * cos( (1:3500)*2*pi*50/360) + 50 * cos( (1:3500)*2*pi*60/360) ;
subplot(3,1,1)
plot(x);
y = filter(b,a,x) ;
subplot(3,1,2)
plot(y) 

% Phase Correction ....
y1 = filtfilt(b,a,x); 
subplot(3,1,3)
plot(y1) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Compute the minimum filter order of an elliptic filter with the desired response characteristics. 
% n = ellipord (wp, ws, rp, rs)
% If wp and ws are scalars, then wp is the passband cutoff frequency and ws is the stopband edge frequency. 
% If ws is greater than wp, the filter is a low-pass filter. If wp is greater than ws, the filter is a high-pass filter.

% If wp and ws are vectors of length 2, then wp defines the passband interval and ws defines the stopband interval.
% If wp is contained within ws (ws1 < wp1 < wp2 < ws2), the filter is a band-pass filter.
% If ws is contained within wp (wp1 < ws1 < ws2 < wp2), the filter is a band-stop or band-reject filter. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

close all
clear all


% low-pass
Ord = ellipord( 40/180, 45/180, 6, 40 ) 

% low-pass
Ord = ellipord( 40/180, 41/180, 6, 40 ) 

% low-pass
Ord = ellipord( 40/180, 41/180, 0.01, 40 ) 

% high-pass
Ord = ellipord( 41/180, 40/180, 0.01, 40 ) 

% Band-pass
Ord = ellipord( [41 69]/180, [40 70]/180, 6, 40 ) 

% Band-reject
Ord = ellipord( [40 70]/180, [45 65]/180, 6, 40 ) 



