Matlab code to plot the FFT of the windowed segments of ECG signal

ECE/BIOM 537: Biomedical Signal Processing

Colorado State University

Student: Minh Anh Nguyen

minhanhnguyen@q.com

• 1Split your assigned ECG waveform into individual PQRST segments. Multiply each segment by a window. Plot the FFT of the windowed segments and inspect a couple to make sure the windows make sense (you should include one or two of these plots in your report). Once the segments are properly windowed, line them up if needed and average all the FFT segments together, and include this in your report.

UPDATE: Compare this against  the full  length Fourier transform. Make comments about both spectra about the features.

-         The FFT of the windowed or one segment has less noise than the full length Fourier transform.

-         The dominant peak value of the full length Fourier transforms has higher values (.0534 dB) than FFT of one segment 0.2002dB. The spike in the frequency spectrum corresponds to dominant of frequency is 4.5 Hz in the full length Fourier transform while the dominant of frequency of the FFT of one segment is 3.9 Hz.

-         When I multiple each segment by a window, the ECG signal flip; therefore the fft result is different from the original ECG signal. The dominant peak value is .0092dB and the spike frequency spectrum is 4.5 Hz.

-         I multiple one single PQRST segment with the window to compare confirm my result. I noted that the ECG signal of multiple a segment with window is flipped. For this reason, I assumed that my output for multiple all segment with window is correct.

• Create a 4 Hz sinusoidal waveform (ie, a waveform of the form A*sin(w0*t) where w0 = 2*pi*f, A is amplitude, and f is the sinusoidal frequency) and simulate the waveform for some time. Combine that with a 1 Hz sinusoid by amplitude modulation, and plot this. Also, plot the frequency spectrum of the new waveform and detect the peaks. 

Remember, AM modulation can be achieved by A*sin(w0*t) where A = B*cos(w*t), instead of a scalar.

 

close all;

clear all;
clc;
fs = 250 % find the sampling rate or frequency
T = 1/fs;% sampling rate or frequency
window = 120; % 2 min 0r 120 second
fignum = 0;
load('J:\BIOM_Signal_processing\exam1\ECGsignalTest_1.mat') % contains hr_sig and fs
% Make time axis for ECG signal
tx = [0:length(hr_sig)-1]/fs;
fignum = fignum + 1;
figure(fignum)
%subplot (211)
plot(tx,hr_sig)
xlabel('Time (s)')
ylabel('Amplitude (mV)')
title('Zoom into ECG signal')
xlim([30.3,31]) % Used to zoom in on single ECG waveform
%load('\\tsclient\E\BIOM_Signal_processing\exam1\ECGsignalTest_1.mat')
%disp('Contents of workspace after loading file:')
%whos
%% copy the data and put into excel
y1=xlsread('J:\BIOM_Signal_processing\exam1\ecg_1.xls');
% find the length of the data per second
N = length(y1); %length of signal
ls = size(y1);
t = (0 : N-1) / fs;% sampling period or time vector
fignum = fignum + 1; %% keep track of figures
figure(fignum)
subplot (211), plot(t,y1);
title ('plot of the original of ECG signal')
xlabel ('time (sec)')
ylabel ('Amplitute (mv)')
grid on;
%Compute the spectrum of the ECG and provide remarks on the spectral features of the ECG ( see reference “ECG Statistics, Noise, Artifacts, and Missing Data”).
%%%  DFT to describe the signal in the frequency
NFFT = 2 ^ nextpow2(N);
Y = fft(y1, NFFT) / N;
f = (fs / 2 * linspace(0, 1, NFFT / 2+1))'; % Vector containing frequencies in Hz
amp = ( 2 * abs(Y(1: NFFT / 2+1))); % Vector containing corresponding amplitudes
%figure;
subplot(212), plot (f, amp);
title ('plot single-sided amplitude spectrume of the ECG signal');
xlabel ('frequency (Hz)');
ylabel ('|y(f)|');
grid on;

%%% Create a period
figure;
y1new = y1(180:300);
N1 = length(y1new);
%hold off
subplot (2,1,1), plot (y1new);
title ('plot one  typical PQRST segments of the signal amplitude spectrume of the orginal ECG signal')
xlabel ('time (msec)')
ylabel ('Amplitute  (mv)')
grid on;

%%%  DFT to describe the signal in the frequency of 1 window
NFFT1 = 2 ^ nextpow2(N1); % FFT transform length
Y2 = fft(y1new, NFFT1) / N1; % FFT, Discard unnecessary part of spectrum and scale amplitude
f1 = (fs / 2 * linspace(0, 1, NFFT1 / 2+1))'; % Vector containing frequencies in Hz
amp1 = ( 2 * abs(Y2(1: NFFT1 / 2+1))); % Vector containing corresponding amplitudes
%figure;
subplot (2,1,2), plot (f1, amp1);
title ('single-sided amplitude spectrume of one PQRST segments of ECG signal');
xlabel ('frequency (Hz)');
ylabel ('|y(f)|');
grid on;

xx=zeros(1,length(y1));
Xtwz = zeros(1,N); % pre-allocate STFT output array
M = length(y1new);           % M = window length, N = FFT length
zp = zeros(N-M,1);       % zero padding (to be inserted)
xoff = 0;                % current offset in input signal x
Mo2 = (M-1)/2;           % Assume M odd for simplicity here
for m=1:N
    xt = y1(xoff+1:xoff+M); % extract frame of input data
    xtw = y1new.* xt;         % apply window to current frame
    xtwz = [xtw(Mo2+1:M); zp; xtw(1:Mo2)]; % windowed, zero padded
    Xtwz(m,:) = fft(xtwz); % STFT for frame m
    %xoff = xoff + R;       % advance in-pointer by hop-size R
end

figure;
plot (xtw);

figure;
plot (Xtwz);




%% window
% Stride through the input with 50% overlap
% calculate windowed length “fftlen” FFT of each block,
% then inverse transform and overlap-add.
for i=1:1024/2:(length(y1)-1024)
%ff=fft(y1new.*y1(i:(i+128-1))',128);
%
% User processing of FFT data would go here...
%
%xx(i:(i+NFFT-1))= xx(i:(i+NFFT-1))+ ifft(ff,NFFT);
end

%% create a sub-sample of original ECG
y1_sub = y1(5000:7600);
tsub = 1:length(y1_sub);
figure;
subplot (2,2,1), plot(tsub,y1_sub);
title ('locate ST-interval (Angina pectoris) of the original ECG signal')
xlabel ('time (msec)')
ylabel ('Amplitute (mv)')
grid on

%% create a sub-sample of original ECG for 1 ST
y1_1 = y1(5042:5250);
tsub_1 = 1:length(y1_1);
N1_1=length(y1_1);
%figure;
subplot (2,2,2), plot(tsub_1,y1_1);
title ('Zoom-in ST-interval (Angina pectoris) of the original ECG signal')
xlabel ('time (msec)')
ylabel ('Amplitute (mv)')
grid on


%Compute the spectrum of Angina pectoris the ECG and provide remarks on the spectral features of the ECG .
%%%  DFT to describe the signal in the frequency
NFFT1_1 = 2 ^ nextpow2(N1_1);
Y2_1 = fft(y1_1, NFFT1_1) / N1_1;
f1_1 = (fs / 2 * linspace(0, 1, NFFT1_1 / 2+1))'; % Vector containing frequencies in Hz
amp1_1 = ( 2 * abs(Y2_1(1: NFFT1_1 / 2+1))); % Vector containing corresponding amplitudes
%figure;
subplot (2,2,3), plot (f1_1, amp1_1);
title ('single-sided amplitude spectrume of Angina pectoris segments of ECG signal');
xlabel ('frequency (Hz)');
ylabel ('|y(f)|');
grid on;

%Create a 4 Hz sinusoidal waveform (ie, a waveform of the form A*sin(w0*t) where w0 = 2*pi*f, A is amplitude, and f is the sinusoidal frequency)

clc;
clear all;
close all;

fm =1;
f= 4; %4 Hz sinusoidal
w0 = 2*pi*f;
Ts = 1/f;                     % sampling period
Tmax = 2.0;                    % signal duration
B= 1;
ka=1; %modulation coefficient

%t = [0:Ts:Tmax];               % time vector
t=0:0.001:1;
set(0,'defaultlinelinewidth',2);
A = B*cos(2*pi*f*t);
S=A.*sin(w0*t);%waveform Signal

AM=A.*(1+ka*S); %Amplitude Modulated wave

subplot(3,1,1);%Plotting frame divided in to 3 rows and this fig appear at 1st
plot(t,S);
xlabel('Time');
ylabel('Amplitude');
title('A*sin(2*pi*f*t)waveform');
grid on;

subplot(3,1,2); %plotting the  wave
plot(t,A);
xlabel('Time');
ylabel('Amplitude');
title('A = B*cos(w*t) waveform');
subplot(3,1,3); %plotting the amplitude modulated wave
plot(t,AM);
xlabel('Time');
ylabel('Amplitude');
title('AM modulation signal');


%Compute the spectrum of Angina pectoris the ECG and provide remarks on the spectral features of the ECG .
%%%  DFT to describe the signal in the frequency
N1 = length (S);
NFFT1 = 2 ^ nextpow2(N1);
Y2 = fft(S, NFFT1) / N1;
f1 = (f / 2 * linspace(0, 1, NFFT1 / 2+1))'; % Vector containing frequencies in Hz
amp1 = ( 2 * abs(Y2(1: NFFT1 / 2+1))); % Vector containing corresponding amplitudes
figure;
 subplot(211), plot(f1, amp1);
title ('single-sided amplitude spectrume of  A*sin(2*pi*f*t)waveform ');
xlabel ('frequency (Hz)');
ylabel ('|y(f)|');
grid on;

%Compute the spectrum of Angina pectoris the ECG and provide remarks on the spectral features of the ECG .
%%%  DFT to describe the signal in the frequency
N1_1 = length (AM);
NFFT1_1 = 2 ^ nextpow2(N1_1);
Y2_1 = fft(AM, NFFT1_1) / N1_1;
f1_1 = (f / 2 * linspace(0, 1, NFFT1_1 / 2+1))'; % Vector containing frequencies in Hz
amp1_1 = ( 2 * abs(Y2_1(1: NFFT1_1 / 2+1))); % Vector containing corresponding amplitudes
%figure;
 subplot(212); plot(f1_1, amp1_1);
title ('single-sided amplitude spectrume of AM modulation signal ');
xlabel ('frequency (Hz)');
ylabel ('|y(f)|');
grid on;