%
% 		 Course:	Biomedical Signal Processing
% Post-Graduate-Studies:        Biomedical Engineering-BΙM
%				University of Patras
%   		   Date:	05-06-2016
% 
%               Part 05:	Auto-correlation
%               		Cross-correlation
%				Mean Heart rate extimation
%				types of Noise: white Gaussian, white uniform


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    Auto-correlation Theory  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


clear all
close all
x= [ones(1,10), zeros(1,10) ] ;
A = xcorr(x) ; 
subplot(2,1,1), plot(x)
subplot(2,1,2), plot(A)
pause()
x= [x x x x x x x x x x x ] ;
A = xcorr(x) ; 
subplot(2,1,1), plot(x)
subplot(2,1,2), plot(A)
pause()
A = xcorr(x,'unbiased') ; 
subplot(2,1,1), plot(x)
subplot(2,1,2), plot(A)
pause()


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear all
close all
x= [ones(1,10), cos(2*pi*0.3*(1:20)) ] ;
x= [x x x x x x x x x x x ] ;
A = xcorr(x,'unbiased') ; 
subplot(2,1,1), plot(x)
subplot(2,1,2), plot(A)
pause()


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    Auto-correlation in noise - white & uniform  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear all
close all
x= [ones(1,10), cos(2*pi*0.3*(1:20)) ] ;
x= [x x x x x x x x x x x ] ;
subplot(3,1,1), plot(x)
xl = length(x) ;
x = x + rand(1,xl) ;
A = xcorr(x,'unbiased') ; 
subplot(3,1,2), plot(x)
subplot(3,1,3), plot(A)
pause()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    Auto-correlation in very high noise - white & uniform  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear all
close all
x= [ones(1,10), cos(2*pi*0.3*(1:20)) ] ;
x= [x x x x x x x x x x x ] ;
subplot(3,1,1), plot(x)
xl = length(x) ;
x = x + 5 * rand(1,xl) ;
A = xcorr(x,'unbiased') ; 
subplot(3,1,2), plot(x)
subplot(3,1,3), plot(A)
pause()
figure, plot(A(50:(end-50)))
pause()


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    Auto-correlation in ECG signals  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all
load ECGx1 ;
x = x - mean(x) ;
A = xcorr(x,'unbiased') ; 
subplot(2,1,1), plot(x)
subplot(2,1,2), plot(A)
pause()


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    Auto-correlation in noisy ECG signals  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all
load ECGx1 ;
x = x - mean(x) ;
x = x + 300 * (rand(1,length(x))-0.5) ; 
A = xcorr(x,'unbiased') ; 
subplot(2,1,1), plot(x)
subplot(2,1,2), plot(A)
B = fft(A) ;
B(1:3) = 0 ;
B(end-3:end) = 0 ; 
Bi = real(ifft(B));
figure, plot(Bi(1000:end-1000))
pause()


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    Cross-correlation in ECG signals  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all
load ECGxy ;
x1 = x(1,:) ;
x2 = x(2,:) ;
x1 = x1 - mean(x1) ;
A1 = xcorr(x1,'unbiased') ; 
x2 = x2 - mean(x2) ;
A2 = xcorr(x2,'unbiased') ; 
subplot(2,2,1), plot(x1)
subplot(2,2,2), plot(A1)
subplot(2,2,3), plot(x2)
subplot(2,2,4), plot(A2)
figure
Axy = xcorr(x1, x2, 'unbiased') ; 
subplot(2,1,1), plot(Axy)
Axy = Axy - mean(Axy) ;
AAxy = xcorr(Axy, 'unbiased') ; 
subplot(2,1,2), plot(AAxy(1000:end-1000))
pause()


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    Cross-correlation in very noisy ECG signals  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all
load ECGxy ;
x1 = x(1,:) ;
x2 = x(2,:) ;
x1 = x1 - mean(x1) ;
x1 = x1 + 1000*(rand(1,length(x1))-0.5) ;
A1 = xcorr(x1,'unbiased') ; 
x2 = x2 - mean(x2) ;
x2 = x2 + 1000*(rand(1,length(x1))-0.5) ;
A2 = xcorr(x2,'unbiased') ; 
subplot(2,2,1), plot(x1)
subplot(2,2,2), plot(A1)
subplot(2,2,3), plot(x2)
subplot(2,2,4), plot(A2)
figure
Axy = xcorr(x1, x2, 'unbiased') ; 
subplot(2,1,1), plot(Axy)
Axy = Axy - mean(Axy) ;
AAxy = xcorr(Axy, 'unbiased') ; 
subplot(2,1,2), plot(AAxy(1000:end-1000))
pause()


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    Cross-correlation in very noisy ECG signals: Noise reduction by FFT band-pass filtering  
%     Failed due to 'unbiased' 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all
load ECGxy ;
x1 = x(1,:) ;
x2 = x(2,:) ;
x1 = x1 - mean(x1) ;
x1 = x1 + 1000*(rand(1,length(x1))-0.5) ;
X1 = fft(x1) ;
X1(183:end-183) = 0 ;
xn1 = real(ifft(X1)) ;
A1 = xcorr(xn1,'unbiased') ; 

x2 = x2 - mean(x2) ;
x2 = x2 + 1000*(rand(1,length(x1))-0.5) ;
X2 = fft(x2) ;
X2(183:end-183) = 0 ;
xn2 = real(ifft(X2)) ;
A2 = xcorr(xn2,'unbiased') ; 
subplot(3,2,1), plot(x1)
subplot(3,2,2), plot(xn1)
subplot(3,2,3), plot(x2)
subplot(3,2,4), plot(xn2)
subplot(3,2,5), plot(A1(500:end-500))
subplot(3,2,6), plot(A2(500:end-500))
figure
Axy = xcorr(xn1, xn2, 'unbiased') ; 
subplot(2,1,1), plot(Axy)
Axy = Axy - mean(Axy) ;
AAxy = xcorr(Axy, 'unbiased') ; 
subplot(2,1,2), plot(AAxy(1000:end-1000))
pause()


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    Cross-correlation in very noisy ECG signals: Noise reduction by FFT band-pass filtering  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all
load ECGxy ;
x1 = x(1,:) ;
x2 = x(2,:) ;
x1 = x1 - mean(x1) ;
x1 = x1 + 1000*(rand(1,length(x1))-0.5) ;
X1 = fft(x1) ;
X1(183:end-183) = 0 ;
xn1 = real(ifft(X1)) ;
A1 = xcorr(xn1) ; 

x2 = x2 - mean(x2) ;
x2 = x2 + 1000*(rand(1,length(x1))-0.5) ;
X2 = fft(x2) ;
X2(183:end-183) = 0 ;
xn2 = real(ifft(X2)) ;
A2 = xcorr(xn2) ; 
subplot(3,2,1), plot(x1)
subplot(3,2,2), plot(xn1)
subplot(3,2,3), plot(x2)
subplot(3,2,4), plot(xn2)
subplot(3,2,5), plot(A1(500:end-500))
subplot(3,2,6), plot(A2(500:end-500))
figure
Axy = xcorr(xn1, xn2) ; 
subplot(2,1,1), plot(Axy)
Axy = Axy - mean(Axy) ;
AAxy = xcorr(Axy) ; 
subplot(2,1,2), plot(AAxy(1000:end-1000))
pause()


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    Cross-correlation in ECG signals in Huge noise : Noise reduction by FFT band-pass filtering  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all
close all
load ECGxy ;
x1 = x(1,:) ;
x2 = x(2,:) ;
x1 = x1 - mean(x1) ;
x1 = x1 + 5000*(rand(1,length(x1))-0.5) ;
X1 = fft(x1) ;
X1(183:end-183) = 0 ;
xn1 = real(ifft(X1)) ;
A1 = xcorr(xn1) ; 

x2 = x2 - mean(x2) ;
x2 = x2 + 5000*(rand(1,length(x1))-0.5) ;
X2 = fft(x2) ;
X2(183:end-183) = 0 ;
xn2 = real(ifft(X2)) ;
A2 = xcorr(xn2) ; 
subplot(3,2,1), plot(x1)
subplot(3,2,2), plot(xn1)
subplot(3,2,3), plot(x2)
subplot(3,2,4), plot(xn2)
subplot(3,2,5), plot(A1(500:end-500))
subplot(3,2,6), plot(A2(500:end-500))
figure
Axy = xcorr(xn1, xn2) ; 
subplot(2,1,1), plot(Axy)
Axy = Axy - mean(Axy) ;
AAxy = xcorr(Axy) ; 
subplot(2,1,2), plot(AAxy(1000:end-1000))
pause()


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    Uniform & Gaussian noise and spectrum distribution ....  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

close all
clear all
x = rand( 1, 100000 ) - 0.5 ;
X = fft(x) ;
plot(abs(X))

close all
clear all
x = randn( 1, 100000 );
X = fft(x) ;
plot(abs(X))

