# Spectral Representation Digital Signals

Sampling, Aliasing, and Spectral Representation of Digital Signals

Part 1: Simulated Signals

1. Generate a sin signal consisting of 5 Hz, 30 Hz, 50 Hz, 64 Hz, and 70 Hz frequency components using the following specifications:
2. Sampling frequency of 256 Hz;
3. Sampling time that would generate a signal of 1024 data points;
4. Signal amplitude at 5 Hz of 2 mV; at 30 Hz of 1 mV; at 50 Hz of 2 mV; 64 Hz at 4 mV; at 70 Hz of 3 mV.
1. Plot the signal as a function of time in the expanded view (0 to 0.5 seconds).
1. Compute and display the magnitude spectrum of your signal. Use the following specifications and recommended functions:
1. Using the DSP System Toolbox FFT object (not the fft.m function from the Signal Processing toolbox!), compute the magnitude spectrum of the given signal and plot as a function of cyclic frequency. Be sure to appropriately scale the magnitudes such that the values correspond to the true magnitudes of the given signal components.
2. Make sure to scale your results appropriately to display the peak magnitude of your signal. Limit the frequency axis of your plot to display frequencies between 0 and 75 Hz.
1. Re-sample your signal with a sampling frequency of 128 Hz.
• Do not change the frequency components and their respective amplitudes.
1. Plot the signal as a function of time in the expanded view (0 to 0.5 seconds).
1. Compute and display the magnitude spectrum of your signal using the FFT object function.
1. Re-sample your signal with a sampling frequency of 128 Hz using the following specifications:
• Introduce a phase of to 64 Hz frequency component in your signal.  Do not change any other parameters.
1. Plot the signal as a function of time in the expanded view (0 to 0.5 seconds).
1. Compute and display magnitude spectrum of your signal.

Part II: Biomedical Data Analysis

You are given an Electroencephalogram (EEG) data file EEG_tbi.txt obtained from a male patient who suffered a non-acute traumatic brain injury (TBI) as a result of a blast force in a non-disclosed military location in Iraq (data courtesy of Jordan Neuroscience Inc.).  As a result of the injury, the patient was experiencing ongoing seizures, as well as other clinical symptoms consistent with a severe non-acute TBI.  As part of the medical intervention, the patient was placed in a medically induced coma and transported to the military hospital in Germany for further treatment and observation.  The patient remained in coma, and the continuous EEG recordings were conducted over several days to monitor the patient’s neurological status.  To maintain the patient’s lung recruitment while comatose, the patient was artificially ventilated using a High Frequency Oscillation Ventilation (HFOV) system.  The HFOV system is designed to deliver small tidal volumes (1-4mL/kg) at frequencies in the range of 3-15 Hz.

The typical settings for an adult patient include a bias flow of 40 L/min, an inspiratory time of 33%, and the frequency of operation is based on most recent arterial blood gas: pH<7.10 = 4 Hz, pH 7.10-7.19 = 5 Hz, pH 7.20-7.35 = 6 Hz, and pH>7.35 = 7 Hz (Fessler, H.E. at al., “A protocol for high-frequency oscillatory ventilation in adults: results from a roundtable discussion.” Crit. Care Med. 2007 Jul; 35(7):1649-54. PubMed PMID: 17522576).

The operating frequency range of the HFOV system overlaps with the frequency range of a typical non-convulsive EEG seizure (3-8 Hz), and could potentially make the clinical assessment of EEG in severe TBI patients difficult.  The HFOV artifact, however, would not typically affect all EEG channels equally.  As such, the identification of the affected channel/channels is important, as this channel could simply be excluded from the assessment.  If exclusion is not possible (e.g. too many affected channels), a digital filter could be designed to minimize the contribution of the HFOV artifact to the signal while maintaining all relevant EEG activity as much as possible.  One of the most effective methods for the identification of the instrumentation-induced frequencies in the signal of interest is the Fast Fourier Transform (FFT) spectral estimation.   Provided that the frequency resolution is acceptable, this technique would identify the sharp frequency peaks associated with the instrumentation noise, such as that of HFOV.

The EEG signal was recorded using the standard 10/20 EEG system with standard silver-chloride surface electrodes.  The data were acquired at a frequency of 200 Hz.  Prior to sampling, the signals were bandlimited to 75 Hz.  The signals represent the time varying voltages in units of micro-Volts.  The EEG signals recorded during this session are stored in columns/channels 4-21.

Data Analysis

Load the EEG data into MATLAB and plot as a function of time.  Use the Time Scope to display your data over multiple channels.  Observe the nature of the signals.  To improve the display quality, adjust the size of the axes for optimal visualization.  Evaluate each of the 4-21 channels visually.  Can you tell which of the channels are affected by the HFOV artifact without further processing?

Using the FFT object function from Part I, compute and plot the magnitude spectra of the given EEG signals.  Visually inspect the spectra and determine which of the given channels are affected by the HFOV artifact.  Identify three channels that are affected by the artifact most strongly. Make sure to plot those channels.

Oral Report

Each student is required to meet with Dr. Imas in lab in Week 7 to present and discuss their findings from laboratory.  In preparation for this discussion, please put together a document discussed below.  Bring the printed document to your meeting with Dr. Imas.  The document must include the following.

1. All required figures with appropriate figure numbers and legends. Include the figures under their corresponding Part 1 and Part 2 headers.
2. Under each figure, include the corresponding MATLAB code used to generate it.
3. If computations of specific numerical parameters were performed, organize the computed parameters in the table/s, label the table/s appropriately and bring with you to the meeting.

Recommended Analysis Questions

Part I: Simulated Signals

1. Given the relevant frequency content of your signal (peak frequencies excluding noise), what sampling frequencies were acceptable for accurate frequency domain representation of your signal? Explain your answer clearly.
1. Hopefully, you were able to identify appropriate sampling frequencies. If more than one appropriate sampling frequency was identified, which one would you rather select, i.e. which one is optimal?  Why?  Provide a reason of why the rejected frequencies may be chosen instead?
1. Which one of the sampling frequencies used resulted in aliased frequency spectrum? Explain your answer. Identify all aliased frequencies, and the frequencies that they originated from.
1. Assume that prior to sampling your signal, you analog lowpass-filtered it to contain 10 and 25 Hz components only. You selected your sampling frequency to be 50 Hz.  Based on your laboratory results, is that an advisable selection?  Why or why not?  Make sure to comment on two cases (signal with and without phase of π/2).

Part II: Biomedical Data Analysis

1. Were you able to identify the noise frequency of the HFOV system? What is the frequency?
1. You were asked to identify three EEG channels that were affected by the HFOV noise more strongly. What are the channels?  What is the magnitude of noise relative to signal in these channels?  Make sure to support this answer with measurements and plots from your results.

Solution

# Part 1: Simulated Signals

Question 1and 2 Figure 1: A signal sampled at 256Hz comprising 5 Hz, 30 Hz, 50 Hz, 64 Hz, and 70 Hz at amplitudes of 2mV, 1mV, 2mV, 4mV and 3mV respectively

Figure 1 shows a summation of waveforms sampled at 256Hz and comprising 5 Hz, 30 Hz, 50 Hz, 64 Hz, and 70 Hz at amplitudes of 2mV, 1mV, 2mV, 4mV and 3mV respectively. There are fast as well as slow variations of the wave pattern. The slow variation represents the low frequency component (5Hz) while the fastest variation represents the highest frequency components of the signal. In time domain, identifying the components of a signal can be difficult.

Code:

% Question 1

% generation of a sin signal consisting of 5Hz, 30Hz, 50Hz, 64Hz and 70Hz

f = [5 30 50 64 70];    % frequencies of the sine wave in Hz

% anonymous function for generation

g = @(f,fs,t, A, phi) sum(repmat(A’,1,length(t)).*sin(2*pi*repmat(f’,1,…

length(t)) .*repmat(t,length(f),1) + repmat(phi’,1,length(t))));

% specifications

fs = 256;                   % sampling frequency

phi = zeros(size(f));       % phase of each of the signals

A  = [2 1 2 4 3]*1e-3;      % signal amplitudes

N  = 1024;                  % number of points in a signal

T  = (N-1)/fs;              % duration of the signal – 1 period of least f

t  = 0:1/fs:T;              % time instants

s  = g(f, fs, t, A, phi);   % generated signal

% Question 2

% plot the signal

figure

plot(t, s*1e3, ‘LineWidth’, 1.5)

set(gca,’XLim’, [0 0.5])

xlabel(‘Time (s)’)

ylabel(‘Amplitude (mV)’)

grid

gridminor

Question 3 Figure 2: Frequency spectrum of the signal in Figure 1

In Figure 2, shows the magnitude spectrum of the signal in Figure 1. For accurate representation of a signal in the frequency spectrum, proper scaling is required to capture actual amplitudes (or more appropriately, magnitudes) of the components of the signal. Also, the component frequencies are accurately reflected. It is noted that there is a small bandwidth visible at each component frequency. This is just because of the number of points used to obtain the Discrete Fourier Transform. To increase this resolution, the number of points can be increased.

Evident in Figure 2 is how easy it is to identify the magnitude and frequency of the components of the signal in Figure 1. Therefore, to resolve a time domain waveform into its components, it is convenient to transform it to its frequency domain form.

Code:

% Question 3

% FFT of the signal using DSP System Toolbox

hfft = dsp.FFT(‘FFTLengthSource’, ‘Property’, ‘FFTLength’, N); % fft object

Y = step(hfft, s’);             % obtain FFT

Y = abs(Y(1:N/2+1)/N);          % magnitude

Y(2:end-1) = 2*Y(2:end-1);      % single sided spectrum

fvec = fs*(0:N/2)’/N;           % frequency points for plotting

% create plot

figure

plot(fvec, Y)

set(gca,’XLim’,[0 75])

xlabel(‘f (Hz)’)

ylabel(‘|Y|’)

title(‘Magnitude Spectrum’)

grid

gridminor

## Question 4and 5 Figure 3: A signal sampled at 128Hz comprising 5 Hz, 30 Hz, 50 Hz, 64 Hz, and 70 Hz at amplitudes of 2mV, 1mV, 2mV, 4mV and 3mV respectively

In Figure 3, a signal similar to what was sampled in Figure one is represented, only that the sampling rate is reduced to 128Hz. This limits the highest frequency component of the signal that can be accurately reproduced or represented by the sampled signal. This follows from the sampling theorem that states a signal can only be accurately reconstructed if its sampling rate is at least twice as high as the highest frequency component in the signal. This sampling rate is the Nyquist rate and the highest frequency component of the signal is the Nyquist frequency. The consequences of this are seen and explained in Figure 4.

Code:

% Question 4

% resampling at 128Hz

fs = 128;                   % new sampling frequency

T   = (N-1)/fs;             % duration of the signal – 1 period of least f

t   = 0:1/fs:T;             % time instants

s  = g(f, fs, t, A, phi);   % generated signal

% Question 5

% plot the signal

figure

plot(t, s*1e3, ‘LineWidth’, 1.5)

set(gca,’XLim’, [0 0.5])

xlabel(‘Time (s)’)

ylabel(‘Amplitude (mV)’)

grid

gridminor

## Question 6 Figure 4: Frequency spectrum of the signal in Figure 3

Figure 4 represents the magnitude spectrum of the signal in Figure 3, sampled at 128Hz. Therefore, the highest frequency that can be reconstructed from the samples is only 64Hz. The 70Hz component becomes a different frequency, within the range 0 to 64Hz, or more accurately, 58Hz (which is a reflection about half the sampling frequency). Since this frequency does not represent the actual frequency component in the signal, it is called an ‘alias’ of that component. Therefore, sampling at a rate lower than the Nyquist rate leads to aliasing. The magnitude of this component, as can be seen, is accurately captured.

Also, to sample the 64Hz component and reconstruct it, care has to be taken so that it is not sampled at its zero crossing, as is the case in this instance, because these are sine waves and the sampling begins at t=0s and then follows at intervals of 1/128s which are all zero crossings of the 64Hz component.  This is why the 64Hz component is missing from the magnitude spectrum of the wave. This problem is corrected in Figure 5 where the 64Hz component has a phase that allows for it to be sampled at locations other than the zero crossings.

Code:

% Question 6

% FFT of the signal using DSP System Toolbox

Y = step(hfft, s’);             % obtain FFT

Y = abs(Y(1:N/2+1)/N);          % magnitude

Y(2:end-1) = 2*Y(2:end-1);      % single sided spectrum

fvec = fs*(0:N/2)’/N;           % frequency points for plotting

% create plot

figure

plot(fvec, Y)

set(gca,’XLim’,[0 75])

xlabel(‘f (Hz)’)

ylabel(‘|Y|’)

title(‘Magnitude Spectrum’)

grid

gridminor

## Question 7 and 8 Figure 5:A signal sampled at 128Hz comprising 5 Hz, 30 Hz, 50 Hz, 64 Hz, and 70 Hz at amplitudes of 2mV, 1mV, 2mV, 4mV and 3mV respectively with 64Hz component having a phase of

In Figure 5, the same signal in Figure 1 is sampled at 128Hz as in Figure 3 but the 64Hz component has a phase of . This allows for the wave to be sampled at the peaks (both negative and positive) rather than the zero crossings. The spectrum of the signal therefore should contain the representation of the 64Hz component, unlike in Figure 4.

Code:

% Question 7

% sampling time of 128Hz but with a phase of pi/2 introduced for 64Hz

phi(f==64) = pi/2;              % phase of pi/2 at 64Hz

s = g(f, fs, t, A, phi);        % signal

% Question 8

% plot the signal

figure

plot(t, s*1e3, ‘LineWidth’, 1.5)

set(gca,’XLim’, [0 0.5])

xlabel(‘Time (s)’)

ylabel(‘Amplitude (mV)’)

grid

gridminor

## Question 9 Figure 6: Magnitude spectrum of the signal if Figure 5

In Figure 6, there is still aliasing caused by the 70Hz component at 58Hz. The 64Hz component is now represented because of the deliberate action of phase shifting by  done on the 64Hz component.

Code:

% Question 9

% FFT of the signal using DSP System Toolbox

Y = step(hfft, s’);             % obtain FFT

Y = abs(Y(1:N/2+1)/N);          % magnitude

Y(2:end-1) = 2*Y(2:end-1);      % single sided spectrum

fvec = fs*(0:N/2)’/N;           % frequency points for plotting

% create plot

figure

plot(fvec, Y)

set(gca,’XLim’,[0 75])

xlabel(‘f (Hz)’)

ylabel(‘|Y|’)

title(‘Magnitude Spectrum’)

grid

gridminor

The complete code for Part 1 is implemented in Part1.m.

# Part 2: Biomedical Data Analysis

A visual evaluation of the time domain signals shows that channels 5, 8, 9, 11, 17 through 21 are affected by HFOV.  Some, like channel 17, have a very high magnitude of the HFOV signal at 5Hz. The most affected channels are channel 5, channel 17 and channel 18. They HFOV signal magnitude is typically more than 10 times greater than in the channels with nearly no HFOV artifact. Figure 7: Channel 5 signal

Channel 5 signal represent in Figure 7 shows consistent higher frequency components superimposed on slower variations. These fast variations are the result of measurement noise due to HFOV signal. The spectral representation of the signal is shown in Figure 8 whereupon it is established that faster variations are due to a peak at 5Hz, with some visible harmonics at about 10Hz and 15Hz. Figure 8: Channel 5 signal magnitude spectrum

Figure 8 shows the spectral representation of channel 5 signal. The peak at about 5Hz represents the fast variations observed in the signal in Figure 7. It is the HFOV signal. Figure 9: Channel 14 signal

Figure 9 is a representation of a portion of channel 14 signal. As can be seen from the signal, the slow variations are prominent. The spectral response is shown in Figure 10. Figure 10: Channel 14 signal magnitude spectrum

Figure 10 shows the spectral response of channel 14 signal. From the response, the effect of HFOV signal is not pronounced. This is therefore one of the channels that can be used to represent EEG signal without the need for removal of HFOV noise. Figure 11: Channel 17 Signal

Figure 11 shows channel 17 signal and there is a very dominant high frequency variation. It is not easy to discern the EEG signal from this channel without further processing. The spectral representation of the signal is shown in Figure 12. Figure 12: Channel 17 signal magnitude spectrum

Figure 12 shows the spectral representation of the signal in Figure 11. As it can be seen, the HFOV component at about 5Hz has a magnitude larger than any other component of the signal. In Figure 8, even with the presence of the HFOV signal, it was not the largest and as such, slow variations were prominent in the signal as shown in Figure 7. This EEG signal cannot be used without further processing to remove the domineering HFOV component. Figure 13: Channel 18 signal

In Figure 13, the fast variations due to the HFOV measurement noise are clearly visible so much so that discerning the desired signal is somewhat difficult. The spectral representation of the signal is shown in Figure 14. Figure 14: Channel 18 signal magnitude spectrum

The spectral representation of channel 18 signal is shown in Figure 14. The HFOV component is shown to have its fundamental component at around 5Hz, with what seem to be the harmonics at around 10Hz and 15Hz. Again, this signal cannot be used to represent the patients EEG without further processing to remove the HFOV components.

Code:

The code for Part 2 is combined due to the repetitive nature of the involved processes.

% load EEG data and obtain other signal parameters

eeg_cols = 4:21;            % eeg columns

eeg = EEG_tbi(:,eeg_cols);  % eeg data

N = size(eeg,1);            % data points

fs = 200;                   % sampling frequency

t = (0:N-1)/fs;             % time instants

% create plots for each EEG signal

for n = 1:size(eeg,2)

% plot the signal

figure(‘Name’, sprintf(‘EEG Signal Channel %2d’, n+3), …

‘NumberTitle’, ‘off’)

plot(t, eeg(:,n), ‘LineWidth’, 1.5)

set(gca,’XLim’, [0 20])

xlabel(‘Time (s)’)

ylabel(‘Amplitude (\mu{V})’)

grid

gridminor

end

% create FFT plots

Nfft = 2^nextpow2(N); % FFT length

% fft object

hfft = dsp.FFT(‘FFTLengthSource’, ‘Property’, ‘FFTLength’, Nfft);

for n = 1:size(eeg, 2)

% FFT of the signal using DSP System Toolbox

s = eeg(:,n);                   % signal of interest

Y = step(hfft, s);             % obtain FFT

Y = abs(Y(1:Nfft/2+1)/Nfft);          % magnitude

Y(2:end-1) = 2*Y(2:end-1);      % single sided spectrum

fvec = fs*(0:Nfft/2)’/Nfft;           % frequency points for plotting

% create the magnitude plot

figure(‘Name’, sprintf(‘EEG Signal Magnitude Spectrum Channel%2d’, …

n+3), ‘NumberTitle’, ‘off’)

plot(fvec, Y)

set(gca,’XLim’,[0 20])

xlabel(‘f (Hz)’)

ylabel(‘|Y|’)

title(‘Magnitude Spectrum’)

grid

gridminor

end

From Figure 8, Figure 12, and Figure 14, the magnitude spectra of channel 5, channel 17 and channel 18 respectively; the HFOV artifact is visible at 5Hz. These are the three channels with the largest magnitude of the artifact. The effect of the HFOV artifact is not pronounce d in the signals from channel 14 as shown in Figure 10.

The complete code for Part II is implemented in the file Part2.m.

….

Part 1

Question 1:

The acceptable sampling frequency was 256Hz. To sample all the desired frequency components of a signal, the sampling frequency should be at least twice the maximum frequency component of the sampled signals. The frequency that is half the sampling rate is the maximum frequency that can be sampled without aliasing and it is called the Nyquist frequency and the corresponding sampling rate is called the Nyquist Rate. In the simulated signal, the Nyquist frequency is 70Hz and the Nyquist Rate is 140Hz.  Therefore, the sampling frequency of 128Hz is not acceptable. This rate is below the Nyquist rate.

Question 2:

Only 256Hz was an acceptable sampling frequency. 128Hz is below the sampling rate and produces an alias of the 70Hz signal in the frequency spectrum at 58Hz (a reflection of 70Hz about fs/2). If 140Hz is used, to avoid sampling the 70Hz signal at its points of zero crossing, the signal ought to be shifted in phase by  for the maxima to be sampled. To be guaranteed of sampling it regardless of the phase, the sampling frequency should be greater than 140Hz.

Question 3:

A sampling frequency of 128Hz results in an aliased spectrum. This is because the rate is lower than the Nyquist rate of 140Hz. The alias appears at 58Hz and it originates from 70Hz.

Question 4:

If the signal is band limited to contain only the 10Hz and the 25Hz components, a sampling of 50Hz is only appropriate if the 25Hz component has a phase shift of  so that the peaks are sampled. Without a phase shift, the signal is sampled at its zero crossing points and is not reflected at all in the frequency spectrum of the combined band-limited signal.

%Part 1 simulates signals

% prepare

clear all

close all

clc

% Question 1

% generation of a sin signal consisting of 5Hz, 30Hz, 50Hz, 64Hz and 70Hz

f = [5 30 50 64 70];    % frequencies of the sine wave in Hz

% anonymous function for generation

g = @(f,fs,t, A, phi) sum(repmat(A’,1,length(t)).*sin(2*pi*repmat(f’,1,…

length(t)) .*repmat(t,length(f),1) + repmat(phi’,1,length(t))));

% specifications

fs = 256;                   % sampling frequency

phi = zeros(size(f));       % phase of each of the signals

A  = [2 1 2 4 3]*1e-3;      % signal amplitudes

N  = 1024;                  % number of points in a signal

T  = (N-1)/fs;              % duration of the signal – 1 period of least f

t  = 0:1/fs:T;              % time instants

s  = g(f, fs, t, A, phi);   % generated signal

% Question 2

% plot the signal

figure

plot(t, s*1e3, ‘LineWidth’, 1.5)

set(gca,’XLim’, [0 0.5])

xlabel(‘Time (s)’)

ylabel(‘Amplitude (mV)’)

grid

grid minor

% Question 3

% FFT of the signal using DSP System Toolbox

hfft = dsp.FFT(‘FFTLengthSource’, ‘Property’, ‘FFTLength’, N); % fft object

Y = step(hfft, s’);             % obtain FFT

Y = abs(Y(1:N/2+1)/N);          % magnitude

Y(2:end-1) = 2*Y(2:end-1);      % single sided spectrum

fvec = fs*(0:N/2)’/N;           % frequency points for plotting

% create plot

figure

plot(fvec, Y)

set(gca,’XLim’,[0 75])

xlabel(‘f (Hz)’)

ylabel(‘|Y|’)

title(‘Magnitude Spectrum’)

grid

grid minor

% Question 4

% resampling at 128Hz

fs = 128;                   % new sampling frequency

T   = (N-1)/fs;             % duration of the signal – 1 period of least f

t   = 0:1/fs:T;             % time instants

s  = g(f, fs, t, A, phi);   % generated signal

% Question 5

% plot the signal

figure

plot(t, s*1e3, ‘LineWidth’, 1.5)

set(gca,’XLim’, [0 0.5])

xlabel(‘Time (s)’)

ylabel(‘Amplitude (mV)’)

grid

grid minor

% Question 6

% FFT of the signal using DSP System Toolbox

Y = step(hfft, s’);             % obtain FFT

Y = abs(Y(1:N/2+1)/N);          % magnitude

Y(2:end-1) = 2*Y(2:end-1);      % single sided spectrum

fvec = fs*(0:N/2)’/N;           % frequency points for plotting

% create plot

figure

plot(fvec, Y)

set(gca,’XLim’,[0 75])

xlabel(‘f (Hz)’)

ylabel(‘|Y|’)

title(‘Magnitude Spectrum’)

grid

grid minor

% Question 7

% sampling time of 128Hz but with a phase of pi/2 introduced for 64Hz

phi(f==64) = pi/2;              % phase of pi/2 at 64Hz

s = g(f, fs, t, A, phi);        % signal

% Question 8

% plot the signal

figure

plot(t, s*1e3, ‘LineWidth’, 1.5)

set(gca,’XLim’, [0 0.5])

xlabel(‘Time (s)’)

ylabel(‘Amplitude (mV)’)

grid

grid minor

% Question 9

% FFT of the signal using DSP System Toolbox

Y = step(hfft, s’);             % obtain FFT

Y = abs(Y(1:N/2+1)/N);          % magnitude

Y(2:end-1) = 2*Y(2:end-1);      % single sided spectrum

fvec = fs*(0:N/2)’/N;           % frequency points for plotting

% create plot

figure

plot(fvec, Y)

set(gca,’XLim’,[0 75])

xlabel(‘f (Hz)’)

ylabel(‘|Y|’)

title(‘Magnitude Spectrum’)

grid

grid minor

# Part 2

Question 1:

From the magnitude spectra of the channels exhibiting HFOV artifacts, the noise frequency of the HFOV signal is identified at 5Hz. In some cases, like that of channel 20, there is also an odd harmonic at 15Hz.

Question 2:

Channel 5, Channel 17 and Channel 18 are the most affected. The signal magnitude of the HFOV artifact in these channels is at least 10 times greater than the it is in the channels that are least affected.

%PART2 processes EEG data

% prepare

clear all

close all

clc

% load EEG data and obtain other signal parameters

eeg_cols = 4:21;            % eeg columns

eeg = detrend(EEG_tbi(:,eeg_cols))*1e-6;  % eeg data

N = size(eeg,1);            % data points

fs = 200;                   % sampling frequency

t = (0:N-1)/fs;             % time instants

% create plots for each EEG signal

for n = 1:size(eeg,2)

% plot the signal

figure(‘Name’, sprintf(‘EEG Signal Channel %2d’, n+3), …

‘NumberTitle’, ‘off’)

plot(t, eeg(:,n)*1e6, ‘LineWidth’, 1.5)

set(gca,’XLim’, [0 20])

xlabel(‘Time (s)’)

ylabel(‘Amplitude (\mu{V})’)

grid

grid minor

end

% create FFT plots

Nfft = 2^nextpow2(N); % FFT length

% fft object

hfft = dsp.FFT(‘FFTLengthSource’, ‘Property’, ‘FFTLength’, Nfft);

for n = 1:size(eeg, 2)

% FFT of the signal using DSP System Toolbox

s = eeg(:,n);                   % signal of interest

Y = step(hfft, s);             % obtain FFT

Y = abs(Y(1:Nfft/2+1)/Nfft);          % magnitude

Y(2:end-1) = 2*Y(2:end-1);      % single sided spectrum

fvec = fs*(0:Nfft/2)’/Nfft;           % frequency points for plotting

% create the magnitude plot

figure(‘Name’, sprintf(‘EEG Signal Magnitude Spectrum Channel%2d’, …

n+3), ‘NumberTitle’, ‘off’)

plot(fvec, Y)

set(gca,’XLim’,[0 20])

xlabel(‘f (Hz)’)

ylabel(‘|Y|’)

title(‘Magnitude Spectrum’)

grid

grid minor

end