FFT Analysis,denoising And Filtering Of ECG Signal

Kshitij Singh answered . 2024-12-21 02:02:23

Function Plot ATM.m

function [val x] =  plotATM(Name)

% usage: plotATM('RECORDm')
%
% This function reads a pair of files (RECORDm.mat and RECORDm.info) generated
% by 'wfdb2mat' from a PhysioBank record, baseline-corrects and scales the time
% series contained in the .mat file, and plots them.  The baseline-corrected
% and scaled time series are the rows of matrix 'val', and each
% column contains simultaneous samples of each time series.
%
% 'wfdb2mat' is part of the open-source WFDB Software Package available at
%    http://physionet.org/physiotools/wfdb.shtml
% If you have installed a working copy of 'wfdb2mat', run a shell command
% such as
%    wfdb2mat -r 100s -f 0 -t 10 >100sm.info
% to create a pair of files ('100sm.mat', '100sm.info') that can be read
% by this function.
%
% The files needed by this function can also be produced by the
% PhysioBank ATM, at
%    http://physionet.org/cgi-bin/ATM
%

% plotATM.m           O. Abdala			16 March 2009
% 		      James Hislop	       27 January 2014	version 1.1

infoName = strcat(Name, '.info');
matName = strcat(Name, '.mat');
Octave = exist('OCTAVE_VERSION');
load(matName);
fid = fopen(infoName, 'rt');
fgetl(fid);
fgetl(fid);
fgetl(fid);
[freqint] = sscanf(fgetl(fid), 'Sampling frequency: %f Hz  Sampling interval: %f sec');
interval = freqint(2);
fgetl(fid);

if(Octave)
    for i = 1:size(val, 1)
       R = strsplit(fgetl(fid), char(9));
       signal{i} = R{2};
       gain(i) = str2num(R{3});
       base(i) = str2num(R{4});
       units{i} = R{5};
    end
else
    for i = 1:size(val, 1)
      [row(i), signal(i), gain(i), base(i), units(i)]=strread(fgetl(fid),'%d%s%f%f%s','delimiter','\t');
    end
end

fclose(fid);
val(val==-32768) = NaN;

for i = 1:size(val, 1)
    val(i, :) = (val(i, :) - base(i)) / gain(i);
end

x = (1:size(val, 2)) * interval;
%plot(x', val');

for i = 1:length(signal)
    labels{i} = strcat(signal{i}, ' (', units{i}, ')'); 
end

%legend(labels);
%%xlabel('Time (sec)');
% grid on

end

Main.m

%%
clc;
[val, x] = plotATM('100m');%%PlotAtm is function used for reading ecg  file 
%%
h=val(1,:);%% selection of ecg signal for further process
%%
Fs=360;%% sampling frequency given
t=(0:length(h)-1)/Fs;%% time for the signal
figure;
plot(t,h)%% plotting 
title('Original ECG signal in time domain')
%%
u=length(h);%% calculating length of signal
K=abs(fft(h));%%converting into frequency domain
fax_bins = [0 : u-1]; %frequency axis in bins
u_2 = ceil(u/2);
figure;
plot(fax_bins(1:u_2), K(1:u_2))
xlabel('Frequency (Bins)')
ylabel('Magnitude');
title('ECG in frequency domain signal (Single-sided Magnitude spectrum (bins))');
axis tight
ylim([0 200])
%%
%%
rd = 9;
fl = 21;

smtlb = sgolayfilt(h,rd,fl);

figure;
subplot(2,1,1)
plot(t,h)

%%

title('Original in time response')
grid

subplot(2,1,2)
plot(t,smtlb)

title('Filtered using  Savitzky-Golay in time response')
grid
%%
rd1 = 9;
fl2 = 21;

smtlb1 = sgolayfilt(K(1:u_2),rd1,fl2);%Smooth the signal by applying a Savitzky-Golay filter of polynomial order 9 to data frames of length 21.
figure;
subplot(2,1,1)
plot(K(1:u_2))
title('Original in frequency response')
grid

subplot(2,1,2)
plot(smtlb1)

title('Filtered in Savitzky-Golay frequency response')
grid
%% butterworth filter
o=5;%% coffiecient 
[b,a]=butter(o,[50 150]/360); %%[b,a] = butter(n,Wn) returns the transfer function coefficients of an nth-order lowpass digital Butterworth filter with normalized cutoff frequency Wn.
butter_filter=filter(b,a,K(1:u_2));
figure;
subplot(2,1,1)
title('Original EcG  in frequency domain')
plot(K(1:u_2))
subplot(2,1,2)
plot(butter_filter)
title('filtering using butterworth in frequency domain')
%%
o1=1;
[b1,a1]=butter(o1,[1 10]/180,'Bandpass');
butter_filter1=filter(b1,a1,h);
figure;
subplot(2,1,1)
plot(t,h)
subplot(2,1,2)
plot(t,butter_filter1)
title('filtering using butterworth in time domain')
%%

%%
for i=1:length(h)
SIR=h(i)/(h(i)-smtlb(i));
end
fprintf('Signal to interference when Savitzky-Golay is applied is--- %f\n',SIR);
%%
for i=1:length(h)
SIR2=h(i)/(h(i)-butter_filter1(i));
end
fprintf('Signal to interference when butterworth filter is applied is--- %f\n',SIR2);
%%
mse=0;
for i=1:length(h)
mse=mse+(smtlb(i)-h(i))^2;
end
mse=mse/length(h);
fprintf('mean squared error original ecg signal and filtered signal(Savitzky-Golay) %f\n',mse);
%%
%%
mse1=0;
for i=1:length(h)
mse1=mse1+(butter_filter1(i)-h(i))^2;
end
mse1=mse1/length(h);
fprintf('mean squared error original ecg signal and filtered signal(butterworth Filter) %f\n',mse1);


%%
num=0;
den=0;
for i=1:length(h)
den=den+(smtlb(i)-h(i))^2;
end
for i=1:length (h)
num=num+h(i)^2;
end
SNR = 20*log10(sqrt(num)/sqrt(den));
fprintf('signal to noise ratio between original ecg signal and filtered signal(Savitzky-Golay) %f db\n',SNR);
%%
num1=0;
den1=0;
for i=1:length(h)
den1=den1+(butter_filter1(i)-h(i))^2;
end
for i=1:length (h)
num1=num1+h(i)^2;
end
SNR = 20*log10(sqrt(num1)/sqrt(den1));
fprintf('signal to noise ratio between original ecg signal and filtered signal(butterworth filter) %f db\n',SNR);
%%

Not satisfied with the answer ?? ASK NOW