How to do 1D signal analysis?

I just need an example:

%% =========================
%  [2.1] FFT calculation & plotting
%  [2.5] Band power computation
%  [1.1] Numerical integration (trapz)
%  Custom PSD function (FFT-based periodogram)
% =========================
function [f, PSD] = myPSD(x, Fs)
x = x(:);                 % force column vector
x = x - mean(x);          % remove DC offset

N = length(x);
X = fft(x);

K = floor(N/2);           % one-sided spectrum length index
X = X(1:K+1);

% One-sided PSD estimate
PSD = (1/(Fs*N)) * abs(X).^2;

% Double interior bins (not DC / Nyquist)
if K > 1
    PSD(2:end-1) = 2*PSD(2:end-1);
end

% Frequency vector (matches PSD length exactly)
f = (0:K)' * (Fs/N);
end

%% Compute PSDs for original signals
[f1, PSD1] = myPSD(sig1, Fs);
[f2, PSD2] = myPSD(sig2, Fs);

%% Plot PSDs (same axes is what the lab asks)
% NOTE: If plotting linear PSD, ylabel should NOT be dB/Hz.
figure;
plot(f1, PSD1, 'LineWidth', 1.1); hold on;
plot(f2, PSD2, 'LineWidth', 1.1);
grid on;
xlim([0 200]);
xlabel("Frequency (Hz)");
ylabel("Power/Frequency (linear)");
title("PSD of Signal 1 and Signal 2");
legend("sig1","sig2");

% Optional dB version (more readable visually)
figure;
plot(f1, pow2db(PSD1), 'LineWidth', 1.1); hold on;
plot(f2, pow2db(PSD2), 'LineWidth', 1.1);
grid on;
xlim([0 200]);
xlabel("Frequency (Hz)");
ylabel("Power/Frequency (dB/Hz)");
title("PSD of Signal 1 and Signal 2 (dB)");

%% 20-100 Hz Bandpower using trapz
% [2.5] Band power computation
% [1.1] Numerical integration (trapz)
idx1 = (f1 >= 20) & (f1 <= 100);
idx2 = (f2 >= 20) & (f2 <= 100);

bandpower1 = trapz(f1(idx1), PSD1(idx1));
bandpower2 = trapz(f2(idx2), PSD2(idx2));

fprintf('Original bandpower (20-100 Hz): sig1 = %.6g, sig2 = %.6g\n', bandpower1, bandpower2);

%% =========================
%  [1.5] Down-sampling
%  Anti-aliasing before downsampling (Butterworth lowpass)
% =========================
Fs_down = 200;                 % new sampling rate (Nyquist = 100 Hz)
downFactor = Fs / Fs_down;     % should be integer for downsample()

% Anti-alias lowpass filter (remove content above new Nyquist)
fc = 90;                       % use < 100 Hz for safety
order = 4;
Wn = fc / (Fs/2);              % normalized cutoff

[b, a] = butter(order, Wn, 'low');

sig1_filt = filtfilt(b, a, sig1);   % zero-phase filtering
sig2_filt = filtfilt(b, a, sig2);

% Downsample filtered signals
sig1_down = downsample(sig1_filt, downFactor);
sig2_down = downsample(sig2_filt, downFactor);

%% PSDs of anti-aliased + downsampled signals
[f1_down, PSD1_down] = myPSD(sig1_down, Fs_down);
[f2_down, PSD2_down] = myPSD(sig2_down, Fs_down);

figure;
plot(f1_down, PSD1_down, 'LineWidth', 1.1); hold on;
plot(f2_down, PSD2_down, 'LineWidth', 1.1);
grid on;
xlim([0 100]);
xlabel("Frequency (Hz)");
ylabel("Power/Frequency (linear)");
title("PSD of Downsampled Signals");
legend("sig1 down","sig2 down");

%% Bandpower (20-100 Hz) after anti-aliasing + downsampling
idx1_down = (f1_down >= 20) & (f1_down <= 100);
idx2_down = (f2_down >= 20) & (f2_down <= 100);

bandpower1_down = trapz(f1_down(idx1_down), PSD1_down(idx1_down));
bandpower2_down = trapz(f2_down(idx2_down), PSD2_down(idx2_down));

fprintf('Downsampled+AA bandpower (20-100 Hz): sig1 = %.6g, sig2 = %.6g\n', bandpower1_down, bandpower2_down);

%% =========================
%  [1.4] Signal rectification & enveloping
% =========================
% Lab asks: rectify, then lowpass at 2 Hz to create heartbeat envelope
sig1_rect = abs(sig1_down);
sig2_rect = abs(sig2_down);

fc_env = 2;                    % envelope lowpass cutoff
order_env = 4;
Wn_env = fc_env / (Fs_down/2);

[b_env, a_env] = butter(order_env, Wn_env, 'low');

env1 = filtfilt(b_env, a_env, sig1_rect);
env2 = filtfilt(b_env, a_env, sig2_rect);

t1 = (0:length(env1)-1)'/Fs_down;
t2 = (0:length(env2)-1)'/Fs_down;

figure;
plot(t1, env1); grid on;
xlabel("Time (s)"); ylabel("Envelope amplitude");
title("Envelope of Signal 1");

figure;
plot(t2, env2); grid on;
xlabel("Time (s)"); ylabel("Envelope amplitude");
title("Envelope of Signal 2");

%% Rescale envelopes to [0,1] (helps peak detection)
env1_scaled = rescale(env1);
env2_scaled = rescale(env2);

figure;
plot(t1, env1_scaled); grid on;
xlabel("Time (s)"); ylabel("Scaled envelope");
title("Scaled Envelope of Signal 1");

figure;
plot(t2, env2_scaled); grid on;
xlabel("Time (s)"); ylabel("Scaled envelope");
title("Scaled Envelope of Signal 2");

%% =========================
%  [1.2] Peak identification
% =========================
[pks1, locs1] = findpeaks(env1_scaled, 'MinPeakProminence', 0.1);
[pks2, locs2] = findpeaks(env2_scaled, 'MinPeakProminence', 0.1);

% Optional visualization of peaks
figure;
plot(t1, env1_scaled); hold on; grid on;
plot(locs1/Fs_down, pks1, 'o');
xlabel("Time (s)"); ylabel("Scaled envelope");
title("Signal 1 Peaks");

figure;
plot(t2, env2_scaled); hold on; grid on;
plot(locs2/Fs_down, pks2, 'o');
xlabel("Time (s)"); ylabel("Scaled envelope");
title("Signal 2 Peaks");

%% =========================
%  [1.2] Peak timing -> Heart rate + beat-to-beat timing
% =========================
tpeaks1 = locs1 / Fs_down;      % seconds
tpeaks2 = locs2 / Fs_down;

IBI1 = diff(tpeaks1);           % inter-beat intervals (s)
IBI2 = diff(tpeaks2);

HR1_bpm = 60 / mean(IBI1);      % average heart rate
HR2_bpm = 60 / mean(IBI2);

% Lab also asks for beat-to-beat timing variance
IBIvar1 = var(IBI1);
IBIvar2 = var(IBI2);

fprintf('HR sig1 = %.2f bpm, HR sig2 = %.2f bpm\n', HR1_bpm, HR2_bpm);
fprintf('Beat-to-beat variance: sig1 = %.6g s^2, sig2 = %.6g s^2\n', IBIvar1, IBIvar2);

%% =========================
%  [1.5] Interpolation (recreate original sig1 from downsampled)
%  Compare methods using RMSE
% =========================
t_down = (0:length(sig1_down)-1)' / Fs_down;
t_orig = (0:length(sig1)-1)' / Fs;

% Try a few methods (lab says experiment with interp1 methods)
methods = ["nearest","linear","spline","pchip","cubic"];
rmseVals = zeros(size(methods));

for k = 1:length(methods)
    sig1_int = interp1(t_down, sig1_down, t_orig, methods(k))';
    rmseVals(k) = sqrt(mean((sig1 - sig1_int).^2, 'omitnan'));
end

% Show RMSEs
disp(table(methods', rmseVals', 'VariableNames', {'Method','RMSE'}));

% Pick best method (lowest RMSE)
[bestRMSE, idxBest] = min(rmseVals);
bestMethod = methods(idxBest);

% Recompute best interpolation for plotting
sig1_best = interp1(t_down, sig1_down, t_orig, bestMethod)';

fprintf('Best interpolation method: %s (RMSE = %.6g)\n', bestMethod, bestRMSE);

% Plot original vs best interpolated (short window for visibility)
figure;
plot(t_orig, sig1, 'LineWidth', 1); hold on; grid on;
plot(t_orig, sig1_best, '--', 'LineWidth', 1);
xlim([0 2]);  % zoom into first 2 s so differences are visible
xlabel("Time (s)"); ylabel("Amplitude");
title("Original sig1 vs Interpolated sig1");
legend("Original","Interpolated");

[Mathieu NOE]

what is your question ?

[Joshua]

%% =========================================================

%% 1.1 Numerical integration & differentiation

%% =========================================================

% PRACTICE QUESTION:

% Given a time vector t (s) and signal x(t), compute:

% 1) dx/dt (numerical derivative)

% 2) the area under x(t) from t = 2 s to t = 5 s

% --- CODED ANSWER ---

% Assume t and x are already given

t = t(:); % force column vector (keeps dimensions consistent)

x = x(:);

dt = mean(diff(t)); % sample interval (assuming uniform sampling)

% Why: needed to convert "change per sample" into "change per second"

% 1) Numerical derivative (same length as x)

dxdt = gradient(x, dt);

% Why: derivative tells you how fast the signal is changing (slope/rate of change)

% Important for identifying rapid transitions, motion, onset changes, etc.

% 2) Numerical integration over a time window [2, 5] s

idx = (t >= 2) & (t <= 5);

% Why: selects only the time region we care about

area_2_5 = trapz(t(idx), x(idx));

% Why: integration adds up the signal over time (area under curve)

% Important for total quantity/accumulated signal/power-like calculations

% Optional plot

figure;

plot(t, x); hold on; grid on

plot(t(idx), x(idx), 'LineWidth', 1.5)

xlabel('Time (s)'); ylabel('x(t)');

title('1.1 Integration + Differentiation');

legend('Signal','Integrated Region');

%% =========================================================

%% 1.2 Peak & trough identification

%% =========================================================

% PRACTICE QUESTION:

% Given a 1D signal x(t), identify peaks and troughs using findpeaks.

% Return peak/trough amplitudes and locations, and plot them on the signal.

% --- CODED ANSWER ---

t = t(:);

x = x(:);

% Find peaks (tune parameters if needed)

[pks, locs_pk] = findpeaks(x, t, 'MinPeakProminence', 0.1);

% Why: finds meaningful local maxima (ignores tiny noisy bumps)

% Important for heart beats, pulse peaks, EMG bursts, etc.

% Find troughs by flipping the signal

[trs_neg, locs_tr] = findpeaks(-x, t, 'MinPeakProminence', 0.1);

trs = -trs_neg;

% Why: MATLAB finds peaks, not troughs, so we flip the signal to reuse findpeaks

% Important when you need minima timing/amplitude too

% Plot

figure;

plot(t, x); hold on; grid on

plot(locs_pk, pks, 'o', 'LineWidth', 1.5)

plot(locs_tr, trs, 'x', 'LineWidth', 1.5)

xlabel('Time (s)'); ylabel('Signal');

title('1.2 Peaks and Troughs');

legend('x(t)','Peaks','Troughs');

%% =========================================================

%% 1.3 Signal quality assessment

%% =========================================================

% PRACTICE QUESTION:

% Given a raw signal x_raw and a filtered version x_filt, estimate SNR (dB)

% assuming:

% signal = x_filt

% noise = demeaned raw - filtered

% --- CODED ANSWER ---

x_raw = x_raw(:);

x_filt = x_filt(:);

% Remove DC from raw signal first

x_dm = x_raw - mean(x_raw, 'omitnan');

% Why: DC offset is a constant shift, not useful signal content

% Important because offsets can distort power/SNR and filtering steps

% Estimate noise as difference between demeaned raw and filtered signal

noise_est = x_dm - x_filt;

% Why: what's removed by filtering is treated as "noise"

% Important for quantifying how much junk was in the raw signal

% Power-based SNR (works well for signal quality assessment)

P_signal = mean(x_filt.^2, 'omitnan');

P_noise = mean(noise_est.^2, 'omitnan');

% Why: power = average squared amplitude (standard way to compare signal vs noise)

SNR_dB = 10*log10(P_signal / P_noise);

% Why: expresses signal quality on a dB scale (easy to compare)

% Higher SNR = cleaner signal

fprintf('Estimated SNR = %.2f dB\n', SNR_dB);

% Optional plot

figure;

plot(x_dm); hold on; plot(x_filt); plot(noise_est);

grid on; title('1.3 Signal Quality Assessment');

legend('Demeaned Raw','Filtered (Signal)','Estimated Noise');

%% =========================================================

%% 1.4 Signal rectification & enveloping

%% =========================================================

% PRACTICE QUESTION:

% Given a filtered EMG/audio-type signal x and sampling frequency Fs:

% 1) Full-wave rectify the signal

% 2) Compute an envelope using a lowpass Butterworth filter

% 3) Plot raw/rectified/envelope

% --- CODED ANSWER ---

x = x(:);

% 1) Full-wave rectification

x_rect = abs(x);

% Why: makes all values positive so oscillations don't cancel out

% Important because EMG/audio oscillates around zero, but we care about "magnitude/activity"

% 2) Lowpass filter for envelope (example cutoff = 3 Hz)

fc_env = 3;

order_env = 4;

Wn_env = fc_env/(Fs/2);

% Why: lowpass keeps the slow trend and removes fast wiggles

% This smooth trend is the "envelope" (overall activation pattern)

[b_env, a_env] = butter(order_env, Wn_env, 'low');

x_env = filtfilt(b_env, a_env, x_rect);

% Why filtfilt: zero-phase filtering (doesn't shift peaks in time)

% Important for timing-based analysis (onsets, peaks, heart/EMG events)

% Time vector (if not provided)

t = (0:length(x)-1)'/Fs;

% 3) Plot

figure;

plot(t, x, 'DisplayName','Filtered signal'); hold on; grid on

plot(t, x_rect, 'DisplayName','Rectified');

plot(t, x_env, 'LineWidth',1.5, 'DisplayName','Envelope');

xlabel('Time (s)'); ylabel('Amplitude');

title('1.4 Rectification and Envelope');

legend;

%% =========================================================

%% 1.5 Interpolation & down-sampling

%% =========================================================

% PRACTICE QUESTION:

% Downsample a signal from Fs to 200 Hz (with anti-aliasing), then interpolate

% it back to the original time vector and compute RMSE.

% --- CODED ANSWER ---

x = x(:);

Fs_new = 200;

D = Fs / Fs_new;

% Why: downsample factor tells MATLAB how many samples to skip

% Important for reducing data size and processing load

% Anti-alias filter BEFORE downsampling (lowpass below new Nyquist)

fc_aa = 90;

[b_aa, a_aa] = butter(4, fc_aa/(Fs/2), 'low');

x_aa = filtfilt(b_aa, a_aa, x);

% Why: removes frequencies that would fold into lower frequencies (aliasing)

% Important because downsampling without this can corrupt the signal

% Downsample

x_down = downsample(x_aa, D);

% Why: actually reduces the sample rate by keeping every D-th sample

% Build time vectors

t_orig = (0:length(x)-1)'/Fs;

t_down = (0:length(x_down)-1)'/Fs_new;

% Interpolate back to original time points

x_interp = interp1(t_down, x_down, t_orig, 'cubic');

% Why: estimates missing points between downsampled samples

% Important for reconstructing/upsampling and comparing methods

% Compute reconstruction error

RMSE = sqrt(mean((x - x_interp).^2, 'omitnan'));

% Why: RMSE quantifies how close the reconstruction is to the original

% Lower RMSE = better interpolation

fprintf('Interpolation RMSE = %.6f\n', RMSE);

% Plot short segment

figure;

plot(t_orig, x); hold on; grid on

plot(t_orig, x_interp, '--');

xlim([0 2]);

xlabel('Time (s)'); ylabel('Amplitude');

title('1.5 Downsample + Interpolate');

legend('Original','Interpolated');

%% =========================================================

%% 2.1 FFT calculation & plotting

%% =========================================================

% PRACTICE QUESTION:

% Compute and plot the one-sided PSD of a signal x sampled at Fs.

% --- CODED ANSWER ---

x = x(:);

x = x - mean(x);

% Why: removes DC offset so the 0 Hz component doesn't dominate the spectrum

% Important for cleaner frequency analysis

N = length(x);

X = fft(x);

% Why: FFT converts time-domain signal into frequency-domain components

% Important for seeing what frequencies are present

K = floor(N/2);

X = X(1:K+1);

% Why: for real signals, the spectrum is symmetric, so we keep only one side

% Important for one-sided PSD plotting (0 to Nyquist)

% One-sided PSD estimate

PSD = (1/(Fs*N)) * abs(X).^2;

if K > 1

PSD(2:end-1) = 2*PSD(2:end-1);

end

% Why: converts FFT magnitude to power per Hz and accounts for removed negative frequencies

% Important for bandpower and frequency-domain comparisons

% Frequency axis (must match PSD length)

f = (0:K)' * (Fs/N);

% Why: maps each PSD bin to its actual frequency in Hz

% Plot in dB/Hz

figure;

plot(f, pow2db(PSD)); grid on

xlabel('Frequency (Hz)');

ylabel('Power/Frequency (dB/Hz)');

title('2.1 One-Sided PSD via FFT');

%% =========================================================

%% 2.2 FIR filtering

%% =========================================================

% PRACTICE QUESTION:

% Design an FIR lowpass filter with cutoff 40 Hz and apply it to x.

% Plot raw and filtered signals.

% --- CODED ANSWER ---

x = x(:);

% FIR = Finite Impulse Response (no feedback)

fir_order = 100;

fc = 40;

Wn = fc/(Fs/2);

% Why: normalized cutoff is required by MATLAB filter design functions

% Important because filter design always uses frequency relative to Nyquist

% Design FIR filter coefficients

b_fir = fir1(fir_order, Wn, 'low');

% Why: creates FIR filter coefficients for a lowpass filter

% Important for removing high-frequency noise while keeping lower frequencies

% Apply FIR filter

% filter(...) is causal; filtfilt(...) removes phase delay

x_fir = filtfilt(b_fir, 1, x);

% Why filtfilt: applies filter forward and backward to avoid time shift

% Important when timing/peak location matters

% Plot

t = (0:length(x)-1)'/Fs;

figure;

plot(t, x); hold on; grid on

plot(t, x_fir, 'LineWidth',1.2)

xlabel('Time (s)'); ylabel('Amplitude');

title('2.2 FIR Lowpass Filtering');

legend('Raw','FIR filtered');

%% =========================================================

%% 2.3 Notch filtering for mains interference removal

%% =========================================================

% PRACTICE QUESTION:

% Remove 60 Hz mains interference from a signal x sampled at Fs.

% --- CODED ANSWER ---

x = x(:);

% Mains interference is usually 60 Hz (in your course context)

f_notch = 60;

BW = 4;

% Why: notch is a narrow filter that removes one specific interference frequency

% Important because power-line noise often contaminates biosignals

% Normalize for iirnotch

wo = f_notch/(Fs/2);

Q = f_notch/BW;

bw = wo/Q;

% Why: iirnotch needs normalized center frequency and bandwidth

% Q controls how narrow/sharp the notch is

% Design and apply notch filter

[b_n, a_n] = iirnotch(wo, bw);

x_notch = filtfilt(b_n, a_n, x);

% Why: removes 60 Hz noise while keeping most of the rest of the signal

% Important for cleaner EMG/EEG and better PSD/features

% Plot (time-domain)

t = (0:length(x)-1)'/Fs;

figure;

plot(t, x); hold on; grid on

plot(t, x_notch);

xlabel('Time (s)'); ylabel('Amplitude');

title('2.3 60 Hz Notch Filtering');

legend('Raw','Notch filtered');

%% =========================================================

%% 2.4 Signal identification

%% =========================================================

% PRACTICE QUESTION:

% You are given a signal x and Fs. Use time and frequency plots to identify

% whether it likely contains low-frequency drift, mains noise (60 Hz), or

% high-frequency activity. (Return a brief interpretation.)

% --- CODED ANSWER ---

x = x(:);

t = (0:length(x)-1)'/Fs;

% Time-domain plot

figure;

plot(t, x); grid on

xlabel('Time (s)'); ylabel('Amplitude');

title('2.4 Signal Identification - Time Domain');

% Why: some artifacts are easier to see in time (drift, bursts, spikes)

% PSD (reuse FFT/PSD workflow)

x_dm = x - mean(x);

N = length(x_dm);

X = fft(x_dm);

K = floor(N/2);

X = X(1:K+1);

PSD = (1/(Fs*N))*abs(X).^2;

if K > 1

PSD(2:end-1) = 2*PSD(2:end-1);

end

f = (0:K)'*(Fs/N);

figure;

plot(f, pow2db(PSD)); grid on

xlabel('Frequency (Hz)'); ylabel('Power/Frequency (dB/Hz)');

title('2.4 Signal Identification - PSD');

% Why: frequency plot helps identify what kind of signal/artifact is present

% (e.g., 60 Hz noise, low-frequency drift, band-limited activity)

% Example "identification" checks (simple automatic flags)

% 1) Check for 60 Hz mains peak

idx60 = (f >= 58 & f <= 62);

has60Hz = any(PSD(idx60) > 5*median(PSD));

% Why: a strong narrow peak around 60 Hz often means mains interference

% 2) Check for low-frequency drift (< 2 Hz power)

idxLow = (f >= 0 & f <= 2);

idxMid = (f > 2 & f <= 20);

lowPower = trapz(f(idxLow), PSD(idxLow));

midPower = trapz(f(idxMid), PSD(idxMid));

hasDrift = lowPower > midPower;

% Why: if very low frequencies dominate, the signal likely has baseline drift

% 3) Simple printed interpretation

if has60Hz

disp('Likely mains interference near 60 Hz detected.');

end

if hasDrift

disp('Likely low-frequency drift/baseline wander present.');

end

if ~has60Hz && ~hasDrift

disp('Signal appears cleaner; inspect dominant bands manually.');

end

%% =========================================================

%% 2.5 Band power computation

%% =========================================================

% PRACTICE QUESTION:

% Given a PSD (f, PSD), compute the bandpower from 20 to 100 Hz.

% --- CODED ANSWER ---

% Assume f and PSD are already computed

f_low = 20;

f_high = 100;

idx_band = (f >= f_low) & (f <= f_high);

% Why: selects only the frequencies inside the band of interest

band_power_20_100 = trapz(f(idx_band), PSD(idx_band));

% Why: bandpower = area under the PSD in that band

% Important for comparing signal energy in specific frequency ranges

fprintf('Bandpower (20-100 Hz) = %.6f\n', band_power_20_100);

% Optional visualize selected band

figure;

plot(f, PSD); hold on; grid on

plot(f(idx_band), PSD(idx_band), 'LineWidth',1.5)

xlabel('Frequency (Hz)'); ylabel('Power/Frequency (linear)');

title('2.5 Band Power (20-100 Hz)');

legend('PSD','Integrated Band');

%% =========================================================

%% 3.1 2D image processing

%% =========================================================

% PRACTICE QUESTION:

% Load a noisy grayscale image, remove salt-and-pepper noise using medfilt2,

% compare to averaging filter, create a non-background mask, and adjust

% contrast using only non-background pixels.

% --- CODED ANSWER ---

% Load image

imOriginal = imread("Section3.png");

% Convert to grayscale if needed

if ndims(imOriginal) == 3

imOriginal = rgb2gray(imOriginal);

end

% Why: many image processing functions expect grayscale for simple intensity operations

% Display original

figure; imshow(imOriginal);

title('Original Image');

% 1) Median filter (good for salt-and-pepper noise)

imMed = medfilt2(imOriginal);

% Why: median filtering removes isolated bright/dark pixels without blurring edges too much

% Important for salt-and-pepper noise specifically

figure;

imshowpair(imOriginal, imMed, "montage", "scaling", "none");

title('Original (left) vs Median Filtered (right)');

% 2) Averaging filter (for comparison)

avgFilt = (1/9) * ones(3,3);

imAvg = filter2(avgFilt, imOriginal);

imAvg = uint8(imAvg);

% Why: averaging smooths by local mean, but it can blur edges/details

% Important to compare filter behavior and choose the right tool

figure;

imshowpair(imOriginal, imAvg, "montage", "scaling", "none");

title('Original (left) vs Averaging Filtered (right)');

% 3) Histogram of median-filtered image

figure;

imhist(imMed);

title('Histogram of Median-Filtered Image');

% Why: histogram shows how pixel intensities are distributed

% Important for deciding if contrast needs stretching

% 4) Create mask for non-background pixels (background assumed black = 0)

mask = (imMed ~= 0);

% Why: excludes black background so contrast adjustment is based on the tissue/image region

% Important because background can dominate histogram and reduce contrast stretch quality

% 5) Compute contrast limits using only non-background pixels

p = imMed(mask);

inLow = double(min(p))/255;

inHigh = double(max(p))/255;

% Why: get intensity limits from only useful pixels, then normalize to [0,1] for imadjust

% 6) Adjust contrast

imAdj = imadjust(imMed, [inLow inHigh], []);

% Why: stretches intensities to use more display range, improving visibility

% Important for clearer MRI detail

% 7) Histogram after adjustment

figure;

imhist(imAdj);

title('Histogram After Contrast Adjustment');

% 8) Final comparison

figure;

imshowpair(imOriginal, imAdj, "montage", "scaling", "none");

title('Original (left) vs Median + Contrast Adjusted (right)');

%% =========================================================

%% 3.2 Multichannel statistical analysis

%% =========================================================

% PRACTICE QUESTION:

% Given a CSV table with 10 subject columns of heart rate data (plus a first

% column that is time/index):

% 1) Compute mean, median, range, std for each subject

% 2) Find which subject has max and min mean HR

% 3) Compute hours each subject spent above (mean + 20 bpm)

% if each row = 5 seconds

% --- CODED ANSWER ---

% Load and convert table

T = readtable("Section2.csv");

A = table2array(T);

% Why: readtable is good for CSVs; table2array makes math/stats easier column-wise

% Select subject columns (assuming columns 2:11 are subjects)

X = A(:,2:11);

% Why: isolates just the 10 heart-rate signals (one column per subject)

% 1) Column-wise statistics (one value per subject)

mean_data = mean(X, 'omitnan');

median_data = median(X, 'omitnan');

range_data = max(X, [], 1) - min(X, [], 1);

std_data = std(X, 0, 'omitnan');

% Why: summarizes each subject's heart-rate behavior (central tendency + spread)

% Important for comparing subjects objectively

% 2) Find max/min mean HR and subject index

[max_mean, subj_max_mean] = max(mean_data);

[min_mean, subj_min_mean] = min(mean_data);

% Why: max/min gives both the value and which subject had it

% Important because lab asks for subject ID via coding (not by eye)

fprintf('Max mean HR: Subject %d (%.2f bpm)\n', subj_max_mean, max_mean);

fprintf('Min mean HR: Subject %d (%.2f bpm)\n', subj_min_mean, min_mean);

% 3) Time above (mean + 20 bpm), row spacing = 5 s

dt = 5;

hours_over = zeros(1,10);

for s = 1:10

threshold = mean_data(s) + 20;

count_above = sum(X(:,s) > threshold);

% Why: sum(logical) counts how many samples exceed the threshold

hours_over(s) = (count_above * dt) / 3600;

% Why: convert sample count -> seconds -> hours

% Important because lab asks for duration in hours, not sample count

end

% Display all subject results in a table

Subject = (1:10)';

resultsTbl = table(Subject, mean_data', median_data', range_data', std_data', hours_over', ...

'VariableNames', {'Subject','MeanHR','MedianHR','RangeHR','StdHR','HoursAboveMeanPlus20'});

% Why: makes results easy to read/check quickly

disp(resultsTbl);

[Image Analyst]

Yep - that's some code. But so what? Are we supposed to do anything about it? Do you have a specific question about it?

This looks like a homework problem. If you have any questions ask your instructor or read the link below to get started:

How do I get help on homework questions on MATLAB Answers? - MATLAB Answers - MATLAB Central

Obviously, if it is homework, we can't give you the full solution because you're not allowed to turn in our code as your own.

Another resource that MathWorks Central offers is the AI Chat Playground (on the blue banner above). You can paste your problem text into there and get code. However to be ethical you need to make sure that your professor is willing to accept AI generated code as your own. If he/she is against that and wants you to create your own code, and you submit code from AI or other people, you could be running a risk.

[Joshua]

How does this code look?

%% =========================================================

%% 1.1 Numerical integration & differentiation

%% =========================================================

% LIKELY TEST STYLE:

% Given a time vector t (s) and voltage signal V(t),

% 1) compute dV/dt in Volts/second

% 2) find the MAXIMUM slope

% 3) return the TIME where the maximum slope occurs

% 4) optionally compute area under the curve over a given time range

%

% KEY THINGS FROM FEEDBACK:

% - write dV/dt, not just dV

% - units must be Volts/second

% - return the TIME, not the sample index

% Assume t and V are already given

t = t(:);

V = V(:);

% Numerical derivative with respect to time

dVdt = gradient(V, t);

% gradient(V,t) directly uses the actual time vector

% so result is in Volts/second

% Maximum slope and the TIME where it occurs

[maxSlope, idxMaxSlope] = max(dVdt);

timeMaxSlope = t(idxMaxSlope);

fprintf('Maximum dV/dt = %.4f V/s\n', maxSlope);

fprintf('Time of maximum slope = %.4f s\n', timeMaxSlope);

% Example integration over a time window if asked

idx = (t >= 2) & (t <= 5);

area_2_5 = trapz(t(idx), V(idx));

fprintf('Area from 2 to 5 s = %.4f V*s\n', area_2_5);

% Optional plot

figure;

plot(t, V, 'b'); hold on; grid on

plot(t, dVdt, 'r')

plot(timeMaxSlope, maxSlope, 'ko', 'MarkerFaceColor','k')

xlabel('Time (s)');

ylabel('Amplitude');

title('1.1 Numerical Differentiation');

legend('V(t)','dV/dt','Max slope');

%% =========================================================

%% 1.2 Peak & trough identification

%% =========================================================

% LIKELY TEST STYLE:

% Given an ECG signal and time vector t,

% identify one consistent peak per heartbeat and calculate BPM.

%

% KEY THINGS FROM FEEDBACK:

% - ECG has multiple peaks in each beat, so plain findpeaks() is not enough

% - use extra arguments like MinPeakDistance / MinPeakHeight

% - do not confuse beats/sec with beats/min

% Assume t and ecg are already given

t = t(:);

ecg = ecg(:);

% Plot first to estimate reasonable parameters

figure;

plot(t, ecg); grid on

xlabel('Time (s)');

ylabel('ECG Amplitude');

title('ECG Signal');

% Example: detect R-peaks only

% Adjust values depending on signal appearance

[pks, locs] = findpeaks(ecg, t, ...

'MinPeakDistance', 0.5, ... % ~ prevents multiple peaks per beat

'MinPeakHeight', mean(ecg) + 0.3*std(ecg));

% Compute BPM

duration = t(end) - t(1);

beats = length(pks);

BPM = (beats / duration) * 60;

fprintf('Number of beats = %d\n', beats);

fprintf('Heart rate = %.2f BPM\n', BPM);

% Plot detected peaks

figure;

plot(t, ecg); hold on; grid on

plot(locs, pks, 'rv', 'MarkerFaceColor','r')

xlabel('Time (s)');

ylabel('ECG Amplitude');

title('1.2 Peak Detection for ECG');

legend('ECG','Detected R-peaks');

% If troughs are specifically asked for:

[trs_neg, locs_tr] = findpeaks(-ecg, t, ...

'MinPeakDistance', 0.5);

trs = -trs_neg;

%% =========================================================

%% 1.3 Signal quality assessment

%% =========================================================

% LIKELY TEST STYLE:

% Given a signal segment and a noise segment, compute SNR in dB.

%

% KEY THINGS FROM FEEDBACK:

% - read carefully what is defined as signal and what is defined as noise

% - know the SNR equation

%

% SNR(dB) = 10*log10(P_signal / P_noise)

% Assume signalSeg and noiseSeg are already given

signalSeg = signalSeg(:);

noiseSeg = noiseSeg(:);

% Compute power of signal and noise

P_signal = mean(signalSeg.^2, 'omitnan');

P_noise = mean(noiseSeg.^2, 'omitnan');

% Compute SNR in dB

SNR_dB = 10*log10(P_signal / P_noise);

fprintf('Signal power = %.6f\n', P_signal);

fprintf('Noise power = %.6f\n', P_noise);

fprintf('SNR = %.2f dB\n', SNR_dB);

% Optional comparison plot

figure;

subplot(2,1,1)

plot(signalSeg); grid on

title('Signal Segment')

subplot(2,1,2)

plot(noiseSeg); grid on

title('Noise Segment')

% IMPORTANT:

% Only use this if the question explicitly defines which part is signal

% and which part is noise. Do NOT assume filtered = signal unless the

% question tells you that.

%% =========================================================

%% 1.3 Signal quality assessment (Lab 2 style)

%% =========================================================

% PRACTICE QUESTION:

% You are given a raw EMG signal x (mV) sampled at Fs = 2000 Hz.

% Following the Lab 2 assumption:

% 1) Remove DC offset from x to get x_dm.

% 2) Create a filtered signal x_filt (bandpass + 60 Hz notch).

% 3) Estimate SNR (dB) assuming:

% signal = x_filt

% noise = x_dm - x_filt

% 4) Plot raw demeaned, filtered, and noise estimate.

%

% (This matches Lab 2 Part A step 8: assume filtered output is signal, and

% noise = unfiltered demeaned minus filtered.) [oai_citation:0‡Lab 2 - EMG Signal Processing & Control.pdf](sediment://file_00000000b1d871f5a702dd2b20ce2d1f)

% --- CODED ANSWER ---

% Assume x is provided. Fs is given by lab/test.

Fs = 2000;

x = x(:);

t = (0:length(x)-1)'/Fs; % time vector (s)

% Why: needed for time plots and selecting time windows if asked

% 1) Remove DC offset

x_dm = x - mean(x, 'omitnan');

% Why: DC offset is a constant shift (not muscle activity) and inflates power/SNR at 0 Hz

% 2) Filter to estimate the "true" EMG signal (bandpass + notch)

% Bandpass to keep typical EMG content and remove drift/high-frequency junk

bp_order = 4;

f_low = 20; % Hz (common EMG lower cutoff)

f_high = 450; % Hz (common EMG upper cutoff)

[b_bp,a_bp] = butter(bp_order, [f_low f_high]/(Fs/2), 'bandpass');

x_bp = filtfilt(b_bp, a_bp, x_dm);

% Why filtfilt: zero-phase filtering so timing/features aren't shifted

% Notch out mains interference at 60 Hz

f_notch = 60; % Hz

BW_notch = 4; % Hz (narrow removal band)

wo = f_notch/(Fs/2);

Q = f_notch/BW_notch;

bw = wo/Q;

[b_n,a_n] = iirnotch(wo, bw);

x_filt = filtfilt(b_n, a_n, x_bp);

% Why: removes narrow 60 Hz contamination that often appears in EMG recordings

% 3) Estimate noise and compute SNR

noise_est = x_dm - x_filt;

% Why (lab assumption): what gets removed by filtering is treated as "noise"

P_signal = mean(x_filt.^2, 'omitnan');

P_noise = mean(noise_est.^2, 'omitnan');

% Why: average power = mean(square), standard for SNR

SNR_dB = 10*log10(P_signal / P_noise);

% Why: dB scale is standard for comparing signal vs noise

fprintf('Estimated SNR = %.2f dB\n', SNR_dB);

% 4) Plot demeaned raw, filtered, and noise estimate (quality visualization)

figure; tiledlayout(1,3,'TileSpacing','compact','Padding','compact');

nexttile; plot(t, x_dm); grid on

title('Raw EMG (demeaned)'); ylabel('mV'); xlim([0 max(t)])

nexttile; plot(t, x_filt); grid on

title('Filtered EMG (signal)'); xlabel('Time (s)'); xlim([0 max(t)])

nexttile; plot(t, noise_est); grid on

title('Estimated noise (raw - filtered)'); xlim([0 max(t)])

%% =========================================================

%% 1.5 Interpolation & down-sampling

%% =========================================================

% LIKELY TEST STYLE:

% Given an original signal and time vector:

% 1) downsample the signal

% 2) build the correct downsampled time vector

% 3) interpolate back to the original time points using interp1

% 4) compute RMSE

%

% KEY THINGS FROM FEEDBACK:

% - interp1 inputs must be set up correctly

% - need BOTH original and downsampled time vectors

% - RMSE must be calculated correctly

% Assume t and x are already given

t = t(:);

x = x(:);

% Example downsampling: keep every 4th sample

D = 4;

x_down = x(1:D:end);

t_down = t(1:D:end);

% Interpolate back onto original time vector

x_interp = interp1(t_down, x_down, t, 'linear');

% Could also use 'spline' or 'cubic' if asked

% Remove any NaNs if interpolation does not cover exact endpoints

valid = ~isnan(x_interp);

% Compute RMSE

RMSE = sqrt(mean((x(valid) - x_interp(valid)).^2));

fprintf('RMSE = %.6f\n', RMSE);

% Plot comparison

figure;

plot(t, x, 'b'); hold on; grid on

plot(t_down, x_down, 'ro')

plot(t, x_interp, 'k--')

xlabel('Time (s)');

ylabel('Amplitude');

title('1.5 Downsampling and Interpolation');

legend('Original','Downsampled points','Interpolated');

% interp1 syntax to remember:

% interp1(x_known, y_known, x_query, method)

% Here:

% x_known = t_down

% y_known = x_down

% x_query = t

%% =========================================================

%% 2.3 Notch filtering for mains interference removal

%% =========================================================

% LIKELY TEST STYLE:

% Remove 60 Hz mains interference from a signal.

%

% KEY THINGS FROM FEEDBACK:

% - many people used too small a bandwidth

% - if noise is still present, increase bandwidth

% - always plot before and after to verify

% Assume x and Fs are already given

x = x(:);

% Plot raw signal spectrum first

N = length(x);

x_dm = x - mean(x);

X = fft(x_dm);

f = (0:floor(N/2))' * (Fs/N);

PSD = (1/(Fs*N)) * abs(X(1:floor(N/2)+1)).^2;

if length(PSD) > 2

PSD(2:end-1) = 2*PSD(2:end-1);

end

figure;

plot(f, PSD); grid on

xlabel('Frequency (Hz)');

ylabel('Power/Frequency');

title('Raw Signal Spectrum');

% Notch filter at 60 Hz

f0 = 60;

BW = 6; % increase bandwidth if 60 Hz still remains

wo = f0/(Fs/2);

bw = BW/(Fs/2);

[b,a] = iirnotch(wo, bw);

x_notch = filtfilt(b, a, x);

% Plot filtered spectrum to confirm 60 Hz removal

x_notch_dm = x_notch - mean(x_notch);

Xn = fft(x_notch_dm);

PSDn = (1/(Fs*N)) * abs(Xn(1:floor(N/2)+1)).^2;

if length(PSDn) > 2

PSDn(2:end-1) = 2*PSDn(2:end-1);

end

figure;

plot(f, PSDn); grid on

xlabel('Frequency (Hz)');

ylabel('Power/Frequency');

title('Filtered Signal Spectrum');

% Time-domain comparison

t = (0:length(x)-1)'/Fs;

figure;

plot(t, x); hold on; grid on

plot(t, x_notch)

xlabel('Time (s)');

ylabel('Amplitude');

title('2.3 Notch Filtering');

legend('Raw','Notch filtered');

% MAIN IDEA:

% If a peak near 60 Hz is still visible after filtering,

% the bandwidth was probably too narrow.

%% =========================================================

%% 2.4 Signal identification

%% =========================================================

% LIKELY TEST STYLE:

% Identify whether a signal is ECG, EMG, or EEG using the time plot

% and frequency content.

%

% KEY THING FROM FEEDBACK:

% - the test signal was ECG because it had a strong peak around heart-rate frequency

% - that kind of strong periodic low-frequency peak would not match EEG or EMG

% Assume x and Fs are already given

x = x(:);

t = (0:length(x)-1)'/Fs;

% Time plot

figure;

plot(t, x); grid on

xlabel('Time (s)');

ylabel('Amplitude');

title('Signal in Time Domain');

% Frequency analysis

x_dm = x - mean(x);

N = length(x_dm);

X = fft(x_dm);

K = floor(N/2);

X = X(1:K+1);

PSD = (1/(Fs*N)) * abs(X).^2;

if K > 1

PSD(2:end-1) = 2*PSD(2:end-1);

end

f = (0:K)' * (Fs/N);

figure;

plot(f, PSD); grid on

xlabel('Frequency (Hz)');

ylabel('Power/Frequency');

title('Signal Spectrum');

% Find dominant frequency

[~, idxPeak] = max(PSD(2:end)); % ignore 0 Hz

idxPeak = idxPeak + 1;

domFreq = f(idxPeak);

fprintf('Dominant frequency = %.2f Hz\n', domFreq);

% Basic interpretation guide

if domFreq >= 0.8 && domFreq <= 2.5

disp('Likely ECG: strong periodic component near heart-rate frequency.');

elseif domFreq > 20

disp('Possibly EMG: more high-frequency content.');

else

disp('Could be EEG or another low-frequency signal; inspect overall pattern.');

end

% QUICK ID RULES:

% ECG:

% - repeating heartbeat pattern in time

% - strong narrow peak around ~1 to 2 Hz

%

% EMG:

% - noisier, broader high-frequency content

% - not one clean repeating peak at heart rate

%

% EEG:

% - lower amplitude, not obvious heartbeat-like periodicity

% - does not usually show one dominant heart-rate peak

Answer 1

[Image Analyst]

Right before

%% Compute PSDs for original signals

[f1, PSD1] = myPSD(sig1, Fs);

[f2, PSD2] = myPSD(sig2, Fs);

you're going to have to come up with a sig1, sig2, and Fs and then run the program. If you want you can just use rand() or a rand vector added to a sine wave vector to give a noisy sinne wave signal. Like

numSamplePoints = 1000;
period = 500;
signalAmplitude = 10;
noiseAmplitude = 3;
maxXValue = 2000;
x = linspace(1, maxXValue, numSamplePoints);
sig1 = signalAmplitude * sin(2 * pi * x / period) + noiseAmplitude * rand(1, numSamplePoints);
plot(x, sig1, 'b-');
grid on;
yline(0, 'LineWidth', 1);
% Now create signal 2 with a different amplitude and phase
signalAmplitude = 14;
sig2 = signalAmplitude * sin(2 * pi * (x - 40) / period) + noiseAmplitude * rand(1, numSamplePoints);
hold on;
plot(x, sig2, 'r-');
legend('sig1', 'y=0', 'sig2')

●

Experiment around with different Fs and other parameters and observe how the fft spectrum changes.